Tung Hoa H Nguyen

Quoc Lam Nguyen

Tuan T Tran

Yvon Besanger

Raphael Caire

Hoa Van

Quoc Tuan L Nguyen

Senior member, IEEE, Yvon Besanger, Senior member, IEEE Raphael Tran

Senior member, IEEE Caire

V H Nguyen

Q T Tran

Integration of SCADA Services and Power-Hardware-in-the-Loop Technique in Cross-Infrastructure Holistic Tests of Cyber-Physical Energy Systems

Keywords: Interoperability, SCADA-as-a-service, Crossinfrastructure, Hardware-in-the-loop, Co-simulation, Real-time simulation, Holistic Testing

published or not. The documents may come

I. INTRODUCTION

I N order to achieve the decarbonisation scenario, the power system is modernized and is embedded with various Information and Communication Technologies (ICT) and automation solutions, including advanced metering infrastructure, big data assessment and cyber-security [START_REF] Steinbrink | Simulation-based validation of smart grids-status quo and future research trends[END_REF]. This juxtaposition of multiple technologies has transformed the power system into a cyber-physical energy system (CPES), featuring a strong interrelation and correlation of the individual domains [START_REF] Van Der Meer | Cyber-Physical Energy Systems Modeling, Test Specification, and Co-Simulation Based Testing[END_REF]. A comprehensive and holistic validation strategy for CPES needs to harmonize distinct behaviours (continuous-discrete event)

between the Informatics and Communication Technologies (ICT) and the energy systems (electrical grid, heating system) and various different time scales, from several µs (electrical transient behaviour) to several days (market management). It requires a suitable complexity of validation environment, which is not always feasible in the framework of a single Research Infrastructure (RI). For large-scale CPES assessment, researchers are coupling multiple interoperable platforms [START_REF] Lehfuss | Coupling of Real-Time and Co-Simulation for the Evaluation of the Large Scale Integration of Electric Vehicles into Intelligent Power Systems[END_REF]- [START_REF] Wang | A Distributed Control Scheme of Microgrids in Energy Internet and Its Multi-Site Implementation[END_REF] into a holistic experiment to allow the mutualisation of equipment (hence, reduction of investment) and the combination of expertise in different domains of CPES.

The classical validation workflow in CPES assessment is based on mainly two approaches: simulation and real hardware testing. Simulation provides the advantages of rapidity, flexibility and versatility while does not risk damaging the equipment. Real hardware testing requires much more time and investments and is hard to reconfigure in case of an adaptation, but it allows the consideration of real behaviour and impact of equipment. These individual methods, up to now, were used in a fragmented manner, allows the consideration of some certain aspects of the big picture with a limited interaction to the others, without any guarantee how the whole system would react to a solution (e.g. a local optimization may not lead to global optimization; how a single device failure may affect the network or how the communication performance influences the system stability, etc.). More advanced approaches to cope with the multi-domain natures of CPES are: co-simulation [START_REF] Steinbrink | Simulation-based validation of smart grids-status quo and future research trends[END_REF], real-time (RT) simulation [START_REF]PES Task Force on Real-Time Simulation of Power and Energy Systems, Real-time Simulation Technologies for Power Systems Design, Testing, and Analysis[END_REF], hardware-in-the-loop (HIL) and power-hardware-in-the-loop (PHIL) [START_REF] Lauss | Characteristics and Design of Power Hardware-in-the-Loop Simulations for Electrical Power Systems[END_REF] or eventually the combination of both techniques [START_REF] Nguyen | Using Power-Hardware-in-the-loop Experiments together with Co-simulation in a holistic approach for cyber physical energy system validation[END_REF]. The coupling of advanced techniques involving hardware (e.g. HIL/PHIL) over a long distance is, however, challenging due to the impact of latency to the accuracy and stability of the experiment, which may eventually damage the equipment under test [START_REF] Nguyen | Real-Time Simulation and Hardware-in-the-Loop Approaches for Integrating Renewable Energy Sources into Smart Grids: Challenges & Actions[END_REF].

An important aspect in the CPES system that is not sufficiently considered in current validation frameworks is the Supervision, Control and Data Acquisition (SCADA) system and its impacts on system performance. In a large scale CPES, the SCADA system can regroup and orchestrate all the ICT and automation solutions. In current CPES validation methods, the SCADA system, in particular, and the telecommunication aspects, in general, are mostly considered by integrating communication simulators [START_REF] Palensky | Co-Simulation of Intelligent Power Systems -Fundamentals, software architecture, numerics, and coupling[END_REF], [START_REF]Interfacing Power System and ICT Simulators: Challenges, State-ofthe-Art and Case Studies[END_REF]. While this solution is cost-effective and allows the consideration of the interrelation of cyber-physical aspects, the SCADA system of CPES is a complex multi-layer system with industrial protocols and information models, which are not always supported by the communication simulators (mostly ICT oriented). Moreover, the realistic telecommunication characteristics such as sampling time of sensors, response time of controller, server capability, transition of different information models and protocol, redundancy of architecture (history and back-up servers) can't be fully included with current communication network simulators. They are, however, important for a correct and realistic assessment of advanced functionalities of CPES (distributed control, cyber-security, big-data analysis, etc.). On the other hand, the real SCADA system is often associated with a fixed physical system and does not provide enough flexibility for testing of CPES in different configurations and environments.

In this paper, we propose an approach for interoperability of SCADA services and RT simulation/HIL/PHIL of (remote) RI in a holistic platform for CPES validation and assessment. The integration of SCADA functionalities to the cross infrastructure holistic tests [START_REF] Nguyen | Integration of SCADA services in crossinfrastructure holistic tests of cyber-physical energy systems[END_REF] via the SCADA-as-a-service architecture [START_REF] Nguyen | SCADA as a service approach for interoperability of micro-grid platforms[END_REF] introduces realistic data and SCADA functionalities in real-time to the CPES validation framework while the integration of RTS/HIL enables the consideration of CPES taking into account multiple domains and time scales, as well as real equipment behaviour. This approach allows the configuration of complex and sophisticated validation framework for CPES using realistic data and functions from the SCADA system in flexible or extreme scenarios managed by RTS/HIL.

The approach is demonstrated via a case-study of voltage and frequency restoration in isolated microgrid, implemented in a coupled infrastructure of PREDIS (Grenoble INP) and PRISMES (CEA Ines) platforms, located 70 km apart; including RT simulation, OPC UA SCADA system and PHIL in a holistic test setup.

The rest of the paper is organized as follows. Section II presents the approach of cross-infrastructure holistic testing for smart grid and the interoperability requirements for such test. The integration between the SCADA system and the RTS/HIL platform (among multiple RI) is considered in section III. In section IV, the proposed method is demonstrated via the crossinfrastructure holistic test-case. Section V concludes the paper and outlines possible future directions.

II. CROSS-INFRASTRUCTURE HOLISTIC TESTING FOR SMART GRID AND INTEROPERABILIY REQUIREMENTS

CPES are multi-domain systems and due to their diverse and complex natures, the implementation of their validation framework is challenging:

• The integrity and the interdependency of the subsystems in different domains need to be well reflected in the validation framework. • A local solution obtained by a single domain test may not be (or compatible with) the global solution. • Expertise from multiple domains is required to fully characterize and to implement a holistic framework. Combining complementary expertise and platforms of RI in a holistic test is, therefore, a judicious solution that does not require heavy investments in time, infrastructure and human resources. To successfully integrate different actors together in a holistic framework (e.g. SCADA system, RT simulation platform and different kinds of equipment), it is necessary and imperative to achieve interoperability among them. When the holistic test involves combining assets from several platforms, an inter-platform interoperability is also required for their connection and integration. Achieving interoperability also facilitates the connection and integration of new platforms to the group and reduces the cost of installation of experimental modules and integration processes.

