A spatial multi-scale reservoir computing framework for power flow analysis in power grids

2.1. Traditional numerical methods

Numerical methods model power grid operational constraints and objective functions as mathematical equations, deriving optimal flow solutions via numerical iterations^[18]. Examples include the Newton-Raphson, Gauss-Seidel, and fast-decoupled power flow methods, etc.^[19]. Additionally, to improve the convergence, Niu et al.^[20] designed methods based on Taylor series expansion and improved Newton methods. Montoya et al.^[21] proposed backward/forward sweep methods based on graph theory and improved algorithms based on rectangular coordinates, improving calculation speed and accuracy under specific conditions. Tang et al.^[22] developed an Alternating Current optimal power flow calculation method based on the Quasi-Newton method. Despite these advancements, fundamentally, the non-linear and non-convex nature of the alternating current optimal power flow (AC-OPF) problem means they are computationally intensive, highly sensitive to initial values, and can easily get trapped in local optima, failing to find a global-like solution. This results in a lack of convergence and robustness, rendering them inadequate for the fast-timescale control of modern complex power systems. To address these convergence bottlenecks and sensitivity to initial conditions, heuristic search strategies were introduced to enhance global optimization capabilities.

2.2. Heuristic search methods

Researchers have also explored heuristic search methods for power flow calculation. These methods use novel search mechanisms to enhance exploration and exploitation, improving convergence and accuracy. They are typically used for multi-objective optimization problems difficult for traditional techniques^[23,24]. The genetic algorithm was among the first heuristic methods applied to power flow^[25]. Others include teaching–learning optimization, metaheuristic algorithms, and differential evolution algorithms^[26,27], etc., all of which are used to solve optimal power flow problems involving both continuous and discrete control variables. However, these methods present trade-offs. While often more robust in avoiding poor local optima compared to numerical methods, they are notoriously computationally expensive and time-consuming^[28]. Their performance is sensitive to the initial encoding of the problem and the specific design of the fitness function, including the penalty mechanisms for handling operational constraints. And they face persistent challenges in generalizing effectively to high-dimensional nonlinear systems. These computational bottlenecks have motivated the emergence of data-driven methods, which aim to learn complex mappings directly from historical data without explicit iterative optimization.

2.3. Data-driven methods

With the rapid advancement of deep learning technologies and the growing complexity of modern power grids, there has been increasing research interest in data-driven approaches. These methods are capable of automatically learning feature patterns directly from data, without the need for precise prior knowledge^[29,30]. Currently, the models applied to power flow are primarily based on deep neural networks (DNNs) or graph neural networks (GNNs), as well as physics-informed neural networks expanded on this basis^[31]. Hansen et al.^[32] developed an end-to-end framework based on line GNNs to balance the power flow in energy grids. Lopez-Garciaet al.^[3] developed a typed GNN with linear complexity that solves the power flow problem through an unsupervised method of minimizing the violation of physical laws. In addition, ensuring feasible solutions is crucial in power flow problems. Some studies verify or refine model outputs using additional optimizers or classifiers^[33], while others enhance accuracy and feasibility by incorporating physical information^[34]. However, these approaches often rely on complex penalty-based loss functions to address the limitations of general models such as DNNs or GNNs^[35], which can increase computational cost, complicate gradient calculation, and sometimes lead to parameter sensitivity or unstable outputs. These observations motivate our proposed SMsRC framework that can organically integrate physical laws with efficient learning paradigms to simultaneously ensure computational speed and solution feasibility.

To highlight the research gap, Table 1 compares existing methods with the proposed SMsRC framework. SMsRC addresses previous limitations by integrating physical constraints and functional priors, balancing efficiency, generalization, and physical consistency.

3.1. Reservoir computing (RC)

RC, as a machine learning paradigm that has gained attention in the field of complex systems in recent years, has its core idea in performing nonlinear transformations with the help of a randomly connected neural network (i.e., the reservoir), and then achieving the output of specific tasks by training the weights of the readout layer^[36].

