Exploring Stability and Accuracy Limits of Distributed Real-Time Power System Simulations via System-of-Systems Cosimulation

Electromagnetic transients (EMT) is the most accurate, but computationally expensive method of analyzing power system phenomena. Thereby, interconnecting several real-time simulators can unlock scalability and system coverage, but leads to a number of new challenges, mainly in time synchronization, numerical stability, and accuracy quantification. This study presents such a cosimulation, based on digital real-time simulators (DRTS), connected via Aurora 8B/10B protocol. Such a setup allows to analyze complex and hybrid system-of-systems whose resulting numerical phenomena and artifacts have been poorly investigated and understood so far. We experimentally investigate the impact of IEEE 1588 precision time protocol synchronization assessing both time and frequency domains. The analysis of the experimental results is encouraging and show that numerical stability can be maintained even with complex system setups. Growing shares of inverter-based renewable power generation require larger and interconnected EMT system studies. This work helps to understand the phenomena connected to such DRTS advanced cosimulation setups.


Exploring Stability and Accuracy Limits of Distributed Real-Time Power System Simulations via
System-of-Systems Cosimulation Luca Barbierato , Enrico Pons , Ettore Francesco Bompard , Vetrivel S. Rajkumar , Peter Palensky, Lorenzo Bottaccioli , and Edoardo Patti Abstract-Electromagnetic transients (EMT) is the most accurate, but computationally expensive method of analyzing power system phenomena.Thereby, interconnecting several real-time simulators can unlock scalability and system coverage, but leads to a number of new challenges, mainly in time synchronization, numerical stability, and accuracy quantification.This study presents such a cosimulation, based on digital real-time simulators (DRTS), connected via Aurora 8B/10B protocol.Such a setup allows to analyze complex and hybrid system-of-systems whose resulting numerical phenomena and artifacts have been poorly investigated and understood so far.We experimentally investigate the impact of IEEE 1588 precision time protocol synchronization assessing both time and frequency domains.The analysis of the experimental results is encouraging and show that numerical stability can be maintained even with complex system setups.Growing shares of inverter-based renewable power generation require larger and interconnected EMT system studies.This work helps to understand the phenomena connected to such DRTS advanced cosimulation setups.

I. INTRODUCTION
S IGNIFICANT scientific effort has recently been conducted on computer-aided power system analysis for the design, development, and testing of future power system designs.This has resulted in multiple domain-specific simulation tools that can capture functionalities and behavior of different power system components with a high-precision accuracy [1].In particular, numerical real-time simulation has arisen as a critical tool for modern-day power system planning, design, and operation [2].Real-time simulation refers to a digital twin of a real-world power system, which is simulated in wall clock time, to mimic the behavior of its real-world physical counterpart.These real-time simulations are conducted using small, discrete, and constant time steps, usually in the order of a few microseconds.The real-time constraints are susceptible to the targeted analysis.For instance, transient stability analysis is carried out with phasor simulations of around 10 ms of time step duration.EMT analysis instead requires a smaller time step duration (i.e., tens of microseconds) to detail the dynamics of large power systems [3].In particular for EMT, innovative multiprocessor architecture (e.g., IBM Power8) and FPGA have been used to deploy a DRTS.DRTSs are solution to hardware-accelerated EMT analysis [4] respecting real-time restrictions.These simulators allow fast analog and digital I/O to connect in a closed-loop real-world equipment, allowing PHIL to test its functionalities in a safe experimental testbed.The major limitation of the commercial DRTS is the significant amount of computational power required to run detailed EMT power system analysis, generating limitations on the dimension of simulated models [4], [5].
This limitation has prompted power systems and ICT experts to collaborate with each other to connect two or more DRTS by leveraging on novel cosimulation techniques [6], communication protocols, and standards [7].These techniques permits This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/breaking up a larger PSUT into smaller submodels, with each being executed on a dedicated DRTS.Such an approach involves the usage of IA, such as the ITM [8], [9].This IA is inspired by PHIL applications and involves representing the DUT as a current source, with the network simulation as a voltage source [10], [11].Thereby, a larger system is decomposed into a set of smaller subsystems, thereby, resulting in a SoS approach.This approach is demonstrated in recent work such as [12], wherein two university labs representing local energy communities have been interconnected.Similarly, in [13] and [14], the ITM is used to interface geographically separated DRTS and PHIL equipment.Furthermore, such setups have also been extensively used for cybersecurity studies and resilience analysis, as noted in [15] and [16].
Typically, such distributed real-time simulations are achieved through high-bandwidth communication protocols (e.g., IEEE802.3) to exchange interface signals (i.e., voltages and currents) between the submodels/subsystems.For instance, multiple joint research experiments, spanning continents, have been carried out to overcome the limitations of individual capabilities of real-time laboratories, such as the infrastructures proposed in [17] and [18].The objective of these experiments is to virtually interconnect DRTSs and PHIL setups in geographically distributed laboratories via the Internet, allowing simulation and testing of future power systems.Thus, a virtual research infrastructure can be formed.
In [19], the authors present a detailed overview and analysis of case studies related to cosimulation application, called Real-Time Coupling of Geographically Distributed Research Infrastructures, subsequently leading to the paradigm of geographically distributed simulations.However, these interconnections provoke numerical instabilities and accuracy issues due to inherent communication latencies between interconnected DRTSs, as noted in [20] and [21].The proposed solutions have been applying time-delay compensation methods, such as the one in [22] and [23].However, these solutions severely limited the scope of EMT studies.In fact, such compensation methods are still geared toward power PHIL setups but with significantly slower dynamics than EMT simulations.Another associated drawback of distributed real-time cosimulations is time synchronization between DRTSs.This leads to misalignment and inaccurate logging of results due to the probabilistic nature of the communication medium [24], [25].Therefore, there is a need for a scalable and accurate solution to interconnect DRTSs for distributed real-time EMT simulations while adhering to real-time constraints.
This article presents a novel application of a locally distributed hybrid digital real-time cosimulation infrastructure that allows a point-to-point connection of two DRTSs via the high-bandwidth communication protocol Aurora 8B/10B.Aurora reduces the consequence of uncertainty brought about by communication latency by adhering to real-time constraints in a cosimulated environment.The proposed solution synchronizes DRTSs' interconnections through the IEEE 1588 PTP [26] to adjust cosimulation execution and logging of synchronized simulation results.Through this setup, we seek to explore the accuracy and stability limits of a distributed real-time cosimulation.
With respect to our previous works [27], [28], the key contributions and novelties proposed in this manuscript are summarized in the following.
1) This solution presents a more general time synchronization strategy that allows to concurrently implement different IEEE 1588 PTP stack configurations following the same GPS clock reference, required because different DRTS could implement different PTP profiles and delay mechanisms.2) Different DRTS in our infrastructure can be interconnected for cosimulating a power system scenario, not only RTDS.In our laboratory setup, we propose RTDS and OPAL-RT, which are the main DRTS commercial solutions for power system analysis.Moreover, the infrastructure is flexible to integrate other DRTS solutions (e.g., Speed Goat and NI PXI) that implement the Aurora and IEEE 1588 PTP. 3) A comprehensive description of the intratime step operations of RTDS and OPAL-RT is given, explaining how the different DRTS solutions manage the model calculation, the CPU, the CPUs/FPGA communications, and the management of the Aurora link to better understand where the communication latency of the cosimulation impact is generated.4) The laboratory setup has been tested on four different configurations involving either RTDS, OPAL-RT, or both, where the two power subsystems have been run separately, namely RTDS-RTDS, OPAL-OPAL, RTDS-OPAL, and OPAL-RTDS.5) The sequence operations of each of the four configurations are described to explain the communication latency generated by the application of our cosimulation infrastructure.This allows a precise understanding of the contribution of the experienced latencies, instead of measuring the latency by running a simulation as in our previous works.6) Time-domain and frequency-domain results present different latencies effects between the different infrastructure configurations.The structure of this article is described as follows.The computing capabilities and system size simulation limits of the DRTS are discussed in Section II.The locally distributed hybrid digital real-time cosimulation infrastructure is described in Section III, along with the laboratory setup implemented for the experiments.The experimental results are shown in Section IV that includes the analytical description of the communication latency in the four proposed configurations, the time-domain accuracy for steady-state and transient scenarios, and the frequency-domain accuracy for the steady-state scenario.Finally, Section V concludes this article.

