Simplified A-Diakoptics for Accelerating QSTS Simulations

Abstract

The spread of distributed energy resources (DERs) across the distribution power system demands complex planning studies based on quasi-static-time-series (QSTS) simulations, requiring a significant amount of computing time to complete, leading planners to look for alternatives to QSTS. Diakoptics based on actors (A-Diakoptics) is a technique for accelerating simulations combining two computing techniques from different engineering fields: diakoptics and the actor model. Diakoptics is a mathematical method for tearing networks, reducing their complexity by using smaller subcircuits that can be solved independently. The actor model is used to coordinate the interaction between these subcircuits and their control actions, given the pervasive inconsistency that can be found when dealing with large-scale models. A-Diakoptics is a technique that simplifies the power flow problem for improving the simulation time performance, leading to faster QSTS simulations. This paper presents a simplified algorithm version of A-Diakoptics for modernizing sequential power simulation tools to use parallel processing. This simplification eliminates critical points found in previous versions of A-Diakoptics, improving the performance of the algorithm and facilitating its implementation to perform QSTS simulations. The performance of the new version of A-Diakoptics is evaluated by its integration into EPRI’s open-source simulator OpenDSS, which uses standard computing architectures and is publicly available.

1. Introduction

Quasi-static-time-series (QSTS) simulation is a valuable tool for assessing the behavior of distribution power systems over time [1]. By performing daily, yearly, and other time-based simulations, power engineers can characterize the impact of time-varying power devices such as solar photovoltaic generation (PV), energy storage, loads, generators, voltage regulators, and shunt capacitors, among other control devices within the distribution system.

The proliferation of distributed energy resources (DERs) has generated the need for studies such as hosting capacity, interconnection studies, microgrids, etc. These studies are QSTS based, and depending on the time-step resolution and simulation duration, the computational burden can be considerable, requiring a significant computing time to complete; this has led some planners and analysts to disregard QSTS simulations as an option for planning and other studies.

On the other hand, current computing architectures are branded by the introduction of multicore computing, a feature that allows computer applications to be accelerated by distributing tasks on multiple cores to work concurrently. This is conceptually straightforward, but it requires applications to be compatible with the parallel programming paradigm.

In the power system simulations domain, there are several techniques for making simulations compatible with parallel processing. Some of these techniques look for reducing the computational burden by dividing the total simulation time in multiple cores using temporal parallelization [2].

Other techniques link the computational burden to the circuit complexity, which is addressed by tearing the interconnected circuit into multiple subcircuits. These new subcircuits are solved separately using multicore computers to obtain the solution for the interconnected circuit. This approach is known as spatial parallelization, which addresses multiple techniques aiming to tear the interconnected model for its processing. These techniques are classified into two large groups, depending on the power flow problem formulation: bordered block diagonal matrix (BBDM) [3] and piecewise methods; the last one is the one in which diakoptics can be found, which is the base for the model presented in this paper called diakoptics based on actors (A-Diakoptics) [4,5]. A more detailed description and bibliography about the differences between the model-partitioning technique groups mentioned above can be found here [4,6].

A-Diakoptics combines two computing techniques from different engineering fields—diakoptics and the actor model. Diakoptics is a mathematical method for tearing networks [7,8]. Initially proposed by Gabriel Kron and later used and modified by other authors, diakoptics is a technique for tearing large physical circuits into several subcircuits to reduce the modeling complexity and accelerate the solution of the power flow problem using a computer network. In previous publications, the authors of this paper offered an extensive literature review on diakoptics as well as provided a populated list of references for the reader to review [4,5,6].

The actor model [9,10] is used to coordinate the interaction between subcircuits. The actor model is an information model for dealing with inconsistency robustness in parallel, concurrent, and asynchronous systems. This information model permits multiple processes to be executed in parallel sharing information using messages. As a result, information consistency can be guaranteed, avoiding common issues when working with parallel processing, such as race conditions and memory underutilization, among others.

A-Diakoptics is a technique that seeks to simplify the power flow problem to achieve a faster solution at each simulation step. Consequently, the total time reduction when performing QSTS will be evident at each simulation step with A-Diakoptics.

This paper presents a simplified algorithm for implementing A-Diakoptics in standard multicore computers. This simplification eliminates critical points found in previous versions of A-Diakoptics, improving the performance of the algorithm and facilitating its implementation to perform QSTS simulations. To evaluate the performance of the new version of A-Diakoptics, this method was integrated into EPRI’s open-source simulator OpenDSS, which can be executed in standard computing architectures and is publicly available.

2. Background

2.1. The Power Flow Problem in OpenDSS