A. Interoperability requirements

In [START_REF] Nguyen | On Conceptual Structuration and Coupling Methods of Co-Simulation Frameworks in Cyber-Physical Energy System Validation[END_REF], we proposed a five-layer interoperability architecture for a holistic CPES test (Fig. 1). A holistic experiment integrates information and data from various sources with different formats (via cosimulation or via platform coupling). Depending on the level of coupling that one wishes to achieve in the framework, it is necessary to achieve interoperability up to a certain abstraction layer. The practical setup and harmonization between different parties involved in the test need to be taken care of accordingly. ). Interoperability at syntactic and semantic levels allows the correct and automatable processing of the scenario. • Finally, the conceptual scenario level describes the multi-domain system configuration (e.g. interconnection power/ICT/heating, etc.), and the desired test scenario (e.g. generator x starts at time t1, load y trips at time t2, etc.). Interoperability at this level allows the participating parties to execute their parts and to contribute to the global test scenario correctly. In general, the integration of SCADA system to the holistic framework involves mostly the three lower layers. While the practical implementation of interoperability at physical and dynamic levels are quite straightforward (i.e. implementing common communication protocols or employing a protocol conversion interface). It is not the case for the three higher abstraction levels. In general, the syntactic and semantic convention of the power system are achieved via employing a common information model (e.g. CIM (IEC 61970/61968) [START_REF]Application integration at electric utilities -System interfaces for distribution management -Part 14: Multispeak -CIM Harmonization[END_REF], Multispeak [START_REF] Gary | MultiSpeak and IEC 61968 CIM: Moving Toward Interoperability[END_REF], IEC 61850 dynamic data model [START_REF]TC 57 -Power system management and associated information exchange[END_REF] or OPC UA address space [START_REF]OPC Unified Architecture -Part 1: Overview and concepts[END_REF], etc.). The existing information models are however domain oriented and may not fully cover the inter-domain links (e.g. between ICT-Power or between the electrical grid and the heating system). This gap can be fixed via the inclusion of a scenario description language (standard TTCN-3 [START_REF] Nguyen | On Conceptual Structuration and Coupling Methods of Co-Simulation Frameworks in Cyber-Physical Energy System Validation[END_REF], TPLAN [START_REF] Institute | Methods for Testing and Specification (MTS); TPLan: A notation for expressing Test Purposes[END_REF] and the 62559 template from SGAM model [START_REF]Smart Grid Coordination Group, Smart Grid Reference Architecture (SGAM)[END_REF]) which also enables the implementation of interoperability at the conceptual level. For instance, there are automated or semi-automated methods for implementing popular information models (e.g. CIM [START_REF] Le | CIM compliant multiplatform approach for cyberphysical energy system assessment[END_REF]), converting and communicating between different information models (e.g. CIM and OPC UA [START_REF] Rohjans | Cimbat -automated generation of cim-based opc ua-address spaces[END_REF]). The combination of power system's information models and scenario description languages is unfortunately not easy to achieve, let alone to automate the process.

B. Practical Implementation of Interoperability

In case one desires to include the power system SCADA system, as proposed in [START_REF] Nguyen | SCADA as a service approach for interoperability of micro-grid platforms[END_REF], it is imperative for the participating parties to exchange and to comprehend correctly the system description and the scenario (i.e. to reach at least syntactic interoperability). This task is however challenging because the traditional SCADA system does not support the concept of the information model. Data is processed and stored syntactically only with advanced systems (i.e. object-oriented data model). In case that the SCADA system does not support information model implementation, additional human effort is necessary for harmonization (e.g. pre-agreed convention of naming scheme, message format and communication protocol, etc.). On the other hand, the integration of more advanced SCADA system to a holistic framework in the syntactic layer requires potentially a harmonization and adaptation of the SCADA's information model to the one employed in the framework, or vice versa.

In fact, the physical implementation of SCADA system of infrastructure depends strongly on:

• the field-level technical choices of the system operator: communication infrastructure (e.g. Ethernet, power line communication-PLC, Wi-Fi, etc.), protocols (e.g. IEC61850, DLMs, etc.), information models (e.g. CIM, Multispeak, etc.) and environments (e.g. Labview, etc.), • the topological interrelation of its elements (centralized, distributed, hybrid) • and the actual need for application functionalities (availability, time steps, stability, etc.).

The suggestions for implementations at different levels of interoperability are presented in a generic manner because the specific technical choices would depend on a real situation (e.g. the SCADA's information model, the employed protocols, etc.) and maybe not trivial in another situation.

Implementing interoperability including SCADA system can be done via several levels. At the most basic level (i.e. technical interoperability), the communication infrastructure among parties needs to be integrated. We focus on the following implementation at the dynamic and syntactic level which are required for the integration of the SCADA.

1) Interoperability of Dynamic Data: At this level, the participating parties exchanges dynamic data (e.g. measurement, set points, etc.). These data may not be related to another by default. The interrelations among them (i.e. syntactic and semantic levels), which are necessary to execute the holistic test accurately, are achieved via pre-agreed convention (e.g. common naming scheme, system architecture diagram) and may require human intervention.

An important issue in implementing interoperability at the dynamic level is the transition among various communication protocols. CPES hosts a wide range of protocols from field devices (e.g. Modbus, DNP3, ZigBee, etc.) to the internet (e.g. https), via Local Area Network (LAN) or long-distance communication (e.g. LoRa, GSM, 4G, etc.). Each protocol specifies a different message format. The complexity is often degraded by the variety of operating systems (OS -e.g. UNIX, LINUX, etc.), which demand different memory order (e.g. 32-64 bits, big-Endian or little-Endian, etc.).

In general, OPC and OPC UA can be used for the transition from platform-dependent and component-based field level protocol to SCADA system and other interoperable applications. Besides acting as a translator, OPC UA can also bridge the gap between the local network and wide area services [START_REF] Nguyen | Synchronization conditions and Real-time constraints in co-simulation and Hardware-in-the-Loop techniques for Cyber-Physical Energy System assessment[END_REF].

2) Interoperability of Information Models: In general, it is ideal to implement a common information model for all participating parties in the holistic scenario. However, the SCADA system of each partner is issued from its system configuration and architecture, as well as the development environment and it would be costly to switch. Even if a platform does not possess a SCADA system supporting information model, its data is often organized with its own naming and format convention, which may differ from one partner to another. Most of the time, interoperability in a holistic scenario is achieved with harmonization of those "semi" information models.

The adaptation can be done at the syntactic level or semantic level.

• If the SCADA system is integrated as "black box" and no topology information is required, then a simple naming scheme and message format harmonization suffices for syntactic level. • However, if the topology of the considered system matters during the integration, then semantic interoperability is required and that information needs to be extracted and imported. For instance, CIM allows the representation of topological information in the form of CIM/XML (eXtensible MarkUp Language)/RDF (Resource Description Framework). This description can be used to configure the semantic communication in the framework or to set up the intermediate server for information exchange [START_REF] Le | CIM compliant multiplatform approach for cyberphysical energy system assessment[END_REF]. These implementations can be done over the base of interoperability at the communication level.