II. SYSTEM SCALABILITY AND MOTIVATION
As briefly introduced in Section I, the main limitation of commercial DRTS (e.g., OPAL-RT and RTDS Technologies) is the significant amount of computational power required to run detailed EMT models, thereby limiting the overall size of power systems that can be accurately simulated.The maximum size of the power system that can be simulated in one real-time simulator cannot be easily and precisely assessed, as it depends mainly on the following: 1) the size of the time step used in the simulation; 2) the model of the real-time simulator, software platform, and solution algorithms; 3) the complexity of the control systems in the model; 4) the number of I/O channels needed.A very rough figure is that on a single RTDS NovaCor chassis with one activated core, the simulated system can have a maximum of 90 single-phase buses.With more activated cores, the maximum network size per chassis is 600 single-phase buses.The actual size then depends on the complexity of the models used for the power system and control system components.A very rough figure for OPAL-RT is that with the eMEGASIM software, it is possible to simulate networks with a maximum of 60-70 single-phase nodes per activated core with a time step of 50 μs.
By coupling several DRTS in a cosimulation infrastructure, it is possible to sum up the computational power of the different simulators, and therefore, simulate larger networks.Unfortunately, it is not possible to compare the simulation results for large networks in a monolithic simulation with the cosimulation, as the monolithic system would not run on a single DRTS due to a lack of computational resources.A detailed discussion on the scalability of such setups is beyond the scope of this work.

III. METHODOLOGY
The proposed locally distributed digital real-time cosimulation infrastructure is a hybrid architecture that makes use of Aurora 8B/10B, a serial fiber optic link protocol, and IEEE1588 PTP, which is commonly used to synchronize internal reference clocks throughout a computer network, to enable the serial interconnection of various commercial DRTS (e.g., OPAL-RT and RTDS technologies).The architectural description of the proposed infrastructure is presented in Fig. 1 and consists of three main layers: the GPS synchronization layer, the DRTS, and the power system cosimulation layer.The rest of this section presents every single layer and its main components and the laboratory setups that are used to test the infrastructure functionalities.