Echo State Network (ESN), as a prominent implementation of RC, features a basic architecture comprising three core parts: the input layer, the reservoir, and the readout layer. Operations related to these three parts mainly include initializing the input matrix $$ {\mathbf W}^{\text {in }} \in \mathbb{R}^{N_{x} \times N_{r}} $$ and the reservoir adjacency matrix $$ {\mathbf W} \in \mathbb{R}^{N_{r} \times N_{r}} $$, updating the reservoir state, and training the readout layer weights $$ {\mathbf W}^{\text {out }} \in \mathbb{R}^{N_{y} \times N_{r}} $$. $$ N_{x} $$, $$ N_{r} $$, and $$ N_{y} $$ represent the input dimension, the number of reservoir nodes, and the output dimension, respectively. Traditionally, $$ {\mathbf W}^{\mathrm{in}} $$ and $$ {\mathbf W} $$ are fixed after initialization and are typically randomly sampled from a uniform distribution, followed by appropriate sparsification and scaling. In this work, the traditional initialization approach is improved through the construction of a functional matrix, which will be discussed in detail later. A key characteristic of ESN is the "Echo State Property", which ensures that the state of the reservoir is determined primarily by the input rather than by the initial conditions^[37]. This is typically achieved by setting the spectral radius $$ \rho $$ of $$ {\mathbf W} $$ to a value close to, but less than 1.

For the problem of the power flow calculation, which is a system of nonlinear equations, the training process of the ESN can be characterized using ^[38,39]:

where $$ {\mathbf r}(t) \in \mathbb{R}^{N_r} $$ is the reservoir state vector, $$ {\mathbf x} \in \mathbb{R}^{N_x} $$ is the input vector, and $$ \alpha \in(0, 1] $$ is the leakage rate. $$ t_{\text {max }}=T $$ is the predefined number of iterations in the reservoir, $$ {\mathbf W}^{\text {out }} $$ is the output weight matrix, and $$ {\mathbf y} \in \mathbb{R}^{N_{y}} $$ is the output vector. Assuming that there are $$ N $$ training samples, the state matrix $$ {\mathbf R}=\left[{\mathbf r}_{1}(T), \ldots, {\mathbf r}_{N}(T)\right]^{\top} $$ and the target matrix $$ {\mathbf Y}=\left[{\mathbf y}_{1}, \ldots, {\mathbf y}_{N}\right]^{\top} $$ can be constructed. $$ J\left({\mathbf W}^{\text {out }}\right) $$ is the objective function for the ridge regression, where $$ \beta $$ is the regularization parameter, and $$ \|\cdot\|_{F} $$ denotes the Frobenius norm.

3.2. Power flow

In power system operations, two primary power flow models are widely used: direct current optimal power flow (DC-OPF) and AC-OPF. Although the DC-OPF model is relatively simple to compute, its precision is compromised due to the variations in reactive power and voltage magnitude. In contrast, the AC-OPF model considers more comprehensive system characteristics, providing a more accurate description of the actual system′s operational state. Therefore, this study adopts AC-OPF as the foundational model.

Under steady-state conditions, buses in a power grid can be classified into three types: PQ buses, PV buses, and slack bus. Each type of bus is characterized by distinct known and unknown parameters. Specifically, PQ buses correspond to load nodes, with known active power $$ P $$ and reactive power $$ Q $$, and unknown voltage magnitude $$ V $$ and stage angle $$ \theta $$; PV buses correspond to generator nodes, with known active power $$ P $$ and voltage magnitude $$ V $$, and unknown reactive power $$ Q $$ and stage angle $$ \theta $$; slack bus is the sole balance node in the system, with known voltage magnitude $$ V $$ and stage angle $$ \theta $$ serving as the global system reference, and unknown active power $$ P $$ and reactive power $$ Q $$.

The AC-OPF problem seeks to minimize the operational cost of the system while adhering to various physical and operational constraints^[40]:

where $$ N_{b} $$ and $$ N_{g} $$ represent the numbers of buses and generators in the system, respectively. The functions $$ f_{P}^{i}\left(P_{i}^{g}\right) $$ and $$ f_{Q}^{i}\left(Q_{i}^{g}\right) $$ denote the cost functions associated with the active and reactive power of the generator $$ i $$, respectively. $$ \mathcal{E} $$ represents the set of all transmission lines. $$ \theta_{i j} $$ is the stage angle of the transmission line $$ (i, j) $$, where $$ \theta_{i j}=\theta_{i}-\theta_{j} $$. $$ G_{i j} $$ and $$ B_{i j} $$ denote the conductance and susceptibility of the transmission line $$ (i, j) $$. $$ P_{i}^{g}\left(Q_{i}^{g}\right) $$ and $$ P_{i}^{d}\left(Q_{i}^{d}\right) $$ represent the generation and load values of active (reactive) power on the bus $$ i $$. Equations (6) and (7) represent the power balance constraints, reflecting Kirchhoff′s current law. Equations (8)-(11) specify the operational limits for generator active power, reactive power, bus voltage, and stage angle differences, ensuring the system operates within a safe and stable range.

3.3. Problem formulation

This study aims to develop a power flow solving framework based on the RC paradigm, integrating physical functional information of the power grid to enrich the priors. The problem can be formally defined as follows:

where $$ \mathcal{F} $$ is the power flow solution function to be designed, $$ {\mathbf X} \in \mathbb{R}^{N \times 2 N_{b}} $$ is the input matrix consisting of the known attribute values of all buses in the power grid, $$ \mathcal{P} $$ represents the additional physical information of the system, and $$ {\mathbf Y} \in \mathbb{R}^{N \times 2 N_{b}} $$ is the output matrix composed of the unknown attribute values. $$ N $$ represents the number of samples and $$ N_{b} $$ is the number of buses.

More specifically, each row of the input matrix $$ {\mathbf X} $$ corresponds to a system state and can be expressed as:

where $$ x_{i, 2 j-1} $$ and $$ x_{i, 2 j} $$ can represent the two known attribute values of the $$ j $$-th bus in the $$ i $$-th sample. For PQ buses, these values are $$ P $$ and $$ Q $$; for PV buses, they are $$ P $$ and $$ V $$; and for the slack bus, they are $$ V $$ and $$ \theta $$.

The structure of the output matrix $$ {\mathbf Y} $$ follows a similar structure but includes the unknown values for all buses:

where $$ y_{i, 2 j-1} $$ and $$ y_{i, 2 j} $$ can represent the two unknown attribute values of the $$ j $$-th bus in the $$ i $$-th sample. For PQ buses, these values are $$ \hat{V} $$ and $$ \hat{\theta} $$; for PV buses, they are $$ \hat{Q} $$ and $$ \hat{\theta} $$; and for slack bus, they are $$ \hat{P} $$ and $$ \hat{Q} $$.

4.1. Parallel learning auxiliary information

This stage focuses on learning two types of essential auxiliary information for the model: functional feature vectors and higher-order neighbors. The process consists of three steps. First, a parallel ESN architecture is used to generate bus-specific prediction outputs. Second, the resulting parameters are compressed into functional vectors that characterize the unique functional role of each bus. Finally, prediction losses between bus pairs are analyzed to identify higher-order neighbors-buses that exert significant influence despite not being directly connected.

4.1.1. Generating bus-specific prediction outputs via ESNs