III. INTEROPERABILITY OF SCADA SERVICES WITH

REAL-TIME SIMULATION AND (P)HIL PLATFORM Interoperability of SCADA services with RTS and (P)HIL in a holistic platform (Fig. 2) bring out several benefits:

• Possibility to consider the scenario with real data from the expected system environment (including uncertainties that are not readily available with a simple simulation model: e.g. PV/wind installation angle, diversity of meteorological condition, number of users in a building, etc.) and advanced SCADA services. • The RTS can complement the physical platform with virtual equipment that is not available. • RT interface between the SCADA system and the RTS allows the replication and the consideration of real equipment in a virtual scenario, which is sometimes extreme or faulty (e.g. for protection or anti-islanding testing), without damaging the real devices. • The (P)HIL capability allows local consideration of real equipment in their deployment setup, interacting with the virtualized CPES (e.g. harmonic injection, anti-islanding capability, etc.).

A. Communication Implementation

The coupling setup between the SCADA system and the RT platform depends on the architecture of the system and the test-scenario. Generally, the RT platform can be coupled to:

• field devices level (via component-dependent protocols), RTU/PLC level (via IEC 61850/Modbus protocol), • OPC level (via OPC or OPC UA protocol)

• or via application level (via UDP/TCP). The actual setup should be chosen as a function of the required services from the SCADA system and with the criterion of minimizing the latency, e.g. if the RT platform requires measurements from a single device (e.g. RTU), it can be coupled to that RTU via Modbus, instead of going through the upper SCADA layers. The choice is also limited by the communication capability of the RT simulator (i.e. Fig. 2: Integration of SCADA system with RTS and PHIL platform.

supported protocols). Several implementation methods can be

• Via a software interface that converts, (dynamic or steady) data from the SCADA side's protocol to the RTS's protocol (e.g. Modbus to https). This approach is easy to implement (dynamic level interoperability), but lacks support for system semantics and requires more efforts when the system becomes large scale. Via a common database (e.g. SQL, Redis, etc.) that serves as an exchange point. This common database can also be configured to arrange the received data to a defined information model (syntactic and semantic level interoperability). • In the particular case when the required SCADA service is provided by an external actor (e.g. PV production forecasting function developed in Matlab providing prediction to the SCADA system). Then the coupling point (RT platform -SCADA) can be set as the point where the SCADA acquires the service (e.g. a .xml file).

B. Synchronization -Latency Compensation

RT simulation, HIL and PHIL require the simulator to satisfy rigorous constraints on timing and determinism, because the "virtual" elements need to keep up with the real-time behaviour, which is particularly challenging due to the capacity of the simulator and due to the potential variability of time steps. Moreover, latency in the Power-hardware-in-the-loop setup may cause the lead to instability and may cause damage to the equipment [START_REF] Nguyen | Using Power-Hardware-in-the-loop Experiments together with Co-simulation in a holistic approach for cyber physical energy system validation[END_REF].

The coupling of the SCADA system and the RTS/HIL platform introduces into the framework the latency between them as well as the limitations of physical equipment (sampling time, response time, etc.). While they do not prevent the functionality of the interface (i.e. information exchange between the actors), the results may consist disturbances (e.g. sample from a later step may arrive before its precedent samples) and may not reflect a realistic scenario (e.g. the emulators running shortened profiles coupled with real devices).

For synchronization, we can choose to set up the RT platform and the SCADA system to exchange with the time steps satisfying the conditions proposed in [START_REF] Nguyen | Synchronization conditions and Real-time constraints in co-simulation and Hardware-in-the-Loop techniques for Cyber-Physical Energy System assessment[END_REF]. This condition ensures that the exchanged value could be assigned to the same sync step of the other actors (i.e. the order of signal is well preserved), while capturing all the variation and providing a "closest to reality" information to the RT/HIL platform. The last criterion to choose the sync step is the requirement of the considered application (e.g. control generally requires much faster sampling time than monitoring).

Finally, several compensation methods can be considered to compensate for the latency between the SCADA system and the RT/HIL platform or between the RI:

• The exchanged signals are "shifted" through timefrequency transformation [START_REF] Guillo-Sansano | Harmonic-by-harmonic time delay compensation method for Phil simulation of low impedance power systems[END_REF]. The advantage of this approach is that it works well with signals with harmonics. • An observer is constructed to estimate the state of the subsystem subjected to data exchange, based on the measurements of previous steps and the space state of the subsystem [START_REF] Cale | Mitigating Communication Delays in Remotely Connected HArdware-in-the-Loop Experiments[END_REF]. • A similar method using linear prediction for latency compensation is registered in [START_REF] Liu | Geographically distributed real-time digital simulations using linear prediction[END_REF]. This method performs linear prediction directly on the exchanged signal and does not require knowledge on the subsystem space state. All of the existing methods share the two shortcomings: firstly, the latency is assumed to be known. The exact estimation of this latency is however problematic, especially in geographically distributed interoperability involving clock synchronization. Secondly, none of the existing methods can handle well signals with micro-variation, which is very common with real measurements. For RMS signal, they all perform reasonably; however, the performance degrades significantly when we go into sub-second scenarios.

In any case, the latency should be compensated only if it is subsequently introduced by the test-setup. When the latency is a part of the desired scenario (e.g. latency between RTU and field devices), then it should be kept as it is for a realistic consideration.

IV. CASE-STUDY: CROSS-INFRASTRUCTURE IMPLEMENTATION OF VOLTAGE AND FREQUENCY RESTORATION IN ISOLATED MICROGRID

In this section, the proposed approach is demonstrated via a case study of advanced voltage and frequency restoration in an isolated microgrid. The case study is implemented on two remote platforms PREDIS-PRISMES (70 km apart). The validation framework comprises the RTS and PHIL platform coupled with the OPC UA SCADA system in PRISMES, the control algorithm and the communication network simulator located in PREDIS platform.

The considered case study involves the restoration of voltage and frequency in a PV supported micro-grid. In this case study, we are interested mainly in several aspects:

• Validation (good functionality) of the control.

• Demonstration of a holistic test-case, where virtual/real equipment can impact each other in an experimental setup closely representing the desired scenario. • Demonstration of our approach to combine resources and technical expertise from different research infrastructures. • Demonstration of the interoperability of SCADA system and RT simulation.

A. Test-Case Description and Cross-Infrastructure Set Up 1) Test-Case Description: We consider in this test-case the microgrid (MG) supported by four PV packs and four battery energy storage systems (BESS), supplying nine loads; as depicted in Fig. 3. The MG operates in the islanded mode without the support from the main grid. The stability of the grid is guaranteed by the cluster of battery energy storage systems (BESSs). Each BESS has a local controller which operates in parallel with the MG central controller (MGCC) (Fig. 4). A typical hierarchical control structure is applied to control the MG under the fluctuation of loads and PVs power. The aim is to maintain the system frequency and voltage at the references value and to ensure the proportional power-sharing among BESSs.

At primary level, the local droop control of each BESS reacts rapidly to balance power between consumers and suppliers against disturbances:

f -f ref = -k P (P -P ref ) (1) V -V ref = -k Q (Q -Q ref ) (2)
where k P and k Q are the droop coefficient; f -f ref and V -V ref are the grid frequency and the voltage deviation respectively; P -P ref and Q -Q ref are the variations in the active and reactive power delivered by the BESSs to compensate such deviations. The primary control may cause the deviation of frequency and voltage even in the steady-state. The secondary control is then activated to bring the frequency and voltage to the