A. GPS Synchronization Layer
This layer guarantees the complete infrastructure time synchronization by exploiting a GPS clock.The GPS is a global navigation satellite system that allows worldwide time synchronization by broadcasting a time reference clock synchronized to the GPS atomic clock technology, in addition to providing geolocation.The GPS clock is a hardware component that receives the GPS signal and synchronizes its internal reference clock to the atomic time received by the GPS satellites without needing a nearby atomic clock, ensuring synchronization among devices all over the world.Once synchronized with the GPS, this equipment offers different local synchronization protocols, such as IRIG-B, 1PPS, and the IEEE 1588 PTP.Our infrastructure uses the IEEE1588 PTP to synchronize both internal reference clocks of the interconnected DRTS.By exploiting the PTP stack functionalities, the GPS clock sets its operation in master mode to control other slaves by broadcasting PTP packets on a LAN.The end user can choose among a set of standard PTP profiles (e.g., default profile) that serve different final uses.Moreover, the PTP stack offers two different delay mechanisms, namely, the E2E and the P2P, that provides PTP functionalities to the two homonymous types of transparent clocks.The E2E mechanism allows the calculation of the overall latency among the master and slave nodes ensuring the fastest synchronization with low precision.Moreover, this delay mechanism ensures the most interoperable version of the PTP standard, reducing the limitations on hardware that can be used in the PTP network.The P2P delay mechanism instead allows a fine-grained and precise latency calculation of each hop of the network path between the master and the slave node.However, it reduces the interoperability of different hardware since they must implement the overall PTP standard protocol stack.
Furthermore, PTP offers the following two different protocols to allow communication among master and slave nodes: 1) Layer 2 (i.e., IEEE802.3 or Ethernet) that ensures a fast synchronization in the same network segment; 2) Layer 3 (i.e., IP) that allows the time synchronization of wider networks.Depending on the LAN configuration of the laboratory, the end user could choose among the two layers.Apart from this configuration, the GPS clock is interconnected to the LAN using two RJ45 Ethernet interfaces, one serving the E2E, and the other the P2P delay mechanism.

B. DRTS Layer
This is the central layer of the proposed architecture.It enables the connection of two DRTS by exploiting a physical point-to-point bidirectional serial communication.This layer takes advantage of Aurora 8B/10B, a point-to-point full-duplex multimode optical fiber serial link with a high bandwidth of 2 to 5 Gb/s.As a matter of fact, Aurora is the fastest on-board protocol available on commercial DRTS, guaranteeing the lowest latency.To enable Aurora protocol, each DRTS must occupy one of the available Aurora SFP ports for the data exchange.
Each DRTS is also equipped with an IEEE1588 PTP synchronization board since the time synchronization among the DRTSs involved in the proposed infrastructure is ensured by the IEEE1588 PTP.PTP boards are set as slave-only PTP nodes.Moreover, they are interconnected with the GPS synchronization layer with an RJ45 Ethernet cable to receive PTP packets coming from the PTP master node (i.e., GPS clock) and align the DRTS internal reference clock to the reference master clock.Finally, each interconnected DRTS implements its own logic to fulfil the real-time simulation of the assigned power subsystem by including simulation blocks that enables the Aurora link and the PTP synchronization in the compiled models.

C. Power System Cosimulation Layer
The application of the logical SoS split of a PSUT is implemented in this layer via the ITM IA [27].ITM is commonly adopted as interface to interconnect a physical DUT to a simulated ROS in PHIL applications.It is the simplest choice to set up a PHIL system and it has been chosen to simply and efficiently split a PSUT into a source and load power subsystems modeled for each DRTS, where each resulting subsystem runs.The ITM IA illustrated in Fig. 2(b) makes use of a controlled voltage source in the load circuit of the power subsystem to replicate the voltage v A measured in the source circuit (i.e., v A ) and a controlled current source in the source circuit of the power subsystem to replicate the current i B measured in the load circuit (i.e., i B ).
Additionally, it applies a latency that is correlated to the delay experienced by the voltage and current value interchanged between DRTS racks 1 and 2 to affect the load power subsystem from the source power subsystem (i.e., T D 1 ) and vice versa (i.e., T D 2 ).The open-loop transfer function of the ITM IA described in the following equation is obtained by exploiting the ITM equivalent block diagram [28]: where Z A is the source impedance; Z B is the load impedance; T D 1 is the latency affecting the voltage signal exchanged between the source and the load circuits; and T D 2 is the latency affecting the current signal exchanged between the load and the source circuits.This equation describes the frequency stability of the ITM IA solution state space.By applying the Nyquist criterion, the stability of the PSUT is ensured when the ratio Z A /Z B is lower than 1.Furthermore, a large latency T D 1 + T D 2 could provoke nonlinearities (i.e., phase shifts) impacting the accuracy of the overall system in time and frequency domains, also in a stable region ensured by the Nyquist criteria.