First, the conventional ESNs are employed, where each of the $$ N_{b} $$ power grid buses corresponds to an independent reservoir (see Figure 1(b)). The weighted adjacency matrix $$ {\mathbf W} $$ of the reservoir is then reconstructed as a Watts-Strogatz network^[41] $$ {\mathbf W}^{\text{s}} $$, based on the original initialization. All reservoirs share a common input matrix $$ {\mathbf W}^{\text {in }} $$ and $$ {\mathbf W}^{\text{s}} $$, resulting in a total of $$ N_{b} $$ reservoirs $$ \mathcal{R}_{\mathrm{pl}}=\left\{{\mathbf r}_{0}, \ldots, {\mathbf r}_{N_{b}-1}\right\} $$. And a mask matrix $$ {\mathbf M}=\left[{\mathbf m}_{0}: \ldots: {\mathbf m}_{N_{b}-1}\right] $$ is introduced to filter the original input, ensuring that each reservoir only receives information relevant to its corresponding bus. Then $$ \mathcal{R}_{\mathrm{pl} } $$ undergoes parallel training. This approach leverages local information from a single bus to attempt to predict the global state of all buses. More precisely, the iterative update of the reservoir state $$ {\mathbf r}_{i} $$ for bus $$ i $$ is given by:

(16) $$\begin{array}{l} {\mathbf r}_{i} (t+1) \\ \;\;\; = \begin{cases} \begin{aligned}[t] & (1-\alpha) {\mathbf r}_{i}(t) \\ & \quad +\alpha \cdot \tanh \left({\mathbf W}^{\mathrm{in}}\left({\mathbf x} \odot {\mathbf m}_{i}\right)+{\mathbf W}^{\text{s}} {\mathbf r}_{i}(t)\right) , \end{aligned} & t=1 ; \\ (1-\alpha) {\mathbf r}_{i}(t)+\alpha \cdot \tanh ({\mathbf W} {\mathbf r}_{i}(t)), & t>1. \end{cases} \end{array}$$

The purpose here is to observe the global state from the perspective of a single local bus with a certain bias.

A set of readout layer weight matrices $$ \mathcal{W}_{\text{out}} = \{{\mathbf W}_{0}^{\text{out}}, \ldots, {\mathbf W}_{N_{b}-1}^{\text{out}}\} $$ can be obtained. $$ {\mathbf W}_{i}^{\text {out }} \in \mathbb{R}^{N_{y} \times N_{r}} $$ corresponding to each bus′s reservoir exhibits a sparse and inconsistent structure, demonstrating functional localization. As illustrated in Figure 2, the readout layer weights corresponding to buses 1, 2, 3, 8, 9, and 12 are displayed. One can observe that, while all reservoirs in the parallel architecture remain consistent, sharing $${\mathbf W}^{\mathrm{in}}$$ and $${\mathbf W}^{\mathrm{s}}$$, different buses activate distinct functional partitions within the reservoirs.

Figure 2. Schematic diagram of functional localization. Visualization of the readout layer weights of the parallel reservoirs corresponding to buses 1, 2, 3, 8, 9, and 12 in the IEEE14 bus system. For this test case, the output dimension is 28, and the number of nodes in each reservoir is 56.

4.2. Reconstructing input matrix and input information

This stage utilizes the two types of auxiliary information from the previous stage to construct two key components: the functional matrix $$ {\mathbf W}^{{\mathbf f}} $$ and the spatial multi-scale information $$ \mathcal{X}^{\text{SMs}} $$. First, the functional matrix is assembled to serve as a functionally-aware input mapping for the subsequent reservoirs. Second, physics-informed equivalent transformations are used to generate $$ \mathcal{X}^{\text{SMs}} $$, ensuring rational feature aggregation across bus types and enriching the model with prior knowledge.

4.3. Multi-scale RC with physically constrained readout layer

In the third stage, the spatial multi-scale RC (SMsRC) is constructed based on the functional matrix and spatial multi-scale information, as shown in Figure 1(e). To ensure the output solution conforms to actual power grid operation rules without needing an additional optimizer, a physical constraint module is integrated with the readout layer of the SMsRC. The previously improved RC, the efficient training of the RC readout layer, and the physical constraint module together enable the SMsRC framework to both ensure the feasibility of the output and maintain efficiency and precision.