δf = k If (f -f ref ) (3)
δV = k IV (V -V ref ) (4)
where k If and k IV are the control parameters. Fig. 4 illustrates the control diagram of the i th BESS with two control layers: the local control for primary control purpose and the central controller which takes in charge of secondary control. The secondary control receives voltage and frequency signals from a common measurement point. In this work, the secondary control is in a MGCC located in a long distance site and will broadcast control signals to all BESS local controllers.

The time delays when transferring signals between control layers may lead to voltage and current fluctuation [START_REF] Dong | Time-Delay Stability Analysis for Hybrid Energy Storage System With Hierarchical Control in DC Microgrids[END_REF]. A MG with the hierarchical control structure might have weak stability or even instability due to the influence of time delay. The stability analysis considering network delay has been presented theoretically in previous works [START_REF] Jiang | Delay-Dependent Stability for Load Frequency Control With Constant and Time-Varying Delays[END_REF] based on small signal stability, Lyapunov theory and linear matrix inequalities technique. The control gains in the central secondary control need to be tuned to adapt network communication condition. In this work, we investigate the time delay dependency in a more practical way with the proposed platform.

2) Cross-infrastructure Integration of SCADA -PHIL: The test case is implemented in two parts (Fig. 5), combining the resources and expertise of two partners: CEA INES (power system -electrical manipulation) and Grenoble INP (control -communication network simulation). 3) and ( 4). In this case study, we consider the real PV production and load consummation in the MG. For that purpose, we create a link with the SCADA system of CEA INES to get these measurements in real-time and to inject them to the MG simulated in the RT simulator OPAL RT. The smart houses INCAS and the PV packs (each with 20-60 Si-Monocrystalline panels and PV inverters) are physical hardware (Fig. 6) and are synchronized to the test via the OPC UA SCADA Panorama E2 of PRISMES platform in CEA. In term of model, those PV and load are emulated as current sources with synchronized active and reactive power.

To examine the behaviour of the control algorithm (voltagefrequency deviation) to real devices, the PV pack and load at bus 6 is then replace with the real load Cinegia 30kVA and the real PV inverter SMA 17kVA connected with a PSI PV emulator. These real equipments are physically connected to the MG via a PHIL interface (Fig. 7). The PHIL Interface follows Ideal Transformation Method (ITM).

Since bus 6 consists of inverter-based source and load and has a single point of coupling (PCC) with respect to the considered grid, it is necessary to validate the anti-islanding capacity of the PV inverter in case of grid disconnection. IEC 62116 [START_REF] Commission | Utility-interconnected photovoltaic inverters -Test procedure of islanding prevention measures IEC 62116:2014[END_REF] requires the PV inverter to check the loss of the grid and stop supplying the power (i.e. to disconnect within 2s) to the closely matched demand by the loads (with the tolerance of ±5% -10%). The trip time is calculated based on the IAC values in an islanded condition. This scenario is particularly hard to detect for traditional passive based protection, due to the difficulty of setting admissible parameter range. For the validation of anti-islanding capability, the PSI PV emulators are then calibrated to be well within ±5% of the demand by the Cinegia (Fig. 8). The cluster is then "balanced" and there is almost no power exchange with the other parts (I exchange ∼ 0).

The communication between two partners is configurated through a Redis server working as a private cloud database. Communication emulators at Grenoble INP observes the real communication between two site and are capable of emulating artificial communication scenarios between the MGCC and the MG (e.g. latency, noise, congestion, etc.) to investigate the influence of communication network on the control performance.

B. Experimental Results

The parameters and droop coefficients of BESS local controllers are shown in Table I. Three cases will be considered corresponding to three different communication network scenarios. We validate firstly the functionality of the control in the actual setup, and then we investigate the capacity of the control in a more extreme communication scenario (high latency). Finally, communication with noise is then considered. BESS-1 and BESS-2 BESS-3 and BESS-4

k P 1 , k P 2 4e-5 k P 3 , k P 4 2e-5 k Q 1 , k Q 2 2e-3 k Q 3 , k Q 4 1e-3 L f 1 , L f 2 0.4mH L f 3 , L f 4 0.4mH
1) Case 1: Normal condition: The latency between the two sides of the platforms is measured to be consistently around 32 ms, with some occasional peaks of 85 ms. The integral control gains are chosen as a trade-off between tolerable communications delay and convergence speed:

k If = 2, k IV = 5
The measurement of the active power of BESSs, frequency and voltage at the inverter outputs are shown in Fig. 9.

Initially, the local droop controllers regulate the power outputs of BESSs to ensure the stability and keep the system at the steady-state. However, there exist the errors of frequency and voltages in comparing with the normal state. After that, we trigger the secondary controller remotely, the controller gradually compensates the deviation and recovers frequency and voltage to desired values. The controller works correctly to deal with the disturbances of loads and the fluctuation of PV systems.

2) Case 2: Bad Communication Scenario: In this case, we add 500ms to the latency between the MGCC and the local entities to investigate the behaviour of the controller in a bad communication scenario (i.e. longer distance or slower protocols). With the same control parameters in Case 1, when the secondary start, the controller failed to stabilize the MG as illustrated in Fig. 10. The selection of the integrator gain value is no longer correct for this environment network. It demonstrates the huge impact of the communication network to the system operation if it is not taken into account right in the design phase.

At 98.38s the voltage surpasses the limit of 15% raise and does not return to the allowed tolerance in the required time (2s -this time span may vary from one grid code to another). The overvoltage triggered the protection at the physical point coupling and disconnected bus 6. Due to this sudden voltage drop, the PV current supplying the load flows through the PCC for a brief moment before the PV inverter stops. The time taken by the inverter to check the loss of the grid and stop supplying the power to the closely matched demand by the loads is 0.108 s (i.e. time before IAC returns to 0 -Fig. 11), which is well within the limit of 2 s required by IEC 62116.

The controller, therefore, needs to be modified to adapt to the new working condition. The gains are reduced as follows:

k If = 0.5, k IV = 0.5 Fig. 12 displays the results with the updated controller. When the secondary control is activated, the frequency and voltages are restored to the references. The power quality is maintained even with a considerable step change of loads or different operational states of BESSs (charge/discharge).

3) Case 3: Communication With Noise: We consider, in this case, the performance of the control algorithm under communication noise. The parameter is then reconfigured as in Case 1. An additive Gaussian noise (variance σ 2 = 0.01) is added to the communication between the controller and the microgrid. The power sharing among the BESS, the voltage and frequency are then shown in Figure 13. It demonstrates that the control can withstand communication noise to a certain extent.

The considered test case demonstrates a holistic approach to consider a CPES in its complete form with all the interaction and interdependency between cyber and physical parts. This consideration led to the early reconfiguration of the control according to the deployment environment, right in the design phase. Indeed, in this demonstrating test case, we are capable to observe clearly how the communication performance influences the electrical microgrid behavior, triggers anti-islanding protection on physical inverter and how we could use the

V. CONCLUSION

In this paper, we proposed a method to combine SCADA services, simulation and Power-hardware-in-the-loop, in a cross-infrastructure manner, to create a holistic validation framework for cyber-physical energy system solutions. Using a combination of these fragments to create a holistic and realistic environment, we can achieve transient validation on physical components via localized PHIL experiment; coupled with the desired CPES generated in real-time simulation with SCADA service coupling. Moreover, the method is designed in a vision of cross-infrastructure interoperability, which allows the combination of infrastructures (and the expertise to manipulate the infrastructure) from multiple partners in to a large-scale but holistic CPES validation framework. It helps to save time and money (for investment and formation) and more importantly, to close the gap between the design and the industrial deployment of a solution.