D. Laboratory Setups
The proposed laboratory setups reduce significantly the latency between two DRTS chosen among the available commercial solutions for electric grid analysis, i.e., OPAL-RT OP5700 (OPAL) and RTDS Technologies NovaCor (RTDS).Following the architectural description in Fig. 1, the GPS synchronization layer is managed by the Meinberg microSync HR102HQ GPS clock with an IEEE 1588-2008 v2 Default Layer 2 profile.This GPS clock offers four configurable Ethernet interfaces.Since RTDS and OPAL synchronization boards implement two different TC, port 2 has been configured to manage E2E TC and port 3 instead to manage P2P TC.The GPS clock is the master of the IEEE1588 PTP network, and the two DRTS synchronization boards are slaves to ensure a proper set up of the synchronization with the reference master clock.
RTDS racks must include the GTSYNC card to deal with the IEEE1588 PTP network.GTSYNC is a peripheral board interconnected with the core rack that enables different synchronization techniques (e.g., IEEE 1588 PTP, 1PPS, and IRIG-B).Moreover, the GTSYNC only accept P2P TC configuration of the IEEE1588 PTP stack.So, both RTDS racks' GTSYNC cards are connected with a RJ45 Ethernet cable to port 3 of the Meinberg GPS clock exploiting an Ethernet switch.During the design of the RSCAD draft, the GTSYNC block is a prerequisite to enable the synchronization of the GTSYNC card.OPAL instead provides the Oregano card that communicates with the core rack through the PCIe bus.Like the GTSYNC, the Oregano card allows the same set of synchronization protocols.However, the TC configuration of the IEEE1588 PTP stack must be E2E.So, both OPAL racks' Oregano cards are connected with an RJ45 Ethernet cable to port 2 of the GPS clock using the Oregano Syn1588 Ethernet switch to reduce latency of the PTP E2E packets.The time regulation schema is insignificant because both DRTSs evolve their simulations with a distinct time-regulating schema that maintains the correct event ordering while adhering to the appropriate real-time restriction.Because the proposed circuit for testing the connections is a basic source-load circuit, the source power circuit will cause the load power circuit simulation to begin.Because the IEEE 1588 PTP stack maintains the correct time synchronization schema, results from DRTS racks are actually aligned using the internal clock time as a reference.
Since each commercial DRTS solution implements its own numerical solver and implementation of the Aurora 8B/10B, different configurations of the DRTS layer are proposed in Fig. 3 to take into account each possible DRTS combination.Fig. 3(a) presents two RTDS NovaCor that have been interconnected with an optical fiber cable of 25 m by leveraging on Aurora protocol.Fig. 3(b) instead present two OPAL racks interconnected by an identical optical fiber cable of 25 m by leveraging on Aurora protocol.Finally, an RTDS rack and an OPAL rack have been coupled together as depicted in Fig. 3(c).This interconnection has been analyzed in both directions, considering rack 1 as the generation source and rack 2 as the load of the proposed power system cosimulation scenario.
Fig. 4 presents the sequence operation diagrams of a time step for both RTDS and OPAL are presented to highlight the management of the Aurora read (RD) and write (WR) operations.In Fig. 4(a), RTDS executes after each time step starting time an instantaneous WR operation followed by an RD.At the end of the RD, received variables are exchanged with the control signal core that runs its operations in parallel with the power system component cores.At the end of the time step, cntrol signal core and power system component cores exchange both the control signal variable and the network signals representing Aurora is enabled in RTDS chassis through RSCAD adding the Aurora block.It permits to select the specific SFP port, the assigned processor, the computation priority, the frame structure (i.e., exchanged control signals with a maximum width of 128), and the sequence number blocking property.In OPAL systems instead, the Aurora protocol is enabled through RT-LAB by selecting the proper SFP communication block that defines the FPGA DataIn and DataOut port numbers and enables a variable width of the exchanged signals with a maximum of 255 doubles exchanged.The Aurora link configuration (e.g., SFP transceiver port) instead is defined in the FPGA bitstream, which is uploaded during the scenario loading operations.
Finally, the power system cosimulation layers implements the simple ITM circuit presented in Fig. 2(b) by splitting the monolithic circuit [see Fig. 2(a)] into a source circuit A and a load circuit B. This simple scenario permits to precisely calculate the latency among the different configurations of the DRTS layer, and the accuracies in time and frequency domains of the cosimulated results by comparing them with the standalone monolithic results.

IV. EXPERIMENTAL RESULTS
This section presents the experimental results of the simple electric circuit applied to the proposed digital real-time cosimulation infrastructure.The circuit in Fig. 2(b) has 1, respectively.To return feedback to the source part A, a metering point retrieves i B current to then send it back through the Aurora link.This value is received by source part A and imposed through the controlled current source i B .The ITM IA circuit has been tested for all the possible combinations, namely, RTDS-RTDS, OPAL-OPAL, RTDS-OPAL, and OPAL-RTDS (see Fig. 3).All the configurations exploit a 25-m long standard full-duplex multimode optical fiber cable as physical media to interconnect both racks.
In the following subsections, we discuss the experimental results on the following.
1) The ITM IA latency calculation that describes the time step contributions to the overall latency experimented by the ITM IA circuit for each of the digital real-time simulation layer configurations.2) The ITM IA time-domain accuracy analysis that compares the standalone voltages signals with the cosimulated ITM IA case in the stable region (Z B = 500 Ω) and near the instability region (Z B = 50.5Ω) for each configuration of the infrastructure.
3) The ITM IA frequency-domain accuracy analysis that, similarly to the time-domain case, compares the standalone frequential contents with the cosimulated solution ones for both regions and all infrastructure configurations.For each section, the time step duration T Sim has been changed from 20 to 500 μs to highlight possible dependencies from the time step duration that could influence the overall latency.