First, the data of each scale in $$ \mathcal{X}^{\mathrm{SMs}} $$ are respectively mapped to the corresponding reservoirs using the same functional matrix $$ {{\mathbf W}^{\text{f}}}^{\prime} $$, thus forming three reservoirs denoted as $$ \mathcal{R}^{\mathrm{SMs}}=\left\{{\mathbf f}_{\mathrm{r} 1}, {\mathbf f}_{\mathrm{r} 2}, {\mathbf f}_{\mathrm{r} 3}\right\} $$, as illustrated in Figure 1(e). Then, these three reservoirs are connected to a single readout layer for an overall linear regression, which is not only more efficient but also compliant with the computational paradigm of the RC. At this time, the state matrix $$ {\mathbf R} $$ in Eq. (3) is replaced by the more information-rich $$ {\mathbf R}^{\text {SMs }} =\left[{\mathbf r}_{1}^{\mathrm{SMs}}(T), \ldots, {\mathbf r}_{N}^{\mathrm{SMs}}(T)\right]^\top $$, with

Now, the spatial multi-scale readout layer weights $$ {\mathbf W}^{\text {SMs}} $$ can be solved by using ridge regression, i.e.,

It should be noted that the values obtained are not the final results, but the intermediate variable $$ {\mathbf Y}_{\mathrm{vt}} $$, which represents the values of $$ \hat{V} $$ and $$ \hat{\theta} $$ for all buses. To ensure the feasibility of the solutions and their compliance with the power grid, additional physical constraints are applied at the readout layer. Specifically, based on the admittance matrix $$ {\mathbf P}_{\mathrm{y}} $$, the variables are reconstructed using the following power flow equations^[40]:

where $$ G $$ is the real part of the admittance matrix, $$ B $$ is the imaginary part, and this process also adheres to the constraints in Eqs. (6)-(11). Meanwhile, the physical constraint module maintains a low coupling with the overall framework. Finally, the unknown values are obtained for all buses, forming the feasible power flow solution $$ {\mathbf Y}_{\mathrm{pf}} $$.

To further assist the explanation and improve readability, Figure 3 provides a comprehensive flowchart that summarizes the key functional modules and their associated equations across Sections 4.1–4.3. This visual overview clarifies how each stage transforms input data through specific mathematical operations to produce the final feasible power flow solution.

Figure 3. Flowchart summarizing the key components and equations in Sections 4.1–4.3. Section 4.1 (Stage A): Parallel learning extracts functional vectors via reservoir state evolution (Eq. 16) and identifies higher-order neighbors via Bus-Pair loss (Eq. 17). Section 4.2 (Stage B): Constructs the functional matrix (Eq. 18) and generates spatial multi-scale input via physics-informed aggregation (Eqs. 19–21). Section 4.3 (Stage C): Trains multi-scale reservoir (Eq. 22) with ridge regression (Eq. 23) and applies physical power flow constraints (Eqs. 24 and 25) to ensure feasible solutions.

4.4. Algorithm pseudocode

The entire SMsRC can be summarized in Algorithm 1. The algorithm takes input sample data and physical information and outputs the corresponding power flow solution. In the first stage, by initializing the global input matrix $$ {\mathbf W}^{\text {in }} $$ and the reservoir adjacency matrix $$ {\mathbf W}^{\mathrm{s}} $$, the reservoir is constructed for each bus in the system and executed in parallel, and then the readout layer matrix and higher-order neighbors of each bus are identified. In the second stage, the input matrix is reconstructed using the functional matrix, and based on the number of known attribute values $$ k $$ for each bus, the spatial multi-scale input information $$ X^{\mathrm{SMs}} $$ is constructed according to Eqs. (19)-(21). In the third stage, the readout layer weights of the spatial multi-scale reservoir are trained, and the solution is reconstructed under physical constraints to obtain the final output.

5.1. Experimental setup

5.2. Main results