In this vision, the advancements of this approach can be summarized as following:

• The validation framework is closer to reality thanks to the synchronization with SCADA service to create digital twins of physical elements (e.g. PV arrays, building, etc.). • The multi-domain nature of CPES and the interaction among them are addressed via co-simulation. • Local and transient behavior of the integrated DRES can be considered in the most realistic condition with real devices (via PHIL interface) and the propagation of its impact at a global level can be considered in the holistic environment.

• Finally, with cross-infrastructure interoperability, expertise and infrastructure from different partners can be combined to examine a large scale holistic CPES. The proposed method improves the realism of the testing environment, enables the possibility of building large scale CPES testing setup, complements the drawbacks of the existing technologies, without requiring heavy investments in time, infrastructure and human resources.

Fig. 1 :

 1 Fig. 1: Interoperability architecture and practical setup requirements.

Fig. 3 :

 3 Fig. 3: The microgrid configuration considered in the test-case.

Fig. 4 :

 4 Fig. 4: Control diagram of a BESS.

Fig. 5 :

 5 Fig. 5: The multi-site platform setup.

Fig. 6 :

 6 Fig. 6: The INCAS houses and the PV packs at PRISMES that were considered in the test-case.

Fig. 7 :

 7 Fig. 7: The PHIL part setup.

Fig. 8 :

 8 Fig.8: PQ Calibrating error for the islanding situation at bus 6.

Fig. 9 :

 9 Fig. 9: Case 1: without additional latency, k If = 2, k IV = 5.

Fig. 10 :

 10 Fig. 10: Case 2: with additional latency, k If = 2, k IV = 5.

Fig. 11 :

 11 Fig. 11: Anti-Islanding activation at bus 6.

Fig. 12 :

 12 Fig. 12: Case 2: with additional latency, k If = 0.5, k IV = 0.5.

Fig. 13 :

 13 Fig. 13: Case 3: Control under communication noise.

Copyright (c) 2020 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This work is partly supported by the European Community's Horizon 2020 Program -under the projects "ERIGrid II" (Grant Agreement No. 870620) and "United Grid" (Grant Agreement No. 773717).


Dinh Nicolas Retière
email: nicolas.retiere@univ-grenoble-alpes.fr

Truc Ha

Jean-Guy Caputo

J.-G Caputo

Spectral graph analysis of the geometry of power flows in transmission networks

Keywords: Spectral graph theory, Laplacian, nodal domain, power system, power flow, vulnerability

Power flows in transmission networks are driven by the structure of the network and the spatial distribution of generators and loads. Understanding the interplay between grid, generators and loads is crucial for efficient and robust planning and management of electrical networks. Using the spectral properties of the graph Laplacian, we show that we can express power flows on the basis of Laplacian eigenvectors and reveal the dominant modes for nodal voltages and branch powers. The bus voltages are dominated by low rank modes. The power in the lines depends on the nodal distribution of generators/loads, the line impedances and the gradient of the eigenvectors across the branches. The most loaded lines and their associated dominant mode are then identified. A modification of the bus powers is finally proposed to better share the load between the lines and hence to decrease the system vulnerability. Results provided for several IEEE transmission test systems prove the relevance of the approach and suggest practical guidelines to improve the operation of power systems.

I. INTRODUCTION

P OWER grids cope with extreme stress due to the in- creasing demand, the financial pressure from electricity market, the growth of interconnections between regions or countries, the rise of distributed energy resources and the spread of digital technologies [?], [?], [?]. Hence, a single failure may lead to a major collapse and the disconnection of millions of customers. Since the 2000s, more than ten blackouts have impacted more than one billion people around the world [?]. They have not only cost millions of US dollars but have also affected other important infrastructure that deliver essential services such as communication, internet, transport, water and emergency services [?]. The analysis and evaluation of vulnerability of power systems have therefore become a key issue for power engineers and researchers [?]. To prevent major system security issue, it is required to reveal vulnerabilities in power grids. But, because of the system's increasing complexity, new methods are necessary to investigate their structural vulnerabilities [?]. These methods are required to localize and rank the elements less susceptible to absorb hazards, due to their high level of stress in normal operation. The factors and inputs potentially dangerous for structural vulnerabilities are also expected to be explained by these approaches.

Power flow calculation is a fundamental tool for the security analysis of power systems. It can be of use to predict security boundaries, to locate weakest elements, to assess loading of lines, transformers or generators, to verify that voltage values respect the autorized limits [?], Since power systems are networks composed of many interacting components, predicting, interpreting, understanding the vulnerabilities are difficult tasks that require many numerical simulations [?]. The coupling of power flow simulations with complex science has been of great benefit for the engineering community to analyze the results and figure out how power flows distribution and vulnerabilities are driven by the complex structure of electrical networks [?], [?]. A number of studies in the literature have already linked the graph topology of a network to its power flows [?]. Many are based on centrality: which nodes or links are most prominent for the flow of the electrical power in the grid. These geometric approaches evaluate the node degree, the betweenness, the clustering coefficient, to characterize a given topology [?], [?], [?], [?]. Unfortunately, in most cases, these approaches fail to predict with accuracy the power distribution in the network and in particular the most loaded lines. Actually, a central vertex for the geometry may not be central for the power flow. The absence in the centrality-based approaches of the distribution of generators and loads is the main reason of their failure to identify the weakest elements. Alternative approaches use spectral graph theory i.e. the study of the spectral properties of the characteristic matrices of the network. An application of spectral graph theory is clustering: the identification of sets of nodes that present similarities. This method is based on the computation of the eigenvector related to the second smallest eigenvalue of the Laplacian or adjacency matrix of the graph [?]. An extension of this procedure was applied to power systems where the network was partitioned into hierarchical clusters, using the k-th first eigenvectors of a graph matrix [?]. Graph partitioning is strongly connected to the dynamic coherency because it allows to identify groups of generators having similar dynamics. The fundamental work [?] draws theoretical links between dynamic coherency and the confinement of eigenvectors to some sets of nodes. Beyond clustering, spectral analysis of the graph Laplacian was applied to understand the links between the grid structure and its synchronization and stability properties [?], [?], [?], [?], [?]. The identification of the key nodes that drive the network's function can also be performed by spectral analysis of the adjacency matrix [?], [?]. Spectral graph theory was also applied to system security. For example, in [?], the nodal domains related to Laplacian eigenvectors provided the localization of potentially dangerous inputs for the system stability. To improve grid robustness, a lowering of the voltage drops was seeked in [?] by maximizing the second lowest eigenvalue of the Laplacian matrix. The second eigenvalue and its related eigenvector was also used to determine the critical edges of a transmission network [?]. Spectral approaches were finally used in some articles to compute the pseudo-inverse of the network's admittance matrix and thus solve the power flow equations. This allowed to easily compute the redistribution of currents when generators, loads or topology change [?], [?].