A. ITM IA Latency Calculation
For each infrastructure configuration, the overall observed latency the ITM IA circuit solution is described using its sequence operation diagrams in Fig. 5.More in-depth, these sequence diagrams describe the contribution of each single time-step, highlighting internal operations executed by each DRTS and the data exchange among the interconnected racks.Depending on the infrastructure configuration, these operations may or may not contribute to the overall latency represented by T D 1 + T D 2 in (1).The generated phase shifts impact both the time-and frequency-domain accuracy of the ITM IA circuit solution.
1) RTDS-RTDS: By exploiting RSCAD software libraries, Aurora blocks are set with the sequence number blocking configuration activated on port 24 for both RTDS racks.The overall calculated latency by the ITM IA application is 5T Sim for all T Sim values.The contributions are depicted in the sequence operation diagram in Fig. 5(a).From t = 0 to t = T Sim , the power system component core (PSCC) of RTDS rack 1 calculates v A and passes its value to the control system core (CSC).From t = T Sim to t = 2T Sim , rack 1 retrieves v A value from the CSC and sends it through the Aurora link by executing a WR operation.In the same time interval, RTDS rack 2 receives v A value from the Aurora link with an RD operation and gives it to the CSC that is in charge to send its value to the PSCC at the end of the time interval.From t = 2T Sim to t = 4T Sim , rack 2 calculates the current i B by applying v A voltage to the controlled voltage source v A , and finally, moves i B from the PSCC to the CSC.These operations take 2T Sim in total.From t = 4T Sim to t = 5T Sim , rack 2 retrieves the i B value from the CSC and executes an Aurora WR operation.In the same time interval, rack 1 receives its value by executing an Aurora RD operation and passes its value to the CSC.To conclude, the effect of i B on v A calculation requires another T Sim to take effect on source part A of the circuit in Fig. 2(b).The ITM IA latency is calculated from the first v A calculation in rack 1 at t = T Sim to the effect of i B to v A calculation at t = 6T Sim , resulting in a total of 5T Sim .
2) OPAL-OPAL: By exploiting RT-LAB software libraries, the SFP block is integrated into both models of the ITM circuit.SFP blocks are set on port SFP00 for both OPAL racks.The calculated latency is 2T Sim for all T Sim values and its contributions are depicted in Fig. 5(b).From t = 0 to t = T Sim , the CPU of the OPAL rack 1 calculates v A and provides its value to the FPGA that executes the Aurora WR operation.In the same time interval, the FPGA of OPAL rack 2 receives v A value by executing the Aurora RD operation.Right at the beginning of the time interval between t = T Sim and t = 2T Sim , rack 2 moves v A value from the FPGA to the CPU.Subsequently, the CPU of rack 2 executes the calculation to retrieve the current i B , moves its value to the FPGA, and sends it through the Aurora link with an Aurora WR operation.To conclude, in this time interval, rack 1 receives i B value and stores it in the FPGA by executing an Aurora RD operation.In the last time intervals from t = 2T Sim to t = 3T Sim , i B moves from the FPGA to the CPU to allow the calculation of the effect on v A on the source circuit.The v A value is updated at the end of this time interval as a result of the circuit numerical solution.So, the ITM IA latency is calculated from the first v A calculation in rack 1 at t = T Sim to the effect of i B to v A calculation at t = 3T Sim , resulting in a latency of 2T Sim .
3) RTDS-OPAL: The source circuit of Fig. 2(b) in the RSCAD draft and load counterpart in the RT-LAB model of the two previous cases are used in this first hybrid configuration.The physical interconnection of the optical fiber link has been changed from port 23 of the RTDS rack 1 to port SFP01 of the OPAL rack 2. The overall latency is 4T Sim for all T Sim values.The contributions are presented in Fig. 5(c).From t = 0 to t = 2T Sim , RTDS rack 1 executes the same operations described in Section IV-A1, presenting as a result v A at t = T Sim .OPAL rack 2 receives v A value by executing an Aurora RD operation at the end of the time interval from t = T Sim to t = 2T Sim .At the very beginning of the time interval from t = 2T Sim to t = 3T Sim , OPAL rack 2 moves v A value from the FPGA to the CPU, and then, executes i B calculation.In the same time interval, the current i B is moved from the CPU to the FPGA to finally execute the Aurora WR operation.Aurora WR lasts till the next time interval from t = 3T Sim to t = 4T Sim , where RTDS rack 1 reads i B current and passes to the CSC its value, concluding the interval by passing through i B to the PSCC.Finally, RTDS rack 1 calculates v A in the last time interval from t = 4T Sim to t = 5T Sim .The ITM IA results in a latency of 4T Sim by considering from the first v A calculation in the RTDS rack 1 at t = T Sim to the effect of i B to v A calculation at t = 5T Sim .

4) OPAL-RTDS:
The source circuit in the RT-LAB draft and load circuit in the RSCAD model of the first two cases are used in this second hybrid configuration.The physical interconnection of the optical fiber link has been changed from port SFP00 of OPAL rack 1 to port 24 of RTDS rack 2. The ITM IA latency results 5T Sim for all T Sim values.As in previous cases, contributions are described by presenting the configuration's sequence operation diagram in Fig. 5(d).From t = 0 to t = T Sim , OPAL rack 1 executes the same operations described in Section IV-A2.However, RTDS rack 2 reads v A value from the Aurora link in the next time interval from t = T Sim to t = 2T Sim .From t = T Sim to t = 5T Sim , the RTDS rack 2 executes the same operations of Section IV-A1.So, the OPAL rack 1 retrieves the i B value at the end of the time interval from t = 4T Sim to t = 5T Sim .In the last interval, the CPU in the OPAL rack 1 receives the FPGA i B value, calculates the v A effect, and finally, presents its value as a result.The overall calculated latency is 5T Sim , considering the first v A calculation in the OPAL rack 1 to the effect of i B to v A calculation at t = 6T Sim .

B. ITM IA Time-Domain Accuracy Analysis
The time-domain analysis demonstrates that our infrastructure obtains good accuracy of the cosimulated solution.We compared the results of our distributed cosimulation system with a monolithic solution that runs the standalone circuit in Fig. 2(a).The results in this subsection are achieved only for a time step duration T Sim = 500 μs (i.e., worst-case scenario).This scenario is considered as a limit case since normal EMT analysis for electric grid scenario exploits a time step duration of around 50 μs.The standalone electric system in Fig. 2(a) runs simultaneously to the ITM scenario to fetch the standalone voltage and current, namely v real A and i real A .Results in Fig. 6 are presented only for voltages v real A (blu line), v A (green line), and v B (orange line) since currents are in phase, as the power factor of a purely resistive circuit is 1.Voltages v real A , v A , and v B have been analyzed for the two Z B values that represent the stable ITM case (i.e., Z B = 500 Ω) and near the instability region of the ITM circuit (i.e., Z B = 50.5Ω).Finally, the instability transient generated by case Z B = 50.5Ω is presented in Fig. 7.The results demonstrate that applying a T Sim lower than 500 μs ensures a good time accuracy for both regions, reproducing with high fidelity the monolithic solution.

(a). v A overlies v real
A confirming that the ITM application is capable of reproducing correctly the standalone simulation with a slight rise of 2.28% of v A voltage peak caused by the round trip latency of the ITM circuit.v B follows v A , delayed of 1500 μs that is equal to 3T Sim .This is T D 1 latency of the ITM open-loop function G ol in (1).This latency is described in Section IV-A1 and is composed of the first three time step contributions.
The case Z B = 50.5Ω instead is presented in Fig. 6(b) and presents major instabilities due to the phase shift generated by the ITM application.In particular, the distortions are generated by the round trip latency 5T Sim equal to 2500 μs and the Z A /Z B ratio equal to 0.9900, near the instability region of the ITM open-loop function G ol in (1).v A initial peak exceed 40.72% compared to v real A .v B follows the v A trend with a latency of 1500 μs.Similarly to the stable ITM case, it results in 3T Sim .In Fig. 7(a), the distortion transient lasts 0.4 s eventually stabilizing with a 7.92% rise compared to the voltage rise of the case Z B = 500 Ω.Finally, v A presents a nonlinear distortion of the sine due to the ITM application near the instability region.

(c). v A overlies v real
A with an insignificant rise of 0.37% compared to the v A voltage peak caused by the smallest round trip latency between the different infrastructure configurations equal to 2T Sim .v B follows v A , delayed of 500 μs that is equal to 1T Sim .This latency is described in Section IV-A2 and is composed of the first time step contribution.The case Z B = 50.5Ω instead is presented in Fig. 6(d) and still presents major instabilities.In this case, the distortions are denser than in Section IV-B1 due to the higher frequency of the phase shift generated by the round trip latency 2T Sim and the Z A /Z B ratio.v A initial peak exceeds 26.56% compared to v real A .v B follows the v A trend with a latency of 500 μs.Similarly to the stable ITM case, it results correctly in 1T Sim .In Fig. 7(b), the distortion transient lasts 0.16 s stabilizing the result with an 1.25% rise compared to v real A .To conclude, v A does not present particular nonlinear distortion.

(e). v A follows v real
A with a 0.83% rise of the v A voltage peak.v B follows v A , delayed of 1000 μs that is equal to 2T Sim confirming the description of the sequence operation diagram in Section IV-A3.The case Z B = 50.5Ω instead is presented in Fig. 6(f) and presents similar instabilities to the previous cases.The distortions generated by the round trip latency 4T Sim and the Z A /Z B ratio produce a v A initial peak with a rise of 16.97% compared to v real A .v B follows the v A trend with a latency of 1000 μs, which results correctly 2T Sim .In Fig. 7(c), the distortion transient lasts 0.32 s stabilizing the result with a 5.14% rise compared to v real A .Finally, v A presents nonlinear distortion smaller than the case presented in Section IV-B1 due to the lower latency experimented by the ITM IA circuit.

4) OPAL-RTDS:
The case Z B = 500 Ω is presented in Fig. 6(g).This case presents results similar to the case described in Section IV-B1.However, v A overlies v real A with a smaller rise equal to 1.47%.v B follows v A , delayed of 1500 μs that is equal to 3T Sim .The case Z B = 50.5Ω instead is presented in Fig. 6(h) and presents the same instabilities as the case presented in Section IV-B1.In particular, v A initial peak exceeds 40.72% compared to v real A .v B follows v A with the latency of 1500 μs, confirming the 3T Sim latency expectation.In Fig. 7(d), the distortion transient lasts 0.4 s stabilizing the result with a 7.92% rise compared to v real A .Finally, v A presents the same nonlinear distortion as the case presented in Section IV-B1.500 Ω) and near the instability region (i.e., Z B = 50.5Ω) of the proposed infrastructure.For the sake of simplicity, only the case RTDS-OPAL is presented, which presents an intermediate latency of 4T Sim .Results are presented in Fig. 8 for time-domain accuracy in the three transient conditions, namely, the voltage source amplitude transient, the load impedance transient, and the voltage source phase shift transient.