Figure 4 compares the MSE performance of different approaches across four power grids. Several conclusions can be drawn. First, the presented SMsRC framework outperforms other methods, particularly when the training set ratio is below about 50%. Even as the proportion of the training set increases, SMsRC maintains a leading position in most cases, with only a few instances where its performance is matched or slightly surpassed by PINN4PF or PIGN4PF. However, the standard deviation (STD) of SMsRC is smaller than that of PINN4PF and PIGN4PF, indicating better stability. Furthermore, as demonstrated in the next section, SMsRC also excels in computational efficiency compared to these methods.

Figure 4. Comparison of MSE performance across four power grids. (A) IEEE39 bus system. (B) IEEE57 bus system. (C) IEEE118 bus system. (D) IEEE300 bus system. Shaded areas indicate the Standard Deviation (STD). Insets show magnified results for training set ratios of 50%–90% for clearer presentation.

Second, under data scarcity conditions, RC paradigm shows significant accuracy advantages over non-RC-based methods. This phenomenon can be explained from both architectural and theoretical perspectives. Architecturally, RC only trains the linear readout layer while keeping the reservoir weights fixed, which fundamentally reduces the number of learnable parameters and thus the sample complexity required for reliable generalization^[52,53]. Theoretically, the reservoir serves as a high-dimensional nonlinear projection that maps input signals into a rich feature space where different system states become linearly separable^[39,54]. This "separation property" enables the simple linear readout to distinguish complex patterns even with limited training samples. Furthermore, the fading memory property (echo state property) of the reservoir ensures that the network state is primarily determined by recent input history rather than requiring extensive data to learn temporal dependencies^[37]. Recent studies have further demonstrated that RC can achieve effective predictions even from partial or incomplete observations^[55], supporting its robustness in data-scarce scenarios. Consequently, while deep learning methods typically require large datasets for backpropagation-based gradient descent to converge reliably, RC circumvents this bottleneck through its fixed reservoir structure and linear training paradigm.

Third, for RC-based methods, when the amount of data is already sufficient, simply increasing the training samples does little to improve the performance of these models. However, the accuracy of GraphRC and ParallelRC is unsatisfactory, particularly when the training set ratio is higher. This can be attributed to their simple focus on individual nodes and their surroundings, while neglecting higher-order and system-wide information. Although this approach may be sufficient for local-related tasks, power flow calculations require a comprehensive understanding of interactions among bus nodes and the system′s overall physical characteristics. Additionally, their incorporation of physical mechanisms is relatively limited.

5.3. Parameter complexity and stable efficiency

First, as summarized in Table 2, this work analyzes the parameter complexity of each method under identical experimental settings, considering factors such as the number of neurons per layer or module and the number of grid buses. The table explicitly presents the parameter count expressions for each method, where $$ n $$ represents the number of buses. It can be observed that, aside from ESN, the basic RC model, SMsRC exhibits superior parameter efficiency compared to the other models.

Then, by fixing the training set ratio at 70%, the comparison of different methods in terms of efficiency with stability and precision is visually presented in Figure 5. It is clear that RC-based methods offer advantages in computational efficiency (i.e., lower running time) and stability (i.e., more consistent MSE values), while physics-informed methods excel in average prediction accuracy. The presented SMsRC framework effectively balances both, particularly standing out as the top performer in terms of both MSE and computational speed on the larger-scale IEEE300 system.

Figure 5. Comparison of efficiency, stability, and accuracy with a fixed training set ratio of 70%. (A) IEEE39 bus system. (B) IEEE57 bus system. (C) IEEE118 bus system. (D) IEEE300 bus system. Each cluster comprises 100 independent experiments, plotting running time (x-axis) against MSE (y-axis) to demonstrate computational speed and solution consistency.

5.4. Ablation experiments

The ablation experiments are conducted on the proposed SMsRC framework, as shown in Figure 6. In these experiments, "w/o FM" denotes the ablation of the functional matrix, meaning the random input matrix initialization method is retained. "w/o MS" represents the ablation of the spatial multi-scale component, which effectively reduces the three reservoirs in Figure 1(e) to a single reservoir while maintaining the same total number of reservoir nodes, and it includes "w/o PI" and "w/o HB". "w/o PI" indicates that the framework does not incorporate physical information and relies solely on sample data, with spatial multi-scale information aggregation limited to a simple summation of neighboring attributes. Last, "w/o HB" signifies that the framework no longer uses higher-order neighbor information of buses, considering only topologically first-order neighbors.

Figure 6. Average MSE of various ablation methods across different datasets. (A) IEEE39 bus system. (B) IEEE57 bus system. (C) IEEE118 bus system. (D) IEEE300 bus system. The ablation settings are defined as: "w/o FM" (without functional matrix), "w/o MS" (without spatial multi-scale mechanism), "w/o PI" (without physical information integration), and "w/o HB" (without higher-order neighbor mining).

In terms of prediction accuracy, the functional matrix and spatial multi-scale are the most critical contributors to the framework, followed by higher-order neighbors and physical information. Additionally, spatial multi-scale and physical information can further enhance stability in data-scarce scenarios. This is evidenced by the steeper curves for the "w/o PI" and "w/o MS" ablation studies, signifying a less stable model transition from data scarcity to normal conditions, and the "w/o MS" setup typically performs worst when the amount of data is minimal.

From a theoretical standpoint, the ablation results validate that each proposed component contributes to enhanced data efficiency in distinct ways. The functional matrix, by encoding deterministic bus-specific roles rather than random weights, provides structured inductive biases that reduce the hypothesis space the model must explore, thereby requiring fewer samples to identify correct mappings^[38]. The spatial multi-scale mechanism aggregates information from higher-order neighbors, effectively increasing the information content available per sample by exploiting the inherent correlations in the power grid topology^[55]. Finally, the integration of physical information acts as domain-specific regularization, constraining the solution space to physically plausible outputs and preventing overfitting in low-data regimes^[56]. These three mechanisms synergistically enhance the framework′s sample efficiency, enabling robust performance even when only limited operational data is available.

5.5. Case analysis

Finally, a case study conducted on the 39-bus system in New England is presented in Figure 7. Specifically, Figure 7B-D shows the absolute errors between the model outputs and the true system values: (b) $$ \theta $$ values for all buses, (c) $$ Q $$ values for the PV buses, and (d) $$ V $$ values for the load (PQ) buses. All of them again address the great performance of the presented SMsRC approach. Figure 7A also highlights the higher-order neighbors of bus 17 identified by SMsRC, i.e., buses 11, 29, 33, and 38. Among them, bus 11 is a potentially influential node that may not be intuitively apparent, while buses 29, 33, and 38 are situated in the three generation areas closest to bus 17. Their strong influence on bus 17 aligns well with practical system behavior.

Figure 7. Case analysis of IEEE39 system located in New England. (A) The geographical distribution and connection relationships of 39 buses (Map data ©OpenStreetMap contributors). Green nodes are the direct neighbors of bus 17, and red nodes are the higher-order neighbors of bus 17 identified by SMsRC. (B-D) denote the absolute errors between the model outputs and the true system values for $$ \theta $$, $$ Q $$ and $$ V $$, respectively.

Authors′ contributions

Methodology, visualization, and writing-original draft: Zhang, H. F.; Zhang, Y. M.; Ding, X.; Ma, C.

Conceptualization, validation, writing-reviewing, and supervision: Xia, Y.; Tse, C. K.

Availability of data and materials

The raw data supporting the conclusions of this article will be made available by the authors upon request.

Financial support and sponsorship

This work is supported by the National Natural Science Foundation of China (62473001, 62476002), the Hong Kong Research Grants Council under GRF 112015/24, and the City University of Hong Kong (Grant No. 9229105).

Conflicts of interest

All authors declared that there are no conflicts of interest.

Ethical approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Copyright

© The Author(s) 2026.