All the previous works describe the network's structure in terms of oscillatory modes. These modes are constituted by a couple of variables λ and Ψ where the eigenvalue λ defines the wavelength and the eigenvector Ψ gives the shape of the mode. Our first contribution is the comprehensive modeling of the power flows on the basis of such modes. Indeed, the DC approximation of power flow equations yields a linear system for the electrical variables [?], involving the graph Laplacian matrix (i.e. bus admittance matrix) and the vector of powers (generators and loads). The network topology information is contained in the Laplacian matrix while the power flow input/output is in the vector of powers. The Laplacian matrix is symmetric, positive so that its eigenvalues are real and one can choose a basis of orthonormal eigenvectors to write and solve the power flow equations. This was partially done in [?], [?] but we go one step further, developing the spectral expressions of the power flow equations without requiring the calculation of the pseudo-inverse of the Laplacian matrix. In addition, we show the importance of the eigenvalues and eigenvectors characterizing the network. Eigenvalues provide length-like scales: the low eigenvalues correspond to large scales and vice-versa. Eigenvectors give a geometrical global picture of the network at these different scales. Writing and solving the DC load flow in the Laplacian spectral basis allows to separate each wave-length and "hear" the network on each eigenmode of the graph Laplacian, much like a drum [?]. We hence provide a geometric picture of the drivers of the electrical variables. Our second main contribution is the determination of the modes that dominate the nodal voltages and branch powers. We also show that localized eigenvectors can cause a large power flow in some lines. We summarize the roles of modes on the branch powers by a transfer function connecting the power in the lines to the spectral distribution of nodal powers. This transfer function is higher at low wave numbers corresponding to energy transfers on large distances in the network. It can however exhibit some local peaks at high numbers too. Finally, we identify the most loaded lines whose removal is most susceptible to impact the system and propose a way to reduce their number. To show the usefulness of the approach we apply our approach to four IEEE transmission test cases.

The article is organized as follows. A nomenclature of the variables and parameters is provided in Section II. Basic results on spectral graph theory are presented in section III. Section IV shows the solution of the DC approximated power flow in the spectral domain of the Laplacian matrix. We obtain the voltage angles using the eigenvalues and eigenvectors of the graph Laplacian. These angles are then used to compute the power in the lines. Sections V details the spectral results for the nodal voltage angles. Section VI shows how the modes affect line power flows and Section VII provides a spectral view of the vulnerabilities of the networks. Conclusions are given in section VIII.

II. NOMENCLATURE d

Graph diameter

e i i-
i [n × 1] Active power of node i P L [m × 1] Branch power flow vector V Set of nodes of the graph Y L [m × m] Line admittance matrix λ k k-th eigenvalue of the Laplacian matrix Λ[n × n]
Pseudo-inverse of the eigenvalue matrix ω

Edges weight ∇[m × n]

Gradient over the graph edges

Ψ k [n × 1] Eigenvector related to λ k Ψ k (i) i-th component of Ψ k θ[n × 1]
Nodal voltage angle vector

III. BASIC RESULTS ON SPECTRAL GRAPH THEORY

A. Graph Laplacian

A graph G is defined as G = (V, E, ω) where V is the vertex set of size n, E is the edge set of size m along with the functions ω : E → R + for the weights. The weighted adjacency A and degree matrices D of the graph are given by:

A(i, j) = ω ij if i, j ∈ E, 0 otherwise. D(i, j) = n k=1 A(i, k) if j = i, 0 otherwise. (1)
The Laplacian matrix is defined by L = D -A.

Given a function f defined on the nodes, the Laplacian can also be defined by the quadratic form [?]:

f T Lf = i,j∈E ω ij (f (i) -f (j)) 2
(2)

B. Eigenvalues and eigenvectors

Ψ is an eigenvector associated to eigenvalue λ of L, then:

LΨ = λΨ (3)
Since the edge weights are all positive, the Laplacian matrix is a symmetric, positive-semidefinite matrix. Its eigenvalues are thus real and positive and they can be sorted out such that:

λ 1 = 0 ≤ λ 2 ≤ . . . ≤ λ n (4)
The first eigenvalue is 0 because L is singular; there is only one because we assume a connected graph. The smaller eigenvalues are related to the global connectivity of the network whereas the larger ones quantify the local behavior [?]. The eigenvectors associated to the different eigenvalues can be chosen to be orthogonal and form an orthonormal basis of R n . Using the quadratic form (??) of the Laplacian, it is easy to see that the eigenvector Ψ 1 associated to λ 1 = 0 has equal entries. Then since any eigenvector

Ψ k (k > 1) is orthogonal to Ψ 1 , it verifies: n i=1 Ψ k (i) = 0 (5)
Note that Ψ k (i) corresponds to the node i of the network.

C. Nodal domains

A nodal domain N k associated with Ψ k is defined as a maximal connected subgraph of nodes with value of same sign

[?]: i, j ∈ N k ⇔ Ψ k (i)Ψ k (j) ≥ 0 (6)
The nodal domain is positive (resp. negative) if the value is positive (resp. negative). It is said to be strong if the product of eigenvector entries is strictly positive. In the case of weak domains, the value of some nodes of the domain may thus be zero. . For any eigenvalue, there exists an eigenvector that can be decomposed in at most k nodal domains. Basically, these results mean that the network may be decomposed in two domains for the second eigenvalue and at most n domains for the largest eigenvalue. We will see in the next section that nodal domains plays a major role in the distribution of power flows.

IV. SPECTRAL SOLUTION OF DC POWER FLOW

A. DC power flow

In steady-state operation, any AC transmission network made of n nodes and m branches can be described by a set of n algebraic equations. They link active and reactive powers, voltage magnitude and angle. They are called power flow equations [?]. DC approximation is a well known technique to simplify power flow equations of meshed transmission systems [?]. Under this approximation, power flow equations are reduced to a linear matrix equation:

Lθ = P ( 7
)
where P is the vector whose components are the powers at each bus. θ is a vector giving bus voltage angles and L is the DC bus admittance matrix. By construction, it is equal to the Laplacian matrix associated with the underlying weighted graph of the network. Edge weights are the susceptances of the lines.

B. Spectral solving

Instead of solving the DC power flow in the standard basis of the vertex space of the graph, the orthonormal basis consituted by the eigenvectors associated the eigenvalues of Laplacian matrix is used. The standard basis is defined as a set of vectors (e 1 , e 2 , . . . e n ) where e i is a column vector of size n associated to vertex i. Its components are all equal to zero except the i-th entry which is equal to 1. P is expressed in terms of standard basis elements by:

P = n i=1 P i e i ( 8
)
where P i is the active power injected on node i.

Projecting P onto the orthonormal basis Ψ 1 , Ψ 2 , ...Ψ n , we get:

P = n k=1 p k Ψ k (9)
where p k is the projection of P on Ψ k . Developing line by line this expression, we obtain:

P i = n k=1 p k Ψ k (i) (10)
The same projection is used for voltage angles such that: Because power systems are balanced in steady-state operation:

θ = n k=1 o k Ψ k (11)
n i=1 P i = 0 (12)
It comes:

n i=1 n k=1 p k Ψ k (i) = 0 (13)
This can be rewritten, permutating the sums:

p 1 n i=1 Ψ 1 (i) + n k=2 p k n i=1 Ψ k (i) = 0 (14)
The second term of the right hand side is zero because of formula (??). Since all the entries Ψ 1 (i) are equal and non zero, we get p 1 = 0. Now, we come back to the DC power flow equation (??) and express the vectors onto the eigenvector basis. This yields:

n k=1 p k Ψ k = L n k=1 o k Ψ k ( 15
)
Due to eigenvector definition, it can be expressed as:

n k=1 p k Ψ k = n k=1 o k λ k Ψ k (16)
Then, since the eigenvectors are orthogonal, we get:

p k = λ k o k (17)
The computation of phase angle coefficients is straightforward. For all k > 1, λ k = 0 and:

o k = p k λ k (18)
The component o 1 is undetermined because the system is balanced (p 1 = 0) and the first eigenvalue is zero. Actually, the term o 1 is the reference angle of the network; we choose it equal to zero for simplicity.

To summarize, we have reduced the solution of the full [n× n] system described by (??), to the evaluation of n independent spectral responses (??). This fundamental simplification results from the orthogonality of the eigenvectors Ψ k .

C. Calculation of the flow on lines

The power flow on the transmission lines of the grid can be expressed as:

P L = Y L ∇θ (19)
where P L is a vector made of branch power flows, Y L is a diagonal matrix representing the line admittances and ∇ is the gradient over the edges of the graph. Given a function f on the vertices, it is defined by:

∇f (e) = f (i) -f (j), e = (i, j) ∈ E (20)
It is worth noting that the gradient can be also defined as the transpose of the link-node incidence matrix of the graph [?].

In terms of spectral variables, the power on line i is expressed by:

P L (i) = Y L (i, i) n k=1 p k λ k ∇Ψ k (i) (21)
Apart from the eigenvalue λ k , the line admittance Y L (i, i) and the power distribution on modes p k , the gradient of eigenvector components over the edges is a key variable for power flow on lines. It is related to the structure of nodal domains.

Inside a domain, the edges connect vertices that have with the same sign. Then the gradients over the edges are expected to be smaller than for edges at the border between two adjacent domains. To justify this, recall that a nodal domains is the maximal subgraph that gather all the vertices with values of same sign. Then, two adjacent domains are of opposite sign and at their border, the gradient over the edges could be large.

V. SPECTRAL ANALYSIS OF VOLTAGE ANGLES

Four test systems are studied in this section. They are all meshed and their number of nodes vary from 14 to 300. Network data and reference results are provided by Matpower free Matlab package [?].

A. Eigenvalues and eigenvectors of the test cases

The eigenvalues of the four test cases are arranged in the increasing order and plotted in the log-log diagram of figure ??. Eigenvalues λ 1 that are equal to zero are not plotted.

For homogeneous graphs (no weigths on the edges), the spectral distribution is bounded from above. We have [?]:

λ n ≤ max(D(i, i) + D(j, j)), (i, j) ∈ E (22)
where D(i, i) is the degree of node i.

The algebraic connectivity λ 2 verifies [?]:

4 nd ≤ λ 2 ( 23
)
where d is the diameter of the graph, meaning the number of edges in the shortest path between the most distant vertices. For such homogeneous systems, as the number of vertices is increased, the spectral distribution extends towards zero. It is bounded from above if the degrees of the nodes do not vary much. Electrical networks are more complex because of weights. We see from figure ?? that the spectral distributions extend towards zero as the number of buses increases. Its upper bound increases slowly with n, perhaps due to the weights.

The nodal domain decomposition for the algebraic connectivity of the four test cases are plotted in figure ??. The domains are indicated by the color and the shape of the nodes markers. The size of the edges is also proportional to the weight (admittance value). As predicted by the theory, we see two nodal domains for each network. However, their structure differs significantly due to the geometry and connectivity of the network. On one hand, for the 118-bus case, the admittance values are similar for all the lines and the network is separated into two balanced domains. On the other hand, for the 300bus case, the positive domain is confined to a fork which is connected to the rest of network by a line of very high admittance.

The graph clustering provided by the second smallest eigenvalue is a first approach to study system vulnerability. It provides the minimum number of branches that should be removed to cut the graph in two independent parts [?], [?]

B. Power flow spectral solution

Spectral calculation of power flows is performed by using equations (??), (??) and (??). The results have been compared with DC power flow computation performed by Matpower. Results are very consistent and confirm that this spectral approach is a way to calculate power flows under DC assumptions. The drawback that lies behind its apparent simplicity is the computation cost of eigenvalues and eigenvectors that can be expensive for very large systems. And the spectral solution is only valid for very simple loads and generators models. Its main advantage is to provide a deep insight on the links between graph geometry and transmitted powers. The absolute values of the spectral components p k and o k are plotted along the mode rank k. The power is distributed over all modes k but modes below 15 dominate the voltages which appear to be controlled by power transfers at longer distances. We will see below that the higher modes are determinant for the loading of the branches.

VI. SPECTRAL ANALYSIS OF THE POWERS ON LINES

Our main purpose is now to understand how the graph Laplacian modes drive the powers on the lines. Equation (??) shows that the powers on the lines depend on the eigenvalues λ k , the gradient of eigenvectors over edges ∇Ψ k (and so the structure of nodal domains), the line admittances values Y L Fig. 6: 14-bus -Network line admittances -Line width proportional to admittance (Ω -1 ) and the location of sources and loads which is parametrized by p k . They are the components of the projection of the nodal power vector P on the eigenbasis as shown in equation (??).

Two cases are now examined. In this section, the single mode excitation, meaning that the total load power is only on one mode, is studied. In the next section, the nominal excitation, where the power vector is defined in each Matpower test case, will be detailed.

A. Single mode excitation

We first study the power on lines when the total power of the network is distributed on a single mode k 0 . This means that all the components p k are equal to 0 except p k0 . The nodal power vector is thus expressed by:

P = p k0 Ψ k0 (24)
A node i is generating power if Ψ k0 (i) is positive and is a load node if Ψ k0 (i) is negative.

The power on line i is then given by:

P L (i) = Y L (i, i) p k0 λ k0 ∇Ψ k0 (i) (25)
The case k 0 = 2 is shown in figure ?? for the 14-bus case. The values of powers flowing through the lines are indicated near the edges. They are normalized by the total load power fed into the network. The line width of the edges is also proportionnal to these values. The values of the power flows result from the combination of the eigenvector gradient and the admittance values as shown in figures ?? and ??. The gradient values drive the power flow on lines 5-6. A high power flow is also reached on lines 7-9 because of its relatively low impedance and high gradient. It is worth mentionning that the vertices 4 and 5 are very well connected (very high admittance value) but belong to the same domain (very low eigenvector gradient). The power flow between these two nodes is thus relatively low.

The same analysis can be carried out for the other modes. A summarized view of the results is presented in figure ?? for the four test cases. The maximal value reached by powers on lines in every mode is plotted. The values are again normalized by the total load power in order to compare the four test networks. A general trend can be observed. The highest values tend to be reached in the lower and higher modes. According to the theory, for lower eigenvalues, the network is decomposed into very few domains. The loads are concentrated in the negative nodal domains and the generators are located in the positive domains. The interconnection lines between these domains may be very loaded due to this uneven location of power elements. For higher eigenvalues, the network may be fragmented into a lot of small-size domains (whose number is at most equal to the rank of the eigenvalue). Some lines might be locally sollicitated to transfer a high amount of power between two neighbouring domains.

An illustration of this phenomenon is given for the 30-bus case. In figure ??, the power flow distribution is shown for λ 29 . The i-th component of the eigenvector related to this eigenvalue is affected to the i-th node. Most of the eigenvector entries are zero except for two neighbouring vertices. At the node 22 with positive value, all the generation is connected whereas the total load is connected at the node 21. Of course, the line between these two vertices is very loaded. The power flow and eigenvector structure is more distributed for the eigenvalue λ 30 as shown in figure ??. This is the reason why the maximum power flow is not as high as in the 29-th mode.

A summarized view of the eigenvectors entries is given for the 118-and 300-bus test cases. For the 118-bus case and lower eigenvalues, the eigenvectors are uniformly distributed over the vertices (see the right-hand side of figure ??). But for higher ranks, the eigenvectors appear to be very localized (high values on few nodes). The structure of the 300-bus eigenvectors is not similar (cf. figure ??). Lower modes are surprisingly localized. This result has already been observed for the nodal domain of λ 2 on figure ??. Its negative domain is restricted to a very few number of vertices. In any case, the "excitation" of localized modes shall be avoided to not cause a local excess of power flows on some lines.

B. The network as a transfer function

Coming back to equation (??), the line power flow can be expressed under a matrix form such as:

P L = Y L ∇ΨΛ -1 p (26)
Where p = [p 1 , p 2 ..., p n ] T and Λ -1 is a diagonal matrix with diagonal entries equal to the inverse of the non-zero eigenvalues and zero for the first eigenvalue.

The network is therefore described by a transfer function given by the matrix Y L ∇ΨΛ -1 . This matrix expresses the link between the input vector p (the spectral distribution of generators and loads) and the output response (the power on lines P L ). Each entry of the transfer function describes the contribution of the mode k on the line i. They are plotted in figure ??. They exhibit high peaks and deep valleys for lower rank modes and are more flattened with some localized peaks or valleys for higher ranks. For lower ranks, the transfer function of a transmission network is large to facilitate the transfer of power over the system. For higher ranks, at local scales, the transfer function is close to zero. This means that generally small amounts of power are transferred locally between nodes. But the excitation of the localized peaks and valleys shall be avoided to limit the loads on lines.

VII. SPECTRAL ANALYSIS OF VULNERABILITIES

A. Determination of the most loaded lines

We now investigate the spectral features of the four nominal test cases and assess how and why some lines are more loaded than others. Nominal means that we use the nodal distribution of generators and loads as defined in the Matpower's test cases. The spectra p = (p 1 , p 2 , . . . p n ) of these nominal distributions are shown in figure ??. They are the input of the network's transfer function in figure ??. Notice that these inputs are more concentrated on the higher than the lower modes.

The absolute values of the resulting line powers are plotted in figure ??. They are normalized by the total load value. Of course, the lines are not evenly loaded. However, note the larger the network, the less loaded are the lines. Indeed, larger meshed transmission networks offer more paths for the power flow, therefore increasing the transmission capacity and limiting the risks of line congestions.

In the nominal cases, all the modes are excited and interplay in the power line calculation as expressed by equations (??) or (??). For example, the most loaded transmission element of the 30-bus nominal test case is the line 22 as shown in figure ??. Its spectral decomposition is given in figure ??. The most dominant mode is associated with λ 10 . This results from the combination of the level of excitation of the 10-th mode (figure ??), the gradient value across the line 22 between the nodes 12 and 13 (figure ??) and the admittance of the line (its value is closed to the mean value of the line admittances of the test case). None of these values are very high but their combination give way to an overloading of line 22 when compared to the other lines.

The previous detailed analysis can be carried out for every test cases. The five most loaded lines and their associated dominant mode are given in table ??. These lines are the

L (i, i) 1 λ k ∇Ψ k (i).
least susceptible to absorb hazards due to their high level of stress in normal operation. The results clearly show that they are dominated by high rank modes, i.e. controlled by power transfers at shorter distances. These modes may also control more than one overloaded line. For a given graph structure, the choice of the spectral distribution of generators and loads parametrized by p k is therefore a crucial variable. It should be adjusted to better share the loads between lines and increase the system robustness.

B. Reduction of the vulnerabilities

An example of reduction of vulnerability is given for the 300-bus test case. As previously depicted, the spectral distribution of nodal powers in the nominal case produces very stressed lines. We have then decided to choose another input vector p in order to offer a more uniform sharing between nominal case and the standard deviation is reduced meaning a more homogeneous sharing. Finally, the maximal load has been divided by more than 6 showing that the over-excitation of unwanted localized modes has been avoided.

VIII. CONCLUSION

We introduced a spectral approach to decompose and solve the DC power flow equations on the basis of the Laplacian eigenmodes whose characteristic length is given by the eigenvalues and geometrical shape is determined by the eigenvectors. This provides a geometric picture of the drivers of the electrical variables.

We find that bus voltages are mainly dominated by large scale modes corresponding to the low ranking eigenvalues. Line flows are affected by the nodal domains associated with the signed structure of the components of the eigenvectors. Another important parameter is the series admittance of the lines, for example low impedance lines may magnify some Finally, we introduce a transfer function connecting the power in the lines P L to the generator/load distribution p. This transfer function is large for low wave numbers corresponding to energy transfers on large length scales of the network. High numbers exhibit some localized peaks.

A practical consequence is that electrical engineers designing power systems should avoid high ranking eigenmodes that are very localized because these can cause overloaded lines (table ??). An example of the modification of the power distribution was suggested to better share the loads between the lines and therefore to reduce the system vulnerability.

The challenge open now for power engineers is to select the optimal power input taking into account the spectral structure of the network and the operational constraints.

Fig. 1 :

 1 Fig.1: Sorted eigenvalues 14-bus ( ) -30-bus (+) -118-bus (•) -300-bus ( * )

  Fig. 2: λ 2 -nodal domain -: negative nodes -•: positive nodes.

Fig. 3 :

 3 Fig. 3: 30-bus -Power flow spectral results |o k | (yellow) -|p k | (blue)

  Fig. 4: 14-bus -Power on lines if nodal power is distributed on mode 2 only (p k = 0 for k = 2) -Line width proportional to carried power

Fig. 7 :

 7 Fig. 7: Maximal power on lines.

Fig. 8 :Fig. 9 :

 89 Fig. 8: 30-bus -Line power distribution -Line width proportional to carried power.

Fig. 10 :

 10 Fig. 10: Network transfer function -x-axis: k ; y-axis: i ; z-axis: YL (i, i) 1 λ k ∇Ψ k (i).

Fig. 11 :

 11 Fig. 11: Spectra p = (p 1 , p 2 , . . . p n ) of nodal powers for the nominal cases.

Fig. 12 :

 12 Fig. 12: Power on lines for the nominal cases.

Fig. 13 : 22 Fig. 14 :

 132214 Fig. 13: 30-bus -Spectra of power on line 22

Fig. 15 :

 15 Fig. 15: 300-bus -Histograms of nominal (blue) and modified (yellow) input vectors p

TABLE I :

 I Most loaded lines (sorted by descending order) and associated dominant mode

	300-bus	118-bus	30-bus	14-bus
	Line	Mode	Line	Mode	Line	Mode	Line	Mode
	216	166	13	15	22	10	1	7
	12	218	15	15	22	14	2	7
	296	173	9	18	22	14	3	7
	8	7	147	100	32	29	7	7
	221	218	81	71	5	14	4	7

modes and thus avoid over-excitation of localized modes. The histograms of both nominal and new input vectors are shown in figure ??. This clearly indicates that the modified vector is more homogeneously distributed over the modes. This yields a better distribution of the powers on lines as shown in table ??. The mean value of the lines loading is lower than in the

TABLE II :

 II Statistics on the power on lines (values in MW)

		Mean	Standard deviation	Maximum
	Nominal case	135	187	1292
	Modified case	84	34	173

This work benefited from the support of the project FRACTAL GRID ANR-15-CE05-007-01 of the French National Research Agency (ANR).