From 8e27bcb36bca0f3d488d0416d42d5deff2e472ca Mon Sep 17 00:00:00 2001 From: Hou Shengren Date: Mon, 30 Mar 2026 15:07:01 +0800 Subject: [PATCH 1/2] refactor: retire network debug and gpu paths --- rl_adn/network/grid.py | 69 ++---- rl_adn/network/utils.py | 478 ++++++---------------------------------- 2 files changed, 80 insertions(+), 467 deletions(-) diff --git a/rl_adn/network/grid.py b/rl_adn/network/grid.py index 53ee256..90ac370 100644 --- a/rl_adn/network/grid.py +++ b/rl_adn/network/grid.py @@ -16,13 +16,12 @@ pre_power_flow_sam_sequential, pre_power_flow_tensor, ) -from rl_adn.network.utils import GPUPowerFlow, generate_network +from rl_adn.network.utils import generate_network try: import psutil except ImportError: psutil = None -from tqdm import trange class GridTensor: @@ -40,7 +39,6 @@ class GridTensor: nodes_frame (pd.DataFrame): DataFrame containing node data. Default is None. lines_frame (pd.DataFrame): DataFrame containing line data. Default is None. numba (bool): Flag to enable or disable Numba JIT compilation. Default is True. - gpu_mode (bool): Flag to enable or disable GPU mode. Default is False. """ def __init__( @@ -56,7 +54,6 @@ def __init__( nodes_frame: pd.DataFrame = None, lines_frame: pd.DataFrame = None, numba=True, - gpu_mode=False, ): self.s_base = s_base @@ -96,14 +93,12 @@ def __init__( self._power_flow_sam_sequential_constant_power_only = None self.is_numba_enabled = False - self.is_gpu_enabled = False if np.all(self.alpha_P) and not np.any(self.alpha_Z) and not np.any(self.alpha_I): self.constant_power_only = True self.start_time_pre_pf_tensor_constant_power_only = perf_counter() - # TODO: Change to sparse inverse. - self._K_ = np.array(-inv(self.Ydd_sparse).todense()) # Reduced version of -B^-1 (Reduced version of _F_) TODO: check it exist .toarray() + self._K_ = -inv(self.Ydd_sparse).toarray() self._L_ = self._K_ @ self.Yds # Reduced version of _W_ self.end_time_pre_pf_tensor_constant_power_only = perf_counter() else: @@ -117,10 +112,6 @@ def __init__( self.disable_numba() self.is_numba_enabled = False - if gpu_mode: - self.gpu_solver = GPUPowerFlow() - self.is_gpu_enabled = True - def enable_numba(self): """ Disables Numba JIT compilation, reverting to standard Python execution. @@ -176,7 +167,6 @@ def reset_start(self): """ Resets the starting voltage values for power flow calculations to default flat start values. """ - # TODO self.v_0 = np.ones((self.nb - 1, 1), dtype="complex128") # Flat start #2D array def _set_number_of_threads(self, threads): @@ -238,18 +228,15 @@ def _make_y_bus(self) -> None: # build Ybus Ybus = Cf.T * Yf + Ct.T * Yt # Full Ybus - # Dense matrices - # TODO - # TODO: This takes a lot of memory. Check if I can save it always as sparse for all methods. self._Ybus = Ybus.toarray() self.Yss = csr_matrix(Ybus[sl[0] - 1, sl[0] - 1], shape=(len(sl), len(sl))).toarray() - self.Ysd = np.array(Ybus[0, 1:].toarray()) # TODO: Here assume the slack is the first one? + self.Ysd = np.array(Ybus[0, 1:].toarray()) self.Yds = self.Ysd.T - self.Ydd = np.array(Ybus[1:, 1:].toarray()) # TODO: This consumes a huge amount of memory + self.Ydd = np.array(Ybus[1:, 1:].toarray()) self._Ybus_sparse = Ybus self.Yss_sparse = csr_matrix(Ybus[sl[0] - 1, sl[0] - 1], shape=(len(sl), len(sl))) - self.Ysd_sparse = Ybus[0, 1:] # TODO: Here assume the slack is the first one? + self.Ysd_sparse = Ybus[0, 1:] self.Yds_sparse = csc_matrix(self.Ysd.T) self.Ydd_sparse = Ybus[1:, 1:] @@ -314,10 +301,6 @@ def _compute_chunks(self, DIMENSION_BOUND, n_nodes, n_steps): the reminder: 2500-2000=500). """ - # DIMENSION_BOUND = 500 * 5_000 - # n_nodes = 4999 - # n_steps = 3000 - TS_MAX = DIMENSION_BOUND // n_nodes if n_steps > TS_MAX: # Chunk it (quotient, reminder) = divmod(n_steps, TS_MAX) @@ -331,8 +314,6 @@ def _compute_chunks(self, DIMENSION_BOUND, n_nodes, n_steps): else: # The requested amount of TS is lower than the bound. So, everything is ok idx = [0, n_steps] - # print(idx) - return idx def _make_big_sparse_matrices(self, S_nom, Ydd_sparse, Yds_sparse): @@ -428,6 +409,7 @@ def run_pf( tolerance: float = 1e-6, algorithm: str = "tensor", sparse_solver: str = "scipy", + show_progress: bool = False, ): """ Run a power-flow solve for the provided active/reactive power inputs. @@ -468,9 +450,6 @@ def run_pf( pf_algorithm = self.run_pf_tensor elif algorithm == "hp-tensor": pf_algorithm = self.run_pf_tensor_hp_laurent - elif algorithm == "gpu-tensor": - pf_algorithm = self.run_pf_tensor - kwargs.update(compute="gpu") else: raise ValueError("Incorrect power flow algorithm selected") @@ -480,6 +459,7 @@ def run_pf( flat_start=flat_start, start_value=start_value, tolerance=tolerance, + show_progress=show_progress, **kwargs, ) @@ -499,17 +479,14 @@ def run_pf_tensor( iterations: int = 100, tolerance: float = 1e-6, flat_start: bool = True, - compute: str = "cpu", + show_progress: bool = False, ) -> dict: if (active_power is not None) and (reactive_power is not None): - print("ok") assert len(active_power.shape) == 2, "Array must be two dimensional." assert len(reactive_power.shape) == 2, "Array must be two dimensional." assert active_power.shape[1] == reactive_power.shape[1] == self.nb - 1, "All nodes must have power values." else: - # active_power = self.P_file[np.newaxis, : reactive_power = self.Q_file[np.newaxis, :] - # print('zhong') self.ts_n = active_power.shape[0] # Time steps to be simulated if flat_start: @@ -530,30 +507,19 @@ def run_pf_tensor( n_steps = S_nom.shape[0] n_nodes = S_nom.shape[1] - if compute == "cpu": - # print("CPU Solver selected") - self._power_flow_tensor_solver = self._power_flow_tensor_constant_power - elif compute == "gpu" and self.is_gpu_enabled is False: - warnings.warn("GPU library not found, falling back to CPU.") - self._power_flow_tensor_solver = self._power_flow_tensor_constant_power - elif compute == "gpu" and self.is_gpu_enabled is True: - # print("GPU Solver selected") - self._power_flow_tensor_solver = self.gpu_solver.power_flow_gpu - - if compute == "cpu": - DIMENSION_BOUND = 500 * 100_000 # 5_000 x 10_000 did work. Empirical value for my machine - else: - DIMENSION_BOUND = 500 * 125_000 # 5_000 x 15_000 did work. Empirical value for my machine + self._power_flow_tensor_solver = self._power_flow_tensor_constant_power + DIMENSION_BOUND = 500 * 100_000 idx = self._compute_chunks(DIMENSION_BOUND, n_nodes=n_nodes, n_steps=n_steps) n_chunks = len(idx) - 1 - t = trange(n_chunks, desc="Chunk", leave=False) - for ii in t: - t.set_description(f"Chunk: {ii + 1} of {n_chunks}", refresh=True) + chunk_iterator = range(n_chunks) + if show_progress and n_chunks > 1: + from tqdm import trange + chunk_iterator = trange(n_chunks, desc="Chunk", leave=False) + for ii in chunk_iterator: ts_chunk = idx[ii + 1] - idx[ii] # Size of the chunk - # TODO self.v_0 = np.ones((ts_chunk, self.nb - 1)) + 1j * np.zeros((ts_chunk, self.nb - 1)) # Flat start @@ -578,7 +544,6 @@ def run_pf_tensor( end_time_pf = perf_counter() else: - # raise ValueError("This should not be running") start_time_pre_pf = perf_counter() self._F_, self._W_ = self._pre_power_flow_tensor( flag_all_constant_impedance_is_zero=self.flag_all_constant_impedance_is_zero, @@ -683,10 +648,8 @@ def run_pf_sam_sequential( reactive_power = self.Q_file if flat_start: - # TODO self.v_0 = np.ones((self.nb - 1, 1), dtype="complex128") # 2D-Vector elif start_value is not None: - # TODO: Check the dimensions of the flat start self.v_0 = start_value # User's start value active_power_pu = active_power / self.s_base # Vector with all active power except slack @@ -714,7 +677,7 @@ def run_pf_sam_sequential( else: start_time_pre_pf = perf_counter() B_inv, C, S_nom = self._pre_power_flow_sam_sequential( - active_power, # TODO: Change the input to S_nom + active_power, reactive_power, s_base=self.s_base, alpha_Z=self.alpha_Z, diff --git a/rl_adn/network/utils.py b/rl_adn/network/utils.py index 7d5fac7..e348e97 100644 --- a/rl_adn/network/utils.py +++ b/rl_adn/network/utils.py @@ -1,13 +1,9 @@ -import os -import sys -from ctypes import CDLL, POINTER, byref, c_bool, c_double, c_int -from time import perf_counter +from __future__ import annotations import matplotlib.pyplot as plt import networkx as nx import numpy as np import pandas as pd -from numpy import ctypeslib def _require_pandapower(): @@ -19,444 +15,98 @@ def _require_pandapower(): return pp, pandapower_topology -def load_library(): - """ - Loads a shared library for GPU-based power flow calculations. - - Returns: - A reference to the loaded shared library function for tensor power flow calculations. - - Description: - This function loads a shared library (.so or .dll file) based on the operating system. It sets up - the necessary function signatures for calling the GPU-based power flow solver. It supports Linux and Windows, - with a specific path for the shared library file. MacOS is not currently supported. - """ - platform = sys.platform - my_functions = None - - if platform == "linux" or platform == "linux2": - # linux (Hard coded in the meantime) - so_file = r"/home/mauricio/PycharmProjects/gpu_tensorpf/shared_library_complex.so" - my_functions = CDLL(so_file) - - elif platform == "darwin": - raise NotImplementedError("MacOS is not currently supported.") - - elif platform == "win32": - # TODO: Fallback if the gpu library is not found. - # Windows. - # os.add_dll_directory(r"C:\Users\20175334\Documents\PycharmProjects\tensorpowerflow\experiments\dll") - # dynamic_library = r"C:\Users\20175334\Documents\PycharmProjects\tensorpowerflow\experiments\dll\shared_library_complex.dll" - - # Hardcode the directory for debugging purposes. - os.add_dll_directory(r"C:\Users\20175334\source\repos\gpu_windows\x64\Debug") - dynamic_library = r"C:\Users\20175334\source\repos\gpu_windows\x64\Debug\gpu_windows.dll" - - my_functions = CDLL(dynamic_library, winmode=0) - else: - raise ValueError("OS not recognized.") - - tensor_power_flow = my_functions.tensorPowerFlow - tensor_power_flow.restype = None - ctypes_dtype_complex = ctypeslib.ndpointer(np.complex64) - tensor_power_flow.argtypes = [ - POINTER(ctypes_dtype_complex), # Matrix S, dimensions: S(m x m) - POINTER(ctypes_dtype_complex), # Matrix K, dimensions: S(m x m) - POINTER(ctypes_dtype_complex), # Matrix V0, dimensions: V0(m x p)\ - POINTER(ctypes_dtype_complex), # Matrix W, dimensions: W(m x 1) - POINTER(c_int), # m - POINTER(c_int), # p - POINTER(c_double), # tolerance - POINTER(c_int), # iterations - POINTER(c_bool), # convergence - ] - - return tensor_power_flow - - -class GPUPowerFlow(object): - """ - A class to handle GPU-based power flow calculations. - - Attributes: - gpu_solver (function): A reference to the loaded shared library function for tensor power flow calculations. - - Description: - This class provides an interface to perform power flow calculations using a GPU-based solver. It includes - a method to execute the power flow calculations with the necessary data transformations and GPU function calls. - """ - - def __init__(self): - self.gpu_solver = load_library() - - def power_flow_gpu( - self, - K: np.ndarray, # (m x m) == (nodes-1 x nodes-1) - L: np.ndarray, # (p x 1) == (time_steps x 1) or just p - S: np.ndarray, # (p x m) == (time_steps x nodes) - v0: np.ndarray, # (p x m) == (time_steps x nodes) - ts: np.ndarray, - nb: int, - iterations: int = 100, - tolerance: float = None, - ): - """ - Executes the GPU-based power flow calculation. - - Args: - K (np.ndarray): Admittance matrix (m x m). - L (np.ndarray): Load matrix (p x 1). - S (np.ndarray): Power matrix (p x m). - v0 (np.ndarray): Initial voltage matrix (p x m). - ts (np.ndarray): Time steps. - nb (int): Number of buses. - iterations (int, optional): Maximum number of iterations. Defaults to 100. - tolerance (float, optional): Convergence tolerance. If None, a default value is used. - - Returns: - Tuple[np.ndarray, int]: The voltage solution matrix and the number of iterations taken. - - Description: - This method wraps the GPU-based power flow solver, handling data preparation and result extraction. - It performs power flow calculations for a given network configuration and operating conditions. - """ - - if tolerance is None: - tolerance_gpu = 1e-10 - else: - tolerance_gpu = tolerance**2 # Heuristic, this match to the CPU tolerance. - - # ====================================================================== - # Reshape/casting and making sure that the complex matrices are 32 bits. - # S_host = S.T.copy() - S_host = S.T - S_host = S_host.astype(np.complex64) - - # K_host = K.copy() - K_host = K - K_host = K_host.astype(np.complex64) - - # V0_host = V0.T.copy() - V0_host = v0.T - V0_host = V0_host.astype(np.complex64) - - # W_host = L.copy() - W_host = L - W_host = W_host.astype(np.complex64) - - m = int(nb - 1) - p = int(ts) - - tolerance_gpu = float(tolerance_gpu) - iterations = int(iterations) - convergence = c_bool() - - # ====================================================================== - # Pointers for the dynamic library function. - m_int = byref(c_int(m)) - p_int = byref(c_int(p)) - iterations_int = byref(c_int(iterations)) - tolerance_int = byref(c_double(tolerance_gpu)) - convergence_int = byref(convergence) - - ctypes_dtype_complex = ctypeslib.ndpointer(np.complex64) - - S_c = S_host.ravel(order="F") - S_ca = S_c.ctypes.data_as(POINTER(ctypes_dtype_complex)) - - K_c = K_host.ravel(order="F") - K_ca = K_c.ctypes.data_as(POINTER(ctypes_dtype_complex)) - - V0_c = V0_host.ravel(order="F") - V0_ca = V0_c.ctypes.data_as(POINTER(ctypes_dtype_complex)) - - W_c = W_host.ravel(order="F") - W_ca = W_c.ctypes.data_as(POINTER(ctypes_dtype_complex)) - - # start = perf_counter() - self.gpu_solver(S_ca, K_ca, V0_ca, W_ca, m_int, p_int, tolerance_int, iterations_int, convergence_int) - # print(f"GPU Dynamic library execution: {perf_counter() - start} sec.") - # print(f"Convergence: {convergence.value}") - v_solution = V0_c.reshape(m, p, order="F") - iter_solution = iterations_int._obj.value - - # Voltage solution is a matrix with dimensions (time_steps x (n_nodes-1))-> Including the transpose. - return v_solution.T, iter_solution - - -def generate_network(nodes, child=3, plot_graph=False, load_factor=2, line_factor=3): - """ - Generates a network graph and corresponding data frames for buses and lines. - - Args: - nodes (int): Number of nodes in the network. - child (int, optional): Number of children for each node in the tree. Defaults to 3. - plot_graph (bool, optional): Whether to plot the network graph. Defaults to False. - load_factor (int, optional): Factor to scale the load values. Defaults to 2. - line_factor (int, optional): Factor to scale the line impedance values. Defaults to 3. - - Returns: - Tuple[pd.DataFrame, pd.DataFrame]: Data frames for nodes and lines in the network. - - Description: - This function generates a network graph using a full binary tree and creates pandas data frames - for nodes and lines with their respective properties. It can also plot the network graph if required. - """ - LINES = nodes - 1 - G = nx.full_rary_tree(child, nodes) +def generate_network(nodes: int, child: int = 3, plot_graph: bool = False, load_factor: int = 2, line_factor: int = 3): + """Generate a synthetic radial feeder and return node/line DataFrames.""" + line_count = nodes - 1 + graph = nx.full_rary_tree(child, nodes) if plot_graph: fig, ax = plt.subplots(1, 1, figsize=(10, 10)) - nx.draw_kamada_kawai(G, node_size=100, with_labels=True, font_size="medium", ax=ax) - - assert nodes == len(G.nodes) - assert LINES == len(G.edges) - - # Generate a pandas dataframe - PCT, ICT, ZCT = 1, 0, 0 - Tb, Pct, Ict, Zct = 0, PCT, ICT, ZCT - nodes_ = pd.DataFrame(list(G.nodes), columns=["NODES"]) + 1 - - active_ns = np.random.normal(50 * load_factor, scale=50, size=nodes).round(3) - reactive_ns = (active_ns * 0.1).round(3) - - power = pd.DataFrame({"PD": active_ns, "QD": reactive_ns}) - nodes_properties_ = pd.DataFrame(np.tile([[Tb, Pct, Ict, Zct]], (nodes, 1)), columns=["Tb", "Pct", "Ict", "Zct"]) - nodes_properties = pd.concat([power, nodes_properties_], axis=1) - nodes_properties = nodes_properties.astype({"Tb": int, "PD": float, "QD": float, "Pct": int, "Ict": int, "Zct": int}) - nodes_properties = nodes_properties[["Tb", "PD", "QD", "Pct", "Ict", "Zct"]] - nodes_properties.loc[0] = 1, 0.0, 0.0, PCT, ICT, ZCT # Slack - nodes_frame = pd.concat([nodes_, nodes_properties], axis=1) - - # R, X = 0.3144, 0.054 - R, X = 0.3144 / line_factor, 0.054 / line_factor - lines = pd.DataFrame.from_records(list(G.edges), columns=["FROM", "TO"]) + 1 # Count starts from 1 - lines_properties = pd.DataFrame(np.tile([[R, X, 0, 1, 1]], (LINES, 1)), columns=["R", "X", "B", "STATUS", "TAP"]) - lines_properties = lines_properties.astype({"R": float, "X": float, "B": int, "STATUS": int, "TAP": int}) - lines_frame = pd.concat([lines, lines_properties], axis=1) - - return nodes_frame, lines_frame - + nx.draw_kamada_kawai(graph, node_size=100, with_labels=True, font_size="medium", ax=ax) + + assert nodes == len(graph.nodes) + assert line_count == len(graph.edges) + + pct, ict, zct = 1, 0, 0 + nodes_frame = pd.DataFrame(list(graph.nodes), columns=["NODES"]) + 1 + active_power = np.random.normal(50 * load_factor, scale=50, size=nodes).round(3) + reactive_power = (active_power * 0.1).round(3) + bus_properties = pd.DataFrame( + { + "Tb": np.full(nodes, 0, dtype=int), + "PD": active_power, + "QD": reactive_power, + "Pct": np.full(nodes, pct, dtype=int), + "Ict": np.full(nodes, ict, dtype=int), + "Zct": np.full(nodes, zct, dtype=int), + } + ) + bus_properties.loc[0] = [1, 0.0, 0.0, pct, ict, zct] -def create_pandapower_net(network_info: dict, branch_info: pd.DataFrame = None, bus_info: pd.DataFrame = None): - """ - Creates a pandapower network from given network information. + resistance, reactance = 0.3144 / line_factor, 0.054 / line_factor + lines = pd.DataFrame.from_records(list(graph.edges), columns=["FROM", "TO"]) + 1 + line_properties = pd.DataFrame( + np.tile([[resistance, reactance, 0, 1, 1]], (line_count, 1)), + columns=["R", "X", "B", "STATUS", "TAP"], + ).astype({"R": float, "X": float, "B": int, "STATUS": int, "TAP": int}) - Args: - network_info (dict): A dictionary containing network parameters and file paths for branch and bus information. + return pd.concat([nodes_frame, bus_properties], axis=1), pd.concat([lines, line_properties], axis=1) - Returns: - pandapower.network: The created pandapower network. - Description: - This function reads network information from provided CSV files and creates a pandapower network - with buses, lines, loads, and an external grid connection. It sets up the network for power flow analysis. - """ +def create_pandapower_net(network_info: dict, branch_info: pd.DataFrame | None = None, bus_info: pd.DataFrame | None = None): + """Create a pandapower network from packaged network metadata and optional DataFrames.""" vm_pu = network_info["vm_pu"] branch_info_file = network_info["branch_info_file"] bus_info_file = network_info["bus_info_file"] pp, _ = _require_pandapower() - if branch_info is None: - branch_info = pd.read_csv(branch_info_file, encoding="utf-8") - else: - branch_info = branch_info.copy(deep=True) - if bus_info is None: - bus_info = pd.read_csv(bus_info_file, encoding="utf-8") - else: - bus_info = bus_info.copy(deep=True) + branch_frame = pd.read_csv(branch_info_file, encoding="utf-8") if branch_info is None else branch_info.copy(deep=True) + bus_frame = pd.read_csv(bus_info_file, encoding="utf-8") if bus_info is None else bus_info.copy(deep=True) net = pp.create_empty_network() - # Add buses - bus_dict = {} - for bus_name in bus_info["NODES"]: - bus_dict[bus_name] = pp.create_bus(net, vn_kv=11.0, name=f"Bus {bus_name}") + bus_lookup = {bus_name: pp.create_bus(net, vn_kv=11.0, name=f"Bus {bus_name}") for bus_name in bus_frame["NODES"]} - # Slack - bus_slack = bus_info[bus_info["Tb"] == 1]["NODES"].values - assert len(bus_slack.shape) == 1 and bus_slack.shape[0] == 1, "Only one slack bus supported" - pp.create_ext_grid(net, bus=bus_dict[bus_slack.item()], vm_pu=vm_pu, name="Grid Connection") + slack_bus = bus_frame[bus_frame["Tb"] == 1]["NODES"].values + if len(slack_bus) != 1: + raise ValueError("Exactly one slack bus is required for pandapower conversion") + pp.create_ext_grid(net, bus=bus_lookup[slack_bus.item()], vm_pu=vm_pu, name="Grid Connection") - # Lines - active_branches = branch_info[branch_info["STATUS"].astype(float) != 0].reset_index(drop=True) - for i, (_idx, (from_bus, to_bus, res, x_react, b_susceptance)) in enumerate(active_branches[["FROM", "TO", "R", "X", "B"]].iterrows()): + active_branches = branch_frame[branch_frame["STATUS"].astype(float) != 0].reset_index(drop=True) + for line_index, (_row_index, values) in enumerate(active_branches[["FROM", "TO", "R", "X", "B"]].iterrows(), start=1): + from_bus, to_bus, resistance, reactance, susceptance = values pp.create_line_from_parameters( net, - from_bus=bus_dict[from_bus], - to_bus=bus_dict[to_bus], + from_bus=bus_lookup[from_bus], + to_bus=bus_lookup[to_bus], length_km=1, - r_ohm_per_km=res, - x_ohm_per_km=x_react, - c_nf_per_km=b_susceptance, + r_ohm_per_km=resistance, + x_ohm_per_km=reactance, + c_nf_per_km=susceptance, max_i_ka=10, - name=f"Line {i + 1}", + name=f"Line {line_index}", ) - # Loads: - for node in bus_info["NODES"]: - pp.create_load(net, bus=bus_dict[node], p_mw=0.02, q_mvar=0.0, name="Load") - # print(f"Create net time: {perf_counter() - start}") - + for node in bus_frame["NODES"]: + pp.create_load(net, bus=bus_lookup[node], p_mw=0.02, q_mvar=0.0, name="Load") return net -def plot_pandapower_net(net): - """ - Plots a pandapower network. - - Args: - net (pandapower.network): The pandapower network to be plotted. - - Description: - This function creates a plot of the given pandapower network, showing buses, loads, PV generations, - and lines. It uses networkx for graph representation and matplotlib for plotting. - """ +def plot_pandapower_net(net) -> None: + """Plot a pandapower network using the topology graph helper.""" _, pandapower_topology = _require_pandapower() - # Create a graph from the pandapower network - G = pandapower_topology.create_nxgraph(net, respect_switches=False) - - # Set node positions based on bus coordinates - pos = {bus: (net.bus_geodata.at[bus, "x"], net.bus_geodata.at[bus, "y"]) for bus in G.nodes} - - # Draw buses - buses = [bus for bus in G.nodes if net.bus.at[bus, "type"] == "b"] - nx.draw_networkx_nodes(G, pos, nodelist=buses, node_color="red", node_size=200, label="Buses") + graph = pandapower_topology.create_nxgraph(net, respect_switches=False) + pos = {bus: (net.bus_geodata.at[bus, "x"], net.bus_geodata.at[bus, "y"]) for bus in graph.nodes} - # Draw loads - loads = [bus for bus in G.nodes if net.bus.at[bus, "type"] == "l"] - nx.draw_networkx_nodes(G, pos, nodelist=loads, node_color="blue", node_size=200, label="Loads") + buses = [bus for bus in graph.nodes if net.bus.at[bus, "type"] == "b"] + loads = [bus for bus in graph.nodes if net.bus.at[bus, "type"] == "l"] + pv_generations = [bus for bus in graph.nodes if net.bus.at[bus, "type"] == "s"] - # Draw PV generations - pv_generations = [bus for bus in G.nodes if net.bus.at[bus, "type"] == "s"] - nx.draw_networkx_nodes(G, pos, nodelist=pv_generations, node_color="green", node_size=200, label="PV Generations") + nx.draw_networkx_nodes(graph, pos, nodelist=buses, node_color="red", node_size=200, label="Buses") + nx.draw_networkx_nodes(graph, pos, nodelist=loads, node_color="blue", node_size=200, label="Loads") + nx.draw_networkx_nodes(graph, pos, nodelist=pv_generations, node_color="green", node_size=200, label="PV Generations") + nx.draw_networkx_edges(graph, pos, width=1.0, alpha=0.5) + nx.draw_networkx_labels(graph, pos, labels={bus: str(bus).split(" ")[-1] for bus in graph.nodes}, font_size=8) - # Draw lines - nx.draw_networkx_edges(G, pos, width=1.0, alpha=0.5) - - # Add labels to the nodes - node_labels = {bus: bus.split(" ")[-1] for bus in G.nodes} - nx.draw_networkx_labels(G, pos, labels=node_labels, font_size=8) - - # Add a legend plt.legend() - - # Display the plot plt.axis("off") plt.show() - - -def net_test(net): - """ - Tests the pandapower network with different power flow algorithms. - - Args: - net (pandapower.network): The pandapower network to be tested. - - Returns: - bool: True if the test passes, False otherwise. - - Description: - This function runs power flow calculations on the provided pandapower network using different algorithms - (Newton-Raphson and BFSW). It compares the calculated voltages with a predefined solution to verify the correctness. - """ - pp, _ = _require_pandapower() - - v_solution = [ - 0.98965162 + 0.00180549j, - 0.98060256 + 0.00337785j, - 0.96828145 + 0.00704551j, - 0.95767051 + 0.01019764j, - 0.94765203 + 0.01316654j, - 0.94090964 + 0.01600068j, - 0.93719984 + 0.01754998j, - 0.93283877 + 0.01937559j, - 0.93073823 + 0.02026054j, - 0.9299309 + 0.02058985j, - 0.92968994 + 0.02068728j, - 0.98003142 + 0.00362498j, - 0.97950885 + 0.00385019j, - 0.97936712 + 0.00391065j, - 0.97935604 + 0.0039148j, - 0.93971131 + 0.01547898j, - 0.93309482 + 0.01739656j, - 0.92577912 + 0.01988823j, - 0.91988489 + 0.02188907j, - 0.91475251 + 0.02362566j, - 0.90888169 + 0.02596304j, - 0.90404908 + 0.02788248j, - 0.89950353 + 0.02968449j, - 0.89731375 + 0.03055177j, - 0.89647201 + 0.03088507j, - 0.89622055 + 0.03098473j, - 0.94032081 + 0.01625577j, - 0.93992817 + 0.01642583j, - 0.93973182 + 0.01651086j, - 0.9301316 + 0.02052908j, - 0.92952481 + 0.02079761j, - 0.92922137 + 0.02093188j, - 0.92912022 + 0.02097663j, - ] - v_solution = np.array(v_solution, dtype="complex128") - - for pf_algorithm in ["nr", "bfsw"]: - print(f"Testing: {pf_algorithm} - Algorithm") - start = perf_counter() - if pf_algorithm == "bfsw": - pp.runpp(net, algorithm=pf_algorithm, numba=False, v_debug=True, VERBOSE=False, tolerance_mva=1e-6) - # print(f"BFSW. Iterations: {net._ppc['iterations']}. PF time: {net._ppc['et']}") - elif pf_algorithm == "nr": - pp.runpp(net, algorithm=pf_algorithm, numba=False, v_debug=True, VERBOSE=False, tolerance_mva=1e-6) - v_real = net.res_bus["vm_pu"].values * np.cos(np.deg2rad(net.res_bus["va_degree"].values)) - v_img = net.res_bus["vm_pu"].values * np.sin(np.deg2rad(net.res_bus["va_degree"].values)) - v_result = v_real + 1j * v_img - print("here starts the complex printing") - print(v_result) - # print(f"NR. Iterations: {net._ppc['iterations']}. PF time: {net._ppc['et']}") - print(f"Total pf time: {perf_counter() - start}.") - - v_real = net.res_bus["vm_pu"].values * np.cos(np.deg2rad(net.res_bus["va_degree"].values)) - v_img = net.res_bus["vm_pu"].values * np.sin(np.deg2rad(net.res_bus["va_degree"].values)) - v_result = v_real + 1j * v_img - assert np.allclose(v_result[1:], v_solution) - print("Test OK.") - - return True - - -def test_create_pandapower_net(network_info=None): - """Test for ``create_pandapower_net``. - - Creates a network using the provided ``network_info`` or a default - 34 node test network. The function prints the configuration and a - summary of the resulting pandapower ``net`` object. - """ - _require_pandapower() - if network_info is None: - data_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "data_sources", "network_data", "node_34")) - network_info = { - "branch_info_file": os.path.join(data_dir, "Lines_34.csv"), - "bus_info_file": os.path.join(data_dir, "Nodes_34.csv"), - "vm_pu": 1.0, - "s_base": 1000, - } - - print("Network configuration:") - for key, value in network_info.items(): - print(f" {key}: {value}") - - net = create_pandapower_net(network_info) - print("\nCreated pandapower network:") - print(net) - print("\nBus data:") - print(net.bus) - print("\nLine data:") - print(net.line) - print("\nLoad data:") - print(net.load) - return net - - -def test_plot_pandapower_net(net): - """Simple wrapper to plot a pandapower network.""" - print("Plotting pandapower network...") - plot_pandapower_net(net) From ad953431cbc954ec06e8a0fd6c210b01475a10e4 Mon Sep 17 00:00:00 2001 From: Hou Shengren Date: Mon, 30 Mar 2026 15:34:31 +0800 Subject: [PATCH 2/2] refactor: normalize network kernel naming --- rl_adn/network/grid.py | 163 +++++++-------- rl_adn/network/numbarize.py | 391 +++++++++++------------------------- 2 files changed, 191 insertions(+), 363 deletions(-) diff --git a/rl_adn/network/grid.py b/rl_adn/network/grid.py index 90ac370..d1ef0d0 100644 --- a/rl_adn/network/grid.py +++ b/rl_adn/network/grid.py @@ -80,29 +80,27 @@ def __init__( self._make_y_bus() self._compute_alphas() self.v_0 = None - self._F_ = None - self._W_ = None + self.tensor_factor_matrix = None + self.tensor_bias_vector = None - # Placeholder for the methods that are pre-compiled with numba. - self._power_flow_tensor_constant_power = None - self._pre_power_flow_tensor = None - self._power_flow_tensor = None - - self._pre_power_flow_sam_sequential = None - self._power_flow_sam_sequential = None - self._power_flow_sam_sequential_constant_power_only = None + # Runtime-selected numerical kernels. + self._constant_power_tensor_solver = None + self._tensor_prefactor_builder = None + self._tensor_solver = None + self._sam_prefactor_builder = None + self._sam_solver = None + self._sam_constant_power_solver = None self.is_numba_enabled = False if np.all(self.alpha_P) and not np.any(self.alpha_Z) and not np.any(self.alpha_I): - self.constant_power_only = True - self.start_time_pre_pf_tensor_constant_power_only = perf_counter() - - self._K_ = -inv(self.Ydd_sparse).toarray() - self._L_ = self._K_ @ self.Yds # Reduced version of _W_ - self.end_time_pre_pf_tensor_constant_power_only = perf_counter() + self.uses_constant_power_model = True + self.constant_power_prefactor_start_time = perf_counter() + self.constant_power_kernel = -inv(self.Ydd_sparse).toarray() + self.constant_power_slack_projection = self.constant_power_kernel @ self.Yds + self.constant_power_prefactor_end_time = perf_counter() else: - self.constant_power_only = False + self.uses_constant_power_model = False if numba: self.enable_numba() @@ -114,27 +112,28 @@ def __init__( def enable_numba(self): """ - Disables Numba JIT compilation, reverting to standard Python execution. + Enable Numba-backed kernels for the available solver paths. """ parallel = True - self._power_flow_tensor_constant_power = power_flow_tensor_constant_power - self._pre_power_flow_tensor = njit(pre_power_flow_tensor, parallel=parallel) - self._power_flow_tensor = njit(power_flow_tensor, parallel=parallel) + self._constant_power_tensor_solver = power_flow_tensor_constant_power + self._tensor_prefactor_builder = njit(pre_power_flow_tensor, parallel=parallel) + self._tensor_solver = njit(power_flow_tensor, parallel=parallel) - self._pre_power_flow_sam_sequential = njit(pre_power_flow_sam_sequential, parallel=parallel) - self._power_flow_sam_sequential = njit(power_flow_sam_sequential, parallel=parallel) - self._power_flow_sam_sequential_constant_power_only = njit(power_flow_sam_sequential_constant_power_only, parallel=parallel) + self._sam_prefactor_builder = njit(pre_power_flow_sam_sequential, parallel=parallel) + self._sam_solver = njit(power_flow_sam_sequential, parallel=parallel) + self._sam_constant_power_solver = njit(power_flow_sam_sequential_constant_power_only, parallel=parallel) def disable_numba(self): + """Fall back to the pure-Python numerical kernels.""" - self._power_flow_tensor_constant_power = power_flow_tensor_constant_power - self._pre_power_flow_tensor = pre_power_flow_tensor - self._power_flow_tensor = power_flow_tensor + self._constant_power_tensor_solver = power_flow_tensor_constant_power + self._tensor_prefactor_builder = pre_power_flow_tensor + self._tensor_solver = power_flow_tensor - self._pre_power_flow_sam_sequential = pre_power_flow_sam_sequential - self._power_flow_sam_sequential = power_flow_sam_sequential - self._power_flow_sam_sequential_constant_power_only = power_flow_sam_sequential_constant_power_only + self._sam_prefactor_builder = pre_power_flow_sam_sequential + self._sam_solver = power_flow_sam_sequential + self._sam_constant_power_solver = power_flow_sam_sequential_constant_power_only @classmethod def generate_from_graph(cls, *, nodes=100, child=2, plot_graph=True, load_factor=2, line_factor=3, **kwargs): @@ -244,7 +243,10 @@ def _make_y_bus(self) -> None: def _compute_alphas(self): """ - Computes alpha values for different load types in the grid. Assume P-1 and Z,I=0 + Initialize ZIP load coefficients. + + RL-ADN currently operates on the constant-power path, so the coefficients are + set to `P=1, I=0, Z=0`. """ self.alpha_P = 1 self.alpha_I = 0 @@ -415,16 +417,8 @@ def run_pf( Run a power-flow solve for the provided active/reactive power inputs. The method accepts either batched tensors or a single-step vector and dispatches - to the selected solver implementation. - "time_algorithm": Total time algorithm. time_algorithm = time_pre_pf + time_pf - "iterations": Total number of iterations to converge. - - "convergence": Boolean indicating: True: Algorithm converged, False: it didn't. - "iterations_log": NOT USED. - "time_pre_pf_log": NOT USED. - "time_pf_log": NOT USED. - "convergence_log": NOT USED. - } + to the selected solver implementation. It returns a dictionary containing the + complex voltage solution, timing statistics, and convergence metadata. """ is_tensor = False @@ -507,7 +501,6 @@ def run_pf_tensor( n_steps = S_nom.shape[0] n_nodes = S_nom.shape[1] - self._power_flow_tensor_solver = self._power_flow_tensor_constant_power DIMENSION_BOUND = 500 * 100_000 idx = self._compute_chunks(DIMENSION_BOUND, n_nodes=n_nodes, n_steps=n_steps) @@ -525,19 +518,19 @@ def run_pf_tensor( S_chunk = S_nom[idx[ii] : idx[ii + 1]] - if self.constant_power_only: - start_time_pre_pf = self.start_time_pre_pf_tensor_constant_power_only + if self.uses_constant_power_model: + start_time_pre_pf = self.constant_power_prefactor_start_time # No pre-computing (Already done when creating the object) - end_time_pre_pf = self.end_time_pre_pf_tensor_constant_power_only + end_time_pre_pf = self.constant_power_prefactor_end_time start_time_pf = perf_counter() - self.v_0, t_iterations = self._power_flow_tensor_solver( - K=self._K_, - L=self._L_, - S=S_chunk, - v0=self.v_0, - ts=ts_chunk, - nb=self.nb, + self.v_0, t_iterations = self._constant_power_tensor_solver( + kernel_matrix=self.constant_power_kernel, + slack_vector=self.constant_power_slack_projection, + nominal_power=S_chunk, + voltage_guess=self.v_0, + time_steps=ts_chunk, + node_count=self.nb, iterations=iterations, tolerance=tolerance, ) @@ -545,28 +538,28 @@ def run_pf_tensor( else: start_time_pre_pf = perf_counter() - self._F_, self._W_ = self._pre_power_flow_tensor( - flag_all_constant_impedance_is_zero=self.flag_all_constant_impedance_is_zero, - flag_all_constant_current_is_zero=self.flag_all_constant_current_is_zero, - flag_all_constant_powers_are_ones=self.flag_all_constant_powers_are_ones, - ts_n=ts_chunk, - nb=self.nb, - S_nom=S_chunk, - alpha_Z=self.alpha_Z, - alpha_I=self.alpha_I, - alpha_P=self.alpha_P, - Yds=self.Yds, - Ydd=self.Ydd, + self.tensor_factor_matrix, self.tensor_bias_vector = self._tensor_prefactor_builder( + all_constant_impedance_zero=self.flag_all_constant_impedance_is_zero, + all_constant_current_zero=self.flag_all_constant_current_is_zero, + all_constant_power_one=self.flag_all_constant_powers_are_ones, + time_steps=ts_chunk, + node_count=self.nb, + nominal_power=S_chunk, + alpha_z=self.alpha_Z, + alpha_i=self.alpha_I, + alpha_p=self.alpha_P, + yds=self.Yds, + ydd=self.Ydd, ) end_time_pre_pf = perf_counter() start_time_pf = perf_counter() - self.v_0, t_iterations = self._power_flow_tensor( - _F_=self._F_, - _W_=self._W_, - v_0=self.v_0, - ts_n=ts_chunk, - nb=self.nb, + self.v_0, t_iterations = self._tensor_solver( + tensor_factor_matrix=self.tensor_factor_matrix, + tensor_bias_vector=self.tensor_bias_vector, + voltage_guess=self.v_0, + time_steps=ts_chunk, + node_count=self.nb, iterations=iterations, tolerance=tolerance, ) @@ -658,14 +651,14 @@ def run_pf_sam_sequential( -1, ) - if self.constant_power_only: + if self.uses_constant_power_model: start_time_pre_pf = perf_counter() # No precomputing, the minimum matrix multiplication is done in the initialization of the object. end_time_pre_pf = perf_counter() start_time_pf = perf_counter() - V, iteration = self._power_flow_sam_sequential_constant_power_only( - B_inv=-self._K_, + V, iteration = self._sam_constant_power_solver( + B_inv=-self.constant_power_kernel, C=self.Yds.flatten(), v_0=self.v_0, s_n=S_nom, @@ -676,25 +669,25 @@ def run_pf_sam_sequential( else: start_time_pre_pf = perf_counter() - B_inv, C, S_nom = self._pre_power_flow_sam_sequential( + B_inv, C, S_nom = self._sam_prefactor_builder( active_power, reactive_power, s_base=self.s_base, - alpha_Z=self.alpha_Z, - alpha_I=self.alpha_I, - Yds=self.Yds, - Ydd=self.Ydd, - nb=self.nb, + alpha_z=self.alpha_Z, + alpha_i=self.alpha_I, + yds=self.Yds, + ydd=self.Ydd, + node_count=self.nb, ) end_time_pre_pf = perf_counter() start_time_pf = perf_counter() - V, iteration = self._power_flow_sam_sequential( - B_inv, - C, - v_0=self.v_0, - s_n=S_nom, - alpha_P=self.alpha_P, + V, iteration = self._sam_solver( + inverse_matrix_b=B_inv, + matrix_c=C, + voltage_guess=self.v_0, + nominal_power=S_nom, + alpha_p=self.alpha_P, iterations=self.iterations, tolerance=self.tolerance, ) diff --git a/rl_adn/network/numbarize.py b/rl_adn/network/numbarize.py index 2446d70..f342c8b 100644 --- a/rl_adn/network/numbarize.py +++ b/rl_adn/network/numbarize.py @@ -1,332 +1,167 @@ +"""Numerical kernels used by the Laurent and SAM power-flow solvers.""" + +from __future__ import annotations + import numpy as np from numba import prange -def pre_power_flow_sam_sequential(active_power, reactive_power, s_base, alpha_Z, alpha_I, Yds, Ydd, nb): - """ - Prepares the matrices for the SAM sequential power flow method. - - Parameters: - active_power (np.ndarray): Array of active power values. - reactive_power (np.ndarray): Array of reactive power values. - s_base (float): Base power value for per-unit conversion. - alpha_Z (np.ndarray): Array of constant impedance values. - alpha_I (np.ndarray): Array of constant current values. - Yds (np.ndarray): Admittance matrix between slack and load buses. - Ydd (np.ndarray): Admittance matrix between load buses. - nb (int): Number of buses in the network. +def pre_power_flow_sam_sequential(active_power, reactive_power, s_base, alpha_z, alpha_i, yds, ydd, node_count): + """Precompute SAM sequential matrices for ZIP-style load models.""" + active_power_pu = active_power / s_base + reactive_power_pu = reactive_power / s_base + nominal_power = (active_power_pu + 1j * reactive_power_pu).reshape(-1) - Returns: - tuple: Tuple containing B_inv, C, and S_nom matrices used in the SAM sequential power flow method. - """ - active_power_pu = active_power / s_base # Vector with all active power except slack - reactive_power_pu = reactive_power / s_base # Vector with all reactive power except slack + if alpha_z.shape != nominal_power.shape: + raise ValueError("alpha_z must match the flattened nominal power shape") + if alpha_i.shape != nominal_power.shape: + raise ValueError("alpha_i must match the flattened nominal power shape") - S_nom = (active_power_pu + 1j * reactive_power_pu).reshape( - -1, - ) - if not np.any(alpha_Z): # \alpha_z is 0 - B_inv = np.linalg.inv(Ydd) + if not np.any(alpha_z): + inverse_matrix_b = np.linalg.inv(ydd) else: - # TODO: Assert that the shapoes of alpha_Z and S_nom are the same - B = np.diag(np.multiply(alpha_Z, np.conj(S_nom))) + Ydd - B_inv = np.linalg.inv(B) + matrix_b = np.diag(np.multiply(alpha_z, np.conj(nominal_power))) + ydd + inverse_matrix_b = np.linalg.inv(matrix_b) - if not np.any(alpha_I): # all \alpha_i are 0 - C = Yds # Constant + if not np.any(alpha_i): + matrix_c = yds else: - C = Yds + np.multiply(alpha_I, np.conj(S_nom)).reshape(nb - 1, 1) # Constant + matrix_c = yds + np.multiply(alpha_i, np.conj(nominal_power)).reshape(node_count - 1, 1) - return B_inv, C, S_nom + return inverse_matrix_b, matrix_c, nominal_power def power_flow_sam_sequential( - B_inv, - C, - v_0, - s_n, - alpha_P, + inverse_matrix_b, + matrix_c, + voltage_guess, + nominal_power, + alpha_p, iterations, tolerance, ): - """ - Performs the SAM sequential power flow calculation. - - Parameters: - B_inv (np.ndarray): Inverse of matrix B. - C (np.ndarray): Matrix C. - v_0 (np.ndarray): Initial voltage values. - s_n (np.ndarray): Power values. - alpha_P (np.ndarray): Array of constant power values. - iterations (int): Maximum number of iterations. - tolerance (float): Convergence tolerance. - - Returns: - tuple: Tuple containing the final voltage values and the number of iterations performed. - """ + """Run the sequential SAM fixed-point update for general ZIP loads.""" iteration = 0 - tol = np.inf - while (iteration < iterations) & (tol >= tolerance): - # Update matrices A and D: - A = np.diag(alpha_P * (1 / np.conj(v_0) ** 2) * np.conj(s_n)) - D = 2 * alpha_P * (1 / np.conj(v_0)) * np.conj(s_n) - - v = B_inv @ (A @ np.conj(v_0) - C - D) - tol = np.max(np.abs(np.abs(v) - np.abs(v_0))) - v_0 = v # Voltage at load buses + voltage_delta = np.inf + while (iteration < iterations) & (voltage_delta >= tolerance): + matrix_a = np.diag(alpha_p * np.reciprocal(np.conj(voltage_guess) ** 2) * np.conj(nominal_power)) + vector_d = 2 * alpha_p * np.reciprocal(np.conj(voltage_guess)) * np.conj(nominal_power) + + voltage_solution = inverse_matrix_b @ (matrix_a @ np.conj(voltage_guess) - matrix_c - vector_d) + voltage_delta = np.max(np.abs(np.abs(voltage_solution) - np.abs(voltage_guess))) + voltage_guess = voltage_solution iteration += 1 - return v_0, iteration # Solution of voltage in complex numbers + return voltage_guess, iteration def power_flow_sam_sequential_constant_power_only( - B_inv, - C, - v_0, - s_n, + inverse_matrix_b, + matrix_c, + voltage_guess, + nominal_power, iterations, tolerance, ): - """ - Performs the SAM sequential power flow calculation for constant power loads only. - - Parameters: - B_inv (np.ndarray): Inverse of matrix B. - C (np.ndarray): Matrix C. - v_0 (np.ndarray): Initial voltage values. - s_n (np.ndarray): Power values. - iterations (int): Maximum number of iterations. - tolerance (float): Convergence tolerance. - - Returns: - tuple: Tuple containing the final voltage values and the number of iterations performed. - """ + """Run the sequential SAM update for constant-power loads only.""" iteration = 0 - tol = np.inf - while (iteration < iterations) & (tol >= tolerance): - # Update matrices A and D: - A = np.diag((1 / np.conj(v_0) ** 2) * np.conj(s_n)) - D = 2 * (1 / np.conj(v_0)) * np.conj(s_n) - # D = D.reshape(-1, 1) - - v = B_inv @ (A @ np.conj(v_0) - C - D) - tol = np.max(np.abs(np.abs(v) - np.abs(v_0))) - v_0 = v # Voltage at load buses + voltage_delta = np.inf + while (iteration < iterations) & (voltage_delta >= tolerance): + matrix_a = np.diag(np.reciprocal(np.conj(voltage_guess) ** 2) * np.conj(nominal_power)) + vector_d = 2 * np.reciprocal(np.conj(voltage_guess)) * np.conj(nominal_power) + + voltage_solution = inverse_matrix_b @ (matrix_a @ np.conj(voltage_guess) - matrix_c - vector_d) + voltage_delta = np.max(np.abs(np.abs(voltage_solution) - np.abs(voltage_guess))) + voltage_guess = voltage_solution iteration += 1 - return v_0, iteration # Solution of voltage in complex numbers + return voltage_guess, iteration def pre_power_flow_tensor( - flag_all_constant_impedance_is_zero, - flag_all_constant_current_is_zero, - flag_all_constant_powers_are_ones, - ts_n, - nb, - S_nom, - alpha_Z, - alpha_I, - alpha_P, - Yds, - Ydd, + all_constant_impedance_zero, + all_constant_current_zero, + all_constant_power_one, + time_steps, + node_count, + nominal_power, + alpha_z, + alpha_i, + alpha_p, + yds, + ydd, ): - """ - Prepares the matrices for the tensor-based power flow method. - - Parameters: - flag_all_constant_impedance_is_zero (bool): Flag indicating if all constant impedances are zero. - flag_all_constant_current_is_zero (bool): Flag indicating if all constant currents are zero. - flag_all_constant_powers_are_ones (bool): Flag indicating if all constant powers are one. - ts_n (int): Number of time steps. - nb (int): Number of buses. - S_nom (np.ndarray): Nominal power values. - alpha_Z (np.ndarray): Array of constant impedance values. - alpha_I (np.ndarray): Array of constant current values. - alpha_P (np.ndarray): Array of constant power values. - Yds (np.ndarray): Admittance matrix between slack and load buses. - Ydd (np.ndarray): Admittance matrix between load buses. - - Returns: - tuple: Tuple containing matrices _F_2, _W_2 used in the tensor power flow method. - """ - - if not flag_all_constant_impedance_is_zero: - _alpha_z_power = np.multiply(np.conj(S_nom), alpha_Z) # (ts x nodes) + """Precompute tensor fixed-point factors for batched Laurent solves.""" + if not all_constant_impedance_zero: + alpha_z_power = np.multiply(np.conj(nominal_power), alpha_z) else: - _alpha_z_power = np.zeros((ts_n, nb - 1)) # (ts x nodes) + alpha_z_power = np.zeros((time_steps, node_count - 1)) - if not flag_all_constant_current_is_zero: - _alpha_i_power = np.multiply(np.conj(S_nom), alpha_I) # (ts x nodes) + if not all_constant_current_zero: + alpha_i_power = np.multiply(np.conj(nominal_power), alpha_i) else: - _alpha_i_power = np.zeros((ts_n, nb - 1)) # (ts x nodes) + alpha_i_power = np.zeros((time_steps, node_count - 1)) - if flag_all_constant_powers_are_ones: - _alpha_p_power = np.conj(S_nom) # (ts x nodes) + if all_constant_power_one: + alpha_p_power = np.conj(nominal_power) else: - _alpha_p_power = np.multiply(np.conj(S_nom), alpha_P) # (ts x nodes) - - _B_inv2 = np.zeros((ts_n, nb - 1, nb - 1), dtype="complex128") - _F_2 = np.zeros((ts_n, nb - 1, nb - 1), dtype="complex128") - _W_2 = np.zeros((ts_n, nb - 1), dtype="complex128") + alpha_p_power = np.multiply(np.conj(nominal_power), alpha_p) - _C2 = _alpha_i_power + Yds.reshape(-1) # (ts x nodes) Sum is broadcasted to all rows of _alpha_i_power + inverse_matrix_b = np.zeros((time_steps, node_count - 1, node_count - 1), dtype="complex128") + tensor_factor_matrix = np.zeros((time_steps, node_count - 1, node_count - 1), dtype="complex128") + tensor_bias_vector = np.zeros((time_steps, node_count - 1), dtype="complex128") + matrix_c = alpha_i_power + yds.reshape(-1) - for i in prange(ts_n): - _B_inv2[i] = np.linalg.inv(np.diag(_alpha_z_power[i]) + Ydd) - _F_2[i] = -_B_inv2[i] * _alpha_p_power[i].reshape(1, -1) # Broadcast multiplication - _W_2[i] = (-_B_inv2[i] @ _C2[i].reshape(-1, 1)).reshape(-1) + for index in prange(time_steps): + inverse_matrix_b[index] = np.linalg.inv(np.diag(alpha_z_power[index]) + ydd) + tensor_factor_matrix[index] = -inverse_matrix_b[index] * alpha_p_power[index].reshape(1, -1) + tensor_bias_vector[index] = (-inverse_matrix_b[index] @ matrix_c[index].reshape(-1, 1)).reshape(-1) - return _F_2, _W_2 + return tensor_factor_matrix, tensor_bias_vector def power_flow_tensor( - _F_, - _W_, - v_0, - ts_n, - nb, + tensor_factor_matrix, + tensor_bias_vector, + voltage_guess, + time_steps, + node_count, iterations, tolerance, ): - """ - Performs the tensor-based power flow calculation. - - Parameters: - _F_ (np.ndarray): Matrix F. - _W_ (np.ndarray): Matrix W. - v_0 (np.ndarray): Initial voltage values. - ts_n (int): Number of time steps. - nb (int): Number of buses. - iterations (int): Maximum number of iterations. - tolerance (float): Convergence tolerance. - - Returns: - tuple: Tuple containing the final voltage values and the number of iterations performed. - """ + """Run the batched Laurent tensor power-flow update.""" iteration = 0 - tol = np.inf - while (iteration < iterations) & (tol >= tolerance): - v_recp_conj = np.reciprocal(np.conj(v_0)) - RT2 = np.zeros((ts_n, nb - 1), dtype="complex128") - for i in prange(ts_n): # This is critical as it makes a lot of difference - RT2[i] = _F_[i] @ v_recp_conj[i] - v = _W_ + RT2 - tol = np.max(np.abs(np.abs(v) - np.abs(v_0))) - v_0 = v + voltage_delta = np.inf + while (iteration < iterations) & (voltage_delta >= tolerance): + reciprocal_voltage = np.reciprocal(np.conj(voltage_guess)) + residual_term = np.zeros((time_steps, node_count - 1), dtype="complex128") + for index in prange(time_steps): + residual_term[index] = tensor_factor_matrix[index] @ reciprocal_voltage[index] + voltage_solution = tensor_bias_vector + residual_term + voltage_delta = np.max(np.abs(np.abs(voltage_solution) - np.abs(voltage_guess))) + voltage_guess = voltage_solution iteration += 1 - return v_0, iteration + return voltage_guess, iteration -def power_flow_tensor_constant_power_numba_parallel_True(K, L, S, v0, ts, nb, iterations, tolerance): - """ - Performs the tensor-based power flow calculation with constant power loads, optimized for parallel execution. - - Parameters: - K (np.ndarray): Matrix K. - L (np.ndarray): Matrix L. - S (np.ndarray): Power values. - v0 (np.ndarray): Initial voltage values. - ts (int): Number of time steps. - nb (int): Number of buses. - iterations (int): Maximum number of iterations. - tolerance (float): Convergence tolerance. - - Returns: - tuple: Tuple containing the final voltage values and the number of iterations performed. - """ - +def power_flow_tensor_constant_power(kernel_matrix, slack_vector, nominal_power, voltage_guess, time_steps, node_count, iterations, tolerance): + """Run the batched Laurent update for constant-power loads.""" iteration = 0 - tol = np.inf - while (iteration < iterations) & (tol >= tolerance): - v = np.zeros((ts, nb - 1), dtype="complex128") # TODO: Test putting this outside of while loop - for i in prange(ts): - v[i] = (K @ (np.conj(S[i]) * (1 / np.conj(v0[i]))).reshape(-1, 1) + L).T - tol = np.max(np.abs(np.abs(v) - np.abs(v0))) - v0 = v # Voltage at load buses - iteration += 1 - - return v0, iteration - - -def power_flow_tensor_constant_power(K, L, S, v0, ts, nb, iterations, tolerance): - """ - Performs the tensor-based power flow calculation for constant power loads. - - Parameters: - K (np.ndarray): Matrix K. - L (np.ndarray): Matrix L. - S (np.ndarray): Power values. - v0 (np.ndarray): Initial voltage values. - ts (int): Number of time steps. - nb (int): Number of buses. - iterations (int): Maximum number of iterations. - tolerance (float): Convergence tolerance. - - Returns: - tuple: Tuple containing the final voltage values and the number of iterations performed. - """ - iteration = 0 - tol = np.inf - S = S.T - v0 = v0.T - - LAMBDA = np.zeros((nb - 1, ts)).astype(np.complex128) - Z = np.zeros((nb - 1, ts)).astype(np.complex128) - voltage_k = np.zeros((nb - 1, ts)).astype(np.complex128) - - while iteration < iterations and tol >= tolerance: - LAMBDA = np.conj(S * (1 / v0)) # Hadamard product ( (nb-1) x ts) - Z = K @ LAMBDA # Matrix ( (nb-1) x ts ) - voltage_k = Z + L # This is a broadcasted sum dim => ( (nb-1) x ts + (nb-1) x 1 => (nb-1) x ts ) - tol = np.max(np.abs(np.abs(voltage_k) - np.abs(v0))) - v0 = voltage_k - iteration += 1 - - S = S.T # Recover the original shape of the power - v0 = v0.T # Recover the original shape of the power - - return v0, iteration - - -def power_flow_tensor_constant_power_new(K, L, S, v0, ts, nb, iterations, tolerance): - """ - A new version of the tensor-based power flow calculation for constant power loads, supporting parallel execution. - - Parameters: - K (np.ndarray): Matrix K. - L (np.ndarray): Matrix L. - S (np.ndarray): Power values. - v0 (np.ndarray): Initial voltage values. - ts (int): Number of time steps. - nb (int): Number of buses. - iterations (int): Maximum number of iterations. - tolerance (float): Convergence tolerance. - - Returns: - tuple: Tuple containing the final voltage values and the number of iterations performed. - """ - iteration = 0 - tol = np.inf - # S = S.T - # v0 = v0.T - - LAMBDA = np.zeros((nb - 1, ts)).astype(np.complex128) - Z = np.zeros((nb - 1, ts)).astype(np.complex128) - voltage_k = np.zeros((nb - 1, ts)).astype(np.complex128) - - voltage_k = voltage_k.T - W = L.ravel() - - while iteration < iterations and tol >= tolerance: - LAMBDA = np.conj(S.T * (1 / v0.T)) # Hadamard product ( (nb-1) x ts) - Z = K @ LAMBDA # Matrix ( (nb-1) x ts ) - Z = Z.T - for j in prange(ts): # This is a brodcasted sum ( (nb-1) x ts + (nb-1) x 1 => (nb-1) x ts ) - voltage_k[j] = Z[j] + W - - tol = np.max(np.abs(np.abs(voltage_k) - np.abs(v0))) - v0 = voltage_k + voltage_delta = np.inf + transposed_power = nominal_power.T + transposed_voltage = voltage_guess.T + + current_injection = np.zeros((node_count - 1, time_steps), dtype="complex128") + voltage_drop = np.zeros((node_count - 1, time_steps), dtype="complex128") + voltage_solution = np.zeros((node_count - 1, time_steps), dtype="complex128") + + while iteration < iterations and voltage_delta >= tolerance: + current_injection = np.conj(transposed_power * np.reciprocal(transposed_voltage)) + voltage_drop = kernel_matrix @ current_injection + voltage_solution = voltage_drop + slack_vector + voltage_delta = np.max(np.abs(np.abs(voltage_solution) - np.abs(transposed_voltage))) + transposed_voltage = voltage_solution iteration += 1 - return v0, iteration + return transposed_voltage.T, iteration