While the power flow problem is probably the most common problem solved with the program, the OpenDSS is not best characterized as a power flow program. Its heritage is from general-purpose power system harmonics analysis tools. Thus, it works differently than most existing power flow tools. This heritage also gives it some unique and powerful capabilities for modeling complex electrical circuits. The program was originally designed to perform nearly all aspects of distribution planning for distributed generation (DG), which includes harmonics analysis. It is relatively easy to make a harmonics analysis program solve a power flow, while it can be quite difficult to make a power flow program perform harmonics analysis.

The OpenDSS program is designed to perform a basic distribution-style power flow in which the bulk power system is the dominant source of energy. However, it differs from the traditional radial circuit solvers in that it solves networked (meshed) distribution systems as easily as radial systems. It is intended to be used for distribution companies that may also have transmission or subtransmission systems. Therefore, it can also be used to solve small- to medium-sized networks with a transmission-style power flow.

Nearly all variables in the formulation result in a matrix or an array (vector) to represent a multiphase system. Many of the variables are complex numbers representing the common phasor notation used in frequency-domain AC power system analysis.

OpenDSS uses a standard nodal admittance formulation that can be found documented in many basic power system analysis texts. The textbook by Arrillaga and Watson [11] is useful for understanding this because it also develops the admittance models for harmonics analysis similarly to how OpenDSS is formulated. A primitive admittance matrix, Yprim, is computed for each circuit element in the model. These small matrices are used to construct the main system admittance matrix, Y_system, that knits the circuit model together. The solution is mainly focused on solving the nonlinear system admittance equation of the following form:

$$I_{PC}\left(E\right)=Y_{system}E$$

(1)

where I_PC(E) values are the compensation currents from power conversion (PC) elements in the circuit. The currents injected into the circuit from the PC elements, I_PC(E), are functions of voltage, as indicated, and represent the nonlinear portion of the currents from elements such as load, generator, PV system, and storage. There are several ways this set of nonlinear equations could be solved. The most popular way in OpenDSS is a simple fixed-point method that can be written concisely [12].

$$E_{n+1}={\left[Y_{system}\right]}^{-1}{I}_{PC}\left(E_n\right) n=0,1,2 \dots until converged$$

(2)

From (2), it can be inferred that every time, an expression such as the following is found:

$$E={\left[Y_{system}\right]}^{-1}I$$

(3)

The solver in OpenDSS is called an actor in OpenDSS since the introduction of the parallel processing suite [2], and this is the base for the A-Diakoptics analysis.

2.2. The Simplified A-Diakoptics Method

The initial form of A-Diakoptics was developed in 2018 as the result of EPRI’s participation in the SunShot National Laboratory Multiyear Partnership Program (SuNLaMP). SuNLaMP is a government initiative that brought together national laboratories (such as Sandia National Laboratories and National Renewable Energy Laboratory), the private research sector (EPRI), academia (Georgia Institute of Technology), and private industry (CYME International T&D) to investigate novel technologies for accelerating QSTS simulations in power systems using modern computing architectures [13]. From there, the general expression for describing the interactions between the actors (subcircuits) and the coordinator is as follows:

$$E_T=Z_{TT}{I}_{0\left(n-1\right)}-Z_{TC}{Z}_{CC}^{-1}{Z}_{CT}{I}_{0\left(n\right)}$$

(4)

where E_T is the total solution of the system (the voltages in all the nodes of the system). I₀ is the vector containing the currents injected by the PC elements, and the time instants n and n + 1 are discrete time instants to describe the different times in which the vector of currents is calculated [6]. Z_TT is the trees matrix and contains the admittance matrices for all the subcircuits contained in the interconnected system after the partitioning. The form of Z_TT is as follows:

$$Z_{TT}=\left[\begin{array}{cccc}{\left[Y_1\right]}^{-1}& 0& \dots & 0\\ 0& \ddots & & ⋮\\ ⋮& & \ddots & ⋮\\ 0& \dots & \dots & {\left[Y_n\right]}^{-1}\end{array}\right]$$

(5)

$$n=0,1,2\dots \# of sub-circuits$$

The subcircuits contained in Z_TT are not interconnected in any way. The connections are through external interfaces that define the relationship between them. Z_TC, Z_CT, and Z_CC are interfacing matrices for interfacing the separate subsystems using a graph defined by the contours matrix (C), as described in [6]. Z_TTI_0(n−1) corresponds to the solutions delivered when solving the subsystems, and Z_TCZ_CC⁻¹Z_CTI_0(n) is the interconnection matrix to find the total solution to the power flow problem.

The form of Z_TT proposes that multiple OpenDSS solvers can find independent partial solutions that when complemented with another set of matrices, can calculate the voltages across the interconnected power system. However, in this approach, the interconnection matrix is a dense matrix that operates on the injection currents vector provided by the latest solution of the subsystems.

Nowadays, sparse matrix solvers are very efficient, and therefore, operating with a dense matrix is not desirable. The aim of this part of the project is to simplify the expression that defines the interconnection matrix to find a sparse equivalent that can be solved using the KLUSolve [14] module already employed in OpenDSS.

Assuming that the times in which the subsystems are solved, and the interconnection matrix is operated, are the same and fit into the same time window (ideally), A-Diakoptics can be reformulated as

$$E_T=Z_{TT}{I}_{0\left(n\right)}-Z_{TC}{Z}_{CC}^{-1}{Z}_{CT}{I}_{0\left(n\right)}$$

(6)

Additionally, in OpenDSS, the interconnected feeder is solved using

$$E_T={\left[Y_{II}\right]}^{-1}{I}_0$$

(7)

where Y_II is the Y_Bus matrix that describes the interconnected feeder. Equating Equations (6) and (7), the expression is transformed into

$${\left[Y_{II}\right]}^{-1}{I}_{0\left(n\right)}=Z_{TT}{I}_{0\left(n\right)}-Z_{TC}{Z}_{CC}^{-1}{Z}_{CT}{I}_{0\left(n\right)}$$

(8)

This new expression can be taken to the admittance domain since Z_TT is built using the YBus matrices that describe each one of the subsystems created after tearing the interconnected feeder. The new equation is reformulated as

$${\left[Y_{II}\right]}^{-1}{I}_{0\left(n\right)}={\left[Y_{TT}\right]}^{-1}{I}_{0\left(n\right)}-Z_{TC}{Z}_{CC}^{-1}{Z}_{CT}{I}_{0\left(n\right)}$$

(9)

$${\left[Y_{II}\right]}^{-1}={\left[Y_{TT}\right]}^{-1}-Z_{TC}{Z}_{CC}^{-1}{Z}_{CT}$$

(10)

Then, the interconnection matrix can be reformulated as

$$Z_{TC}{Z}_{CC}^{-1}{Z}_{CT}={\left[Y_{XX}\right]}^{-1}={\left[Y_{TT}\right]}^{-1}-{\left[Y_{II}\right]}^{-1}$$

(11)

This new formulation proposes that the interconnection matrix is equal to an augmented representation of the link branches between the subsystems. To support this conclusion, take the simplification proposed by Happ in [15], where:

$$e_{c^{\prime }}=-Z_{CT}{I}_{0\left(n\right)}$$

(12)

$$Z_{CT}=C^T{Z}_{TT}\to e_{c^{\prime }}=-C^T{Z}_{TT}{I}_{0\left(n\right)}$$

(13)

This expression is similar to the partial solution formulation of the Diakoptics equation but includes the contours matrix; then, if the partial solution is called E_T(0), it is possible to define

$$e_{c^{\prime }}=-C^T{E}_{T\left(0\right)}$$

(14)

On the other hand, Z_TC, which is the non-conjugate transposed of Z_CT, is calculated as follows:

$$Z_{TC}=Z_{TT}C$$

(15)

Replacing the new equivalences for Z_CT and Z_TC the equation proposed in (4) is reformulated as follows:

$$E_T=Z_{TT}{I}_{0\left(n\right)}+Z_{TT}{I}_c$$

(16)

Or, in terms of sparse matrices,

$$E_T={\left[Y_{TT}\right]}^{-1}{I}_{0\left(n\right)}+{\left[Y_{TT}\right]}^{-1}{I}_c={\left[Y_{TT}\right]}^{-1}\left(I_{0\left(n\right)}+I_c\right)$$

(17)

where:

$$I_c=-C{Z}_{CC}^{-1}{C}^T{E}_{T\left(0\right)}$$

(18)

In this equation, its partial and complementary solutions use information calculated from the interconnection matrix Z_CC. Z_CC is a small matrix that contains information about the link branches, and the calculation of I_C does not represent a significant computational burden in medium- and large-scale circuits. With this approach, all the matrix calculations are made using sparse matrix solvers, reducing the computational burden.

The architecture for implementing this simplified approach considers the optimization of computational resources for reducing the burden of the operations around the algorithm. This optimization comprehends memory handling, actor coordination, and message structure between actors, and it is described in detail in [16]. Figure 1 describes the implementation of the simplified A-Diakoptics method in OpenDSS 9.3 and later.

Figure 1 presents how the main loops of the OpenDSS algorithm are decomposed into the set of distributed solvers implemented as actors. Each actor (including the coordinator) is executed in its own CPU and counts with its own memory and hardware. All of the actors are executed in parallel and concurrently in an asynchronous parallel computing environment.

Actors remain on standby when not executing any job. A message sent to an actor triggers an event that sets the busy flag for the actor in the global context and leads the actor to perform the job given in the message. Once the actor has finished the job, it updates its state by marking a ready flag in the global context, making it visible to all of the other actors. When the Y_BUS matrices within actors >2 are smaller than the interconnected system, it is expected that the solution time required for solving the subnetworks is substantially lower than that required for solving the interconnected model, as discussed in [5,13]. The calculation of I_C is very fast, given the size of the mathematical problem proposed; at the same time, it solves and avoids moving values between actors by using pointers, adding a negligible overhead to the solution process in parallel.

With the architecture proposed in this implementation of A-Diakoptics, the algorithm becomes available for all QSTS simulation modes. The previous version of the algorithm was only available for yearly simulation. The new implementation is available for the snap, direct, yearly, daily, time, and duty cycle simulation modes in OpenDSS [17].

Centralizing the results on a single memory space (Actor 1) makes it possible to insert monitors and energy meters and request simulation information directly from the model using a single actor, a feature that was not available in the previous version of the A-Diakoptics suite.

The previous considerations are true only if the matrix Z_CC is a fraction of the size of the largest subnetwork in the system, which is the one representing the largest computational burden in the parallel solver, as explained in [16].

3. Results

This section provides a case publicly available at https://sourceforge.net/p/electricdss/code/HEAD/tree/trunk/Version8/Distrib/Examples/ADiakoptics/EPRI_Ckt5-G/ (accessed on 1 January 2022) since October 2021, for illustrating the A-Diakoptics features implemented in OpenDSS v9.2 and later. The features to discuss in this section are as follows:

• Simulation fidelity;
• Simulation performance.

3.1. The Test Case: EPRI Circuit 5

Circuit 5 and the other two circuits (7 and 24) are models EPRI provided to the IEEE test feeder database at the IEEE PES Distribution System Analysis Subcommittee [18]. These feeders are real feeder models from real utilities that have been modified for public use and represent medium- and large-scale power systems. Figure 2 shows Circuit 5 and

Table 1. Circuit 5 technical features.

Feature Name	Value
Number of buses	2998
Number of nodes	3437
Total active power	7.281 MW
Total reactive power	3.584 Mvar

presents the technical features of this model.

Figure 2 indicates the link branches for the two-zone and four-zone partitioning with blue and orange arrows, respectively.

Table 2. Partition statistics for Circuit 5 feeder, two zones.

Name	Value
Circuit reduction (%)	32.24
Maximum imbalance (%)	52.43
Average imbalance (%)	26.21

and

Table 3. Partition statistics for Circuit 5 feeder, four zones.

Name	Value
Circuit reduction (%)	45.94
Maximum imbalance (%)	88.11
Average imbalance (%)	53.75

present the partition statistics for two- and four-zone partitioning, respectively. The partition method for obtaining these statistics was the automated option through MeTIS [19]. For more extensive partitioning (five, six, and seven zones), the manual procedure specifying the link branches was used [17].

The partition statistics reveal a moderate imbalance in the case of two zones that is more pronounced in the case of four zones. The imbalances are expected given the representation of single-phase low-voltage networks within the model, making it more difficult for the tearing algorithm to balance the node distribution between zones. The calculation methods for the partition statistics and their interpretation details are presented in [17].

3.2. The Simulation Fidelity and Performance

Given the number of nodes present in this model, to study the simulation, fidelity monitors were located close to the feeder’s head and at the feeder’s end to check the voltage levels at those two points. The voltages were measured in all phases, and the error was estimated per phase in a yearly simulation (8760 h). Figure 3 and Figure 4 present the results obtained for two and four zones, respectively.

As expected, the simulation fidelity comparison reports that the results obtained in the new implementation differ from the original by a negligible margin. In this simulation, each load has a load profile assigned and, based on the previous experiences with controllers modifying the local Y_BUS matrices at the distributed solvers, the capacitor bank controllers were disabled for the simulation.

The good fidelity when comparing the interconnected model and its partitioned counterparts is a good indicator for this implementation of A-Diakoptics, despite the issues with some controller types. Additionally, it confirms the method’s scalability for allowing a parallel simulation using multicore computing architectures.

The next step was to evaluate the simulation performance. For evaluating the simulation performance, the controllers based on Y_BUS matrix modification were disabled based on the experiences reported in [16].

The yearly simulations for Circuit 5 were conducted on a laptop with an Intel^®® Core™ i7-10750H CPU @ 2.60 GHz, 2592 MHz, 6 Core(s), 12 Logical Processor(s), and 16 GB RAM, which is a standard computer operating under a standard Windows 10 operating system. During the simulation, programs such as antivirus, firewalls, and other business/system-related applications were under execution.

Figure 5 shows the solution times in microseconds for a yearly simulation with the interconnected model and its partitions, using two, four, five, six, and seven zones. The time was measured using the internal process timer included with OpenDSS.

This timer registers the time required for each operation within OpenDSS in real time. Figure 6 displays the simulation time reduction as a percentage of the simulation time for the reference model. Simulation computing time reduces significantly when the number of zones increases; however, the time reduction is the same for the six- and seven-zone partitions, suggesting that the maximum simulation time reduction is reached at that point or that better tearing points need to be located across the feeder to achieve a better reduction.

Figure 7 presents the partition statistics for each case and indicates that the better the circuit reduction, the better the simulation time reduction. Nevertheless, when the circuit reduction is contrasted with the maximum and average imbalances, the parallel solvers are not balanced as expected.

The imbalance level spikes when tearing the model into four zones and goes down slightly when using a larger number of zones. This experimentation signals that the parallel solvers do not necessarily have to be balanced to achieve a considerable simulation time reduction, even though a balanced state is a desirable condition.

Instead, this experiment reveals that simplifying the circuit model complexity requires understanding where the complexity resides in the model, and what challenges exist when locating the tearing points. This criterion can be constructed based on the features of areas that require more exhaustive computer power to be solved given the load models, generators, and topology hosted in the area. This topic is deeply discussed in [16].

4. Conclusions

Spatial parallelization has proven to be a technique for accelerating power system simulations by reducing the circuit model’s complexity. The circuit model’s complexity is reduced by tearing the model into several subnetworks that can be solved faster than the original interconnected model. The partial solutions are put together through a tensor-based interconnection approximation, which delivers an accurate solution that matches the original with a low error margin.

One method for performing spatial parallelization is Diakoptics, which, in combination with a framework for providing inconsistency robustness within a parallel computing architecture, is called A-Diakoptics. The previous implementation of this technique in EPRI’s open-source distribution simulator software OpenDSS provided limited access to this technique, given that the implementation aimed to cover the needs of a particular project.

Through the project presented in this document, the previous implementation of A-Diakoptics in OpenDSS was reformulated using a simplified approach, optimizing critical areas of the algorithm and the mathematical formulation. Many other improvements at the implementation level were performed in this new version of A-Diakoptics to facilitate access to electrical variables through monitors and meters.

This implementation of A-Diakoptics was validated using yearly QSTS simulations in terms of fidelity and performance using as reference the interconnected model. The simulation results revealed the benefits of modernizing sequential power simulation tools into parallel processing, which is the standard for computing architectures nowadays. This paper introduced A-Diakoptics as a method for achieving such modernization and OpenDSS, given that it is distributed as an open-source project using the internet, presents a practical example of how this modernization can be implemented.

A simulation test case was used for illustrating the computational time gains, comparing the sequential and parallel performance on a yearly QSTS. The computing time gains were evaluated for an incremental number of partitions, displaying and discussing the benefits and limits of the algorithm. This test case is publicly available on the internet for the reader to evaluate if desired.

The results presented confirm the performance improvement in QSTS simulations for medium- and large-scale power system models. They also provide guidance on the objectives of circuit tearing for reducing the circuit model complexity to accelerate QSTS simulations.

The test cases presented in this document and the A-Diakoptics capabilities are available in OpenDSS version 9.3, which can be downloaded from the internet. In addition to the test case presented in this paper, five more test cases offering different levels of complexity are available at OpenDSS/Code/[r3354]/trunk/Version8/Distrib/Examples/A-Diakoptics (sourceforge.net) (accessed on 1 January 2022). These cases and their purpose are discussed in [16] and serve as sources of reference for applying A-Diakoptics with OpenDSS in other applications.

Additional to the methods presented in this paper, others using different parallelization techniques and approaches have also been investigated as the result of the SunShot National Laboratory Multiyear Partnership Program (SuNLaMP). The findings, conclusions, and a more extended explanation of these methods can be found in [13]. For future implementations and demonstrations on methods for accelerating simulations, we expect to continue to use OpenDSS as the open-source project for serving as a proof of concept for this type of research, demonstrating viable techniques and methods for modernizing power simulation tools.