C. ITM IA Time-Domain Transient Accuracy Analysis
1) Voltage Source Amplitude Transient: The voltage source v S is modified after 0.5 s from the nominal value 100 to 90 kV.In the stability region in Fig. 8

(a), v A follows v real
A with an 1.47% rise of the v A voltage peak.Near the instability region in Fig. 8(b), v A instead presents similar instabilities to the initial transient.In the transient condition, the distortions generated by the round trip latency 4T Sim and the Z A /Z B ratio produce a v A initial peak with a rise of 5.73% compared to v real A .The distortion transient is negligible and stabilizes instantly the result with a 5.15% rise compared to v real A . 2) Load Impedance Transient: The load impedance Z B is modified after 0.5 s from 200 to 500 Ω in the stability region, and from 50.5 to 60 Ω near the instability region.In the stability region in Fig. 8

(c), v A follows v real
A with an 1.47% rise of the v A voltage peak.Near the instability region in Fig. 8(d

D. ITM IA Frequency-Domain Accuracy Analysis
The frequency-domain analysis allows a complete understanding of the effects of the latency experimented by applying the ITM IA circuit to the proposed cosimulation infrastructure.This analysis is obtained by applying Welch's method for the power spectral density (PSD) estimation applied to voltage signals v real A and v A for both Z B values to compare the monolithic implementation w.r.t. the cosimulated one.As depicted in Section IV-B, the latency experimented from the different configurations of the proposed infrastructure gives rise to nonlinearity in the time-domain results generated by the phase shift caused by the ITM IA application.The highest effect is resulting near the instability region (i.e., Z B = 50.5Ω) for all configurations.
In fact, the phase shift time-domain effect near the instability region is similar to a triangle wave trend summed to the original voltage signal.Applying additive synthesis, the triangle is approximated in time domain summing odd harmonics of the fundamental sine wave of frequency f Δ while multiplying every other odd harmonic by −1 and multiplying the amplitude of the harmonics by 1 over the square of their mode number n described as follows: Because the phase shift effect changes sign for each round trip latency of the open-loop function G ol , we may consider the fundamental sine wave period of the resulting triangular wave T Δ in (1) twice the round trip delay.Then, the fundamental frequency f Δ is equal to the inverse of the period T Δ .So, a spike in f Δ will be always present for the case Z B = 50.5Ω.Following the approximation in the time domain with additive synthesis, the other frequency components of the triangle wave will be the odd harmonics 3f Δ , 5f Δ , and so on.Since the resolution of Welch's method for the PSD estimation is limited by the sampling period T Sim to 1 kHz, some cases will present also odd harmonics depending on the round trip latency.As matter of fact, the higher the round trip latency, the lower is f Δ , and consequently, its odd harmonics.
1) RTDS-RTDS: v A PSD accurately overly v real A peak at f = 50 Hz (i.e., the power supply frequency) for Z B = 500 Ω without any other frequency disturbances, as depicted in Fig. 9(a).So, the standalone frequency content result is correctly replicated.The case Z B = 50.5Ω in Fig. 9(b) instead presents v A PSD with the former peak at f = 50 Hz and three frequency peaks at f = 200, 600, and 1000 Hz.Since the round trip latency calculated in Section IV-A1 is 5T Sim and following the previous assumption, T Δ results 10T Sim and f Δ is exactly 200 Hz that is the fundamental sine wave component of the triangle wave.Consequently, the frequencies of the odd harmonics are 600 Hz, 1000 Hz, and so on.These three frequency components are appreciated in Fig. 9(b).Moreover, they can be noticed also in Fig. 9(a) (i.e.,Z B = 500 Ω) but their effect is mitigated by the magnitude of Z A /Z B equal to 0.1.
2) OPAL-OPAL: v A PSD accurately overly v real A peak at f = 50 Hz (i.e., the power supply frequency) for Z B = 500 Ω in Fig. 9(c).For Z B = 50.5Ω, the following two main components result from the PSD estimation as depicted in Fig. 9(d): 1) f = 50 Hz, which is the former power supply frequency; 2) f Δ = 500 Hz that is the fundamental sine wave of the triangle wave spectrum.The round trip latency of the ITM application for the OPAL-OPAL interconnection is 2T Sim , resulting in a T Δ = 4T Sim that confirms the PSD component of f Δ = 500 Hz.Since the PSD is limited to 1 kHz, the odd harmonics 3f Δ , 5f Δ , and so on, of the triangle wave approximation cannot be appreciated.
3) RTDS-OPAL: The stable ITM application for Z B = 500 Ω correctly reproduce v real A peak at f = 50 Hz in the v A PSD represented in Fig. 9(e).For Z B = 50.5Ω instead two harmonics of f = 250 Hz and f = 750 Hz can be appreciated as well as the former power supply frequency f = 50 Hz, as depicted in Fig. 9(f).Following the additive synthesis hypothesis and considering the round trip latency of the RTDS-OPAL interconnection equal to 4T Sim , T Δ results in 8T Sim and so the calculated fundamental sine frequency of the triangle wave confirms the former peak at f = 250 Hz and its first odd harmonic at f = 750 Hz.
4) OPAL-RTDS: Since for this case, the round trip latency is equal to the case in Section IV-D1 (i.e., 5T Sim ), we can assume similar observations on the frequency-domain accuracy.For Z B = 500 Ω in Fig. 9(g), v A PSD closely reproduces v real A peak at f = 50 Hz without disturbances.For Z B = 50.5Ω, the former power supply peak is reproduced with three harmonic components at f = 200, 600, and 1000 Hz as depicted in Fig. 9(h).With the same assumptions of the case in Section IV-D1, the fundamental sine frequency of the triangle wave is exactly 200 Hz and its visible odd harmonics are the former peaks at f = 600 and 1000 Hz.

V. CONCLUSION
In this work, we proposed a locally distributed digital tealtime cosimulation infrastructure to connect two DRTS by means of Aurora 8B/10B protocol.The infrastructure is capable of reducing the communication latency experienced by the data exchange among DRTS, allowing to run EMT analysis of an SoS cosimulated power system scenario.Moreover, the infrastructure offers IEEE1588 PTP standard as a synchronization method to avoid misalignment of the real-time executions, permitting to compare results coming from every single DRTS for logging and postprocessing purposes.
The infrastructure proposes the PHIL ITM IA as the most simple method to split the PSUT into submodels each one runs by a different interconnected DRTS, following an SoS approach.Similar to a PHIL setup, ITM IA determines a stable and accurate numerical solution of a power system within the following limitation: Z A /Z B 1. Furthermore, using the Aurora protocol helps in bringing down the effect of latency on the ITM IA, and consequently, enhancing its stability and accuracy.The communication latency is calculated exploiting the ICT infrastructure of DRTSs under analysis, highlighting the internal data exchange among CPU and FPGA to correctly manage Aurora 8B/10B communication protocol.The simple PSUT has been run in a scenario with a time step duration of 500 μs (i.e., worst case) to evaluate the accuracy in the time domain of the solution in both regions Z A /Z B 1 and Z A /Z B ≈ 1 of the ITM IA.The proposed infrastructure results are in both regions acceptable in accuracy and reproduce the performance of the standalone electric system.For instance, the amplitude accuracy of the stable cases is around 1.4% with a negligible transient usually shorter than one cycle.In the unstable region, the amplitude accuracy instead suffers from transients that can last up to 27 cycles, depending on the specific case, with a steady state over peak error of almost 5%.Since EMT analysis normally runs with smaller time steps, around 50 μs, we conclude that our infrastructure allows the EMT analysis of a larger scenario, enhancing the scalability of PSUT.It is worth noting that a smaller time step would relax the constraint Z A /Z B 1, allowing to operate in Z A /Z B ≈ 1 region.Future works will involve the application of a real smart grid scenario to the proposed infrastructure and comparing its accuracy on the system dynamics with a standalone simulation.

Fig. 1 .
Fig. 1.Hybrid locally distributed digital real-time cosimulation infrastructure and its three main architectural layers.

Fig. 2 .
Fig. 2. Standalone electric system (a) consisting of an alternate current voltage source u 1 and two impedances Z A and Z B split via (b) SoS approach by exploiting ITM IA.

Fig. 4 .
Fig. 4. Sequence operation diagrams of the FPGA Aurora RD (red blocks) and WR (green blocks) implementations in (a) RTDS NovaCor and (b) OPAL-RT OP5700, highlighting the data exchange between the FPGA and CPU operations (orange blocks).
), v A instead presents similar instabilities to the initial transient but is reduced by the fact that we are moving far away from the instability region.In the transient condition, the distortions generated by the round trip latency 4T Sim and the Z A /Z B ratio produce a v A initial peak with a rise of 7.86% compared to v real A .The distortion transient lasts one cycle stabilizing the result with a 5.05% rise compared to v real A . 3) Voltage Source Phase Shift Transient: In this case, a 90 • phase shift is introduced in the voltage source v S at 0.5 s from the nominal value.In the stability region in Fig.8(e), v A follows v real A with an instantaneous transient peak that increases 6.15% with respect to the nominal value.After the transient in the steadystate condition, v A follows v real A with a 1.43% rise of the v A voltage peak.Near the instability region in Fig.8(f), v A instead presents similar instabilities to the initial transient.During the transient condition, the distortions are generated by the round trip latency 4T Sim and the impedance ratio of the system, i.e., the stiffness.Mathematically, this is Z B /Z A and in this case is approx 1.0, i.e., not stiff.Hence, changes to the load or source magnitude or angle have a significant impact on the resulting voltage.As a result, v A shows an initial peak with a rise of 157.18% compared to v real A .The distortion transient lasts almost 27 cycles, before stabilizing to a new result with an 5.14% rise compared to v real A .

Fig. 8 .
Fig. 8.Comparison of voltage signals in time-domain for the standalone electric system (orange line) vs ITM IA transients (blue line) for T s = 500 µs in RTDS-OPAL configuration in the region Z A /Z B 1 (a), (c), (e) and near the region Z A /Z B ≈ 1 (b), (d), (f), with different imposed conditions, namely the voltage source amplitude transient (a), (b), the load impedance transient (c), (d), and the voltage source phase shift transient (e), (f).