A Three-Level Distributed Architecture for the Real-Time Monitoring of Modern Power Systems

To monitor network operation in real-time, power system operators have developed wide-area monitoring systems (WAMS). However, the centralized communication and information processing architecture of WAMS cannot be extended easily to distribution networks. In this aspect, a three-level distributed network monitoring architecture is proposed in this paper, concerning the dynamic analysis of transmission, primary and secondary distribution networks by exploiting measurements of ambient data and transient responses. In the proposed architecture, operators are responsible for the operation and analysis of their own grid but also can share an overview of the system performance to facilitate their operational coordination. Different online and offline applications are supported within the architecture, including small-signal, transient and frequency stability analysis as well as dynamic equivalencing and real-time inertia estimation. Measurement-based algorithms and models are proposed for each case. Finally, the performance of the developed algorithms has been tested by using a combined transmission and distribution power system model.


I. INTRODUCTION
Today, power systems performance and operation has been significantly influenced by the emergence of renewable energy sources (RESs), introducing new challenges and unforeseen problems at system-wide or local grid level [1]- [3]. At the transmission network (TN) an important issue is the intermittent nature of RESs, resulting into high levels of system state variability as well as to uncertainty in realtime operation, that may cause abnormal frequency deviations and stability problems [4]. Additionally, the complex The associate editor coordinating the review of this manuscript and approving it for publication was Peter Palensky . and fast-response of converter interfaced generators (CIGs) of RES plants have introduced new types of stability issues [5]- [8]. Respectively, distribution networks (DNs) that traditionally considered as passive extensions of the TN operational activity [9], are progressively transformed into active components of the bulk system. Indeed, the advancement of distributed renewable energy sources (DRESs) and the need to actively coordinate their functions and operational activities has led to the formation of microgrids (MGs) and the concept of active distribution networks (ADNs). In addition to MGs and ADNs, the changing nature of system loads, due to higher participation of converter interfaced components, has a growing impact on the overall power system dynamics VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ introducing new stability concerns and challenges for the operation of TNs and DNs [8]- [15]. This evolving scenery necessitates the revision of the traditional operational framework of power systems to consider more sophisticated schemes, e.g., ''cascading responsibilities'' [15]. Within such schemes, transmission and distribution system operators (TSOs and DSOs, respectively) monitor, analyze and control in real-time the operation of their own grids. When needed, system operators exchange information and cooperate with each other to analyze potential implications and/or interactions at the transmission and distribution (T&D) interface and coordinate further their actions [1], [8], [11], [15]. One approach is to develop combined T&D synthetic models to investigate the complex operation of TNs and DNs as well as their interactions [16]. However, such models lack usability beyond practical limits, due to their significant computational burden, difficulties to collect the required data and keep the developed models updated. An additional issue concerning availability of data is that sharing of information may be limited due to confidentiality concerns.
To overcome these weaknesses, the advent of smart grid technologies has enabled the development of online monitoring systems supported by measurement-based techniques to analyze directly from measured data the system dynamics [3], [12], [17]- [23]. Measurement-based techniques allow the close to real-time identification of grid characteristics, enhancing drastically the operational monitoring and control capabilities of power systems [24]. Additionally, measurement-based models can be used to facilitate the information exchange and coordinated operation of TSOs and DSOs as well as to overcome possible confidentiality issues [3], [10].
Within this frame, the concept of wide area measurement/monitoring systems (WAMS) has come forward since the late 1980s as a prominent technology to improve the visibility and situational awareness in both today's and future electrical grids. In WAMS synchronized measurements from several network buses are recorded at high data rates (10 to 100 frames per second for 50 Hz grids) from phasor measurement units (PMUs). They are transmitted through communication networks and gathered at a phasor data concentrator (PDC) or other aggregating platforms, where measurement-based techniques can be applied [25], [26]. In contrast to traditional electric monitoring systems, e.g., Supervisory Control and Data Acquisition (SCADA), the accurate timestamping and fast data rates enable the dynamic performance analysis of the system. However, the conventional WAMS architecture cannot be implemented easily in DNs due to the high manufacturing and installation cost of PMUs, the availability of the associated communication and data storage infrastructure, as well as due to the large number of PMUs required for the monitoring of the extended DNs [27], [28]. Therefore, specialized architectures for the distribution level have been introduced, based on different distribution PMU (D-PMU) concepts [29]- [32], e.g., the frequency disturbance recorder (FDR) [29], the micro-phasor measurement unit (µPMU) [30] and the field-programmable gate array (FPGA) PMU [31]. Apart from these implementations, the synchrophasor capability can be also found nowadays into devices such as relays and digital fault recorders. In the near future it is expected to be embedded into inverters, smart transformers, etc., extending further the WAMS applicability and providing a broad and evolving set of uses [26].
Motivated by the advent of WAMS in power grids and the need to develop measurement-based tools and techniques for TSOs and DSOs, the scope of this paper is to introduce a holistic framework for the dynamic analysis of bulk TNs and DNs. The proposed framework is based on a threelevel monitoring architecture of the transmission, primary and secondary DNs, where measurements are exploited to support a number of online (small-signal, transient stability, inertia estimation and frequency response) and offline (dynamic equivalencing and historical data analysis) applications. All-important techniques for measurement-based analysis are incorporated, namely: a) data logging and event triggering to collect ambient data and detect transient responses, respectively, b) a rule-based algorithm to classify the measured data (ambient/transient) under the most suitable application as well as c) signal processing to improve the quality of the measured data. For transient response analysis measurementbased indexes are used to quantify the severity of the events and facilitate the visualization of data. For the estimation of the system parameters (dynamic equivalent model derivation, mode parameter identification and inertia estimation) a unified autoregressive-moving-average with exogenous inputs (ARMAX)-based methodology by using synchronized ambient data and transient responses is proposed. The contributions of the proposed architecture are summarized in the following aspects: • Three-level architecture: a robust framework supported by several tools to analyze power system dynamics at both the TN and the DN level is introduced.
• Small-signal analysis: an automatic machine-learning multi-signal analysis method for the analysis of TN and DNs is proposed; the core of the methodology is ARMAX modelling.
• Transient stability: most conventional methods focus on the analysis of generator signals; the proposed framework utilizes also new indices to evaluate the system stability in terms of network-oriented analysis by using voltage measurements. This is an important and physically meaningful aspect considering the increasing need to investigate ADN dynamics as well as due to the continuous involvement of CIGs.
• Frequency stability and inertia estimation: the system inertia is evaluated by applying the unified ARMAX modelling approach to both ambient data and transient responses; transient responses are also used to compute frequency stability indices.
• Dynamic equivalencing: the unified ARMAX modelling approach is applied and evaluated to ambient data and transient responses to derive equivalent models of both passive DNs and ADNs.
• Simulation results: investigations of power system dynamics at all grid levels by using a combined T&D network model are performed. The rest of the paper is organized as follows: In Section II the framework of the three-level architecture is introduced. In Section III the data logging/event detection procedures, signal classification and requirements per application are discussed. In Section IV the required background, mathematical formulation of the modelling and analysis procedures as well as the signal processing techniques are described. In Sections V and VI the supported online and offline applications are presented, respectively. The simulation model, the operating scenarios and the obtained results are discussed in Section VII. Finally, Section VIII summarizes the most important remarks of this work.

II. THREE-LEVEL ARCHITECTURE OVERVIEW
Employing a three-level scheme, the proposed framework is designed to process, utilize and archive measurements obtained at the primary DN, secondary DN and TN.
The structure of the proposed architecture is illustrated in Fig. 1. Measuring units widely distributed through all parts of the power system acquire voltage and current waveform data. Since measuring units are equipped with global positioning system (GPS) timestamps, waveforms can be referenced to a perfect 50 Hz (or 60 Hz) signal, thus the magnitude and angle of the corresponding phasor representation is generated [1].
Considering DNs, traditionally measurements at the primary and secondary distribution substations are utilized to derive dynamic equivalent models of network feeders (see Fig. 1). Nevertheless, with the advent of modern synchrophasor technologies, data recordings from different locations of DN are also possible. In particular, considering primary DNs, data recordings from D-PMUs as well as by key DRESs can be utilized; for secondary DNs data can be collected from D-PMUs, measuring units in MGs, DRES CIGs, etc. Moreover, in the bulk TN synchrophasor measurements are deployed. In fact, installed PMUs at transmission substations (T-PMUs) or power plants are utilized; additionally, rotor angle data (if needed) can be also provided by the dedicated measuring infrastructure of power plants. VOLUME 10, 2022 With measuring units installed at multiple locations, data are continuously streamed to a central monitoring unit (CMU) associated with the corresponding grid level, i.e., DN or TN (see Fig. 1) for storing in database and real-time monitoring. Within this context, the proposed architecture supports the analysis, operation and operational coordination of both DSOs and TSOs. Each system operator can process the data received from the upstream or downstream grid, merge them with its own and forward them. This way, DSOs and the TSO share an overview of the network situation, remaining at the same time responsible for maintaining their own grid [15].
Leveraging the proposed architecture, a wide range of potential offline and online applications can be supported to facilitate analysis of power system dynamics. In particular, offline applications include statistical analysis and dynamic equivalent model derivation of DNs (both of the primary and secondary levels). Online analysis consists of several realtime applications, namely, small-signal, transient and frequency (including inertia estimation) stability analysis. In the following sections, the online and offline applications are described, according to the flowchart of Fig. 2.

III. SIGNAL LOGGING/EVENT-DETECTION, CLASSIFICATION AND REQUIREMENTS PER APPLICATION
Ambient and transient measured data of frequency (f ), bus voltage magnitude (V ) and angle (θ), real/reactive power (P/Q) and generator rotor angle (δ) are used for the online and offline applications of the three-level architecture.  sense, ambient data are continuously available and intrude the least possible into power system analysis. On the other hand, transient responses are more rarely available, occurring only after a system disturbance or excitation; however, they contain higher information density providing more accurate parameter estimates and requiring less complex models. Details on the signal requirements per application are discussed on the basis of transient response and ambient data analysis as follows: Transient response detection and classification: Through the recorded data, transient responses are stored as separate files. According to the North American Reliability Corporation (NERC) the minimum disturbance detection criteria can be set on the basis of frequency, rate of change of frequency (RoCof), df /dt, dV /dt and voltage limit violations [33]. The generated disturbance file may include at least 30 cycles of data with at least two cycles of pre-disturbance data. Utilities can comply with the NERC recommendations; however, they can use more stringent detection criteria or longer disturbance files [34]. The recorded instances are forwarded to a physics-based rule/decision tree for event classification as follows [35]: • Small-signal analysis: originally includes the processing of f , V , P, Q and δ (if available) measurements [12], [18], [20], [24] at all grid levels by applying measurement-based mode identification techniques, when a relatively small disturbance is automatically detected (dV /dt < 0.1 pu). Nevertheless, online mode estimation is also feasible by using responses related to frequency or large-disturbance events. In the latter case, it is important to exclude the influence of the non-linearities during the first swing period and focus the analysis on the oscillatory part, since non-linearities cause significant degradation efficacy of the identification techniques [36].
• Transient stability analysis: considering dV /dt > 0.1 pu, e.g., in case of a fault event, the V , θ, P, Q and δ transient responses [23], [24] at all grid levels can be used for the assessment of the system transient stability.
• Dynamic equivalencing: when a voltage disturbance at a distribution substation is detected, the recorded V , f , P and Q data are used to derive offline network equivalent models. As detection criterion, the rate of change of voltage (dV /dt) can be used, assuming a threshold value in the range of 0.5 % to 2.5 %. This threshold has been selected as the lowest voltage disturbance which can be caused by on-load-tap-changer actions [37], [38]. Although such types of disturbances are ideal for dynamic equivalencing [39], dynamic responses also caused by any type of action resulting to a voltage disturbance at the distribution substation can be considered.
• Frequency response analysis: this type of analysis refers to system frequency stability assessment as well as to inertia estimation [21]. It is automatically triggered when an active power imbalance is detected, i.e., change in both f and P without significant change in V measurements. Detection criteria of df /dt and dP/dt assuming a threshold of 0.1 pu can be used for the examined signals, i.e., f and P [21].
Ambient data logging and classification: If no system event occurs, the dynamic equivalencing, inertia estimation and small-signal analysis applications are performed by using the ambient data of the corresponding signals described above, collected at a given rate. Finally, it should be mentioned that the three online real-time applications concern both TSOs and DSOs, while dynamic equivalencing only DSOs.

IV. ANALYSIS AND MODELLING BACKGROUND
This Section provides a concrete description on the background of power system dynamics, design of measurementbased parameter estimation techniques and signal processing algorithms incorporated in the three-level architecture.

A. POWER SYSTEM DYNAMICS BASICS
The dynamic performance of a power system, subject to small perturbations (assume a linearized system), can be approximately described by using the following set of multiple-input and multiple-output (MIMO) model equationṡ where A ∈ R N n ×N n , B u ∈ R N n ×N r , B e ∈ R N n ×N q , C ∈ R N m ×N n , D u ∈ R N m ×N r and D e ∈ R N m ×N q . y ∈ R N m ×1 is the measured system output, which contains measurement noise µ ∈ R N m ×1 , caused by the effect of instruments, communication channels, recording units, etc.; the system inputs are classified to known (exogenous), u ∈ R N r ×1 , and unknown, e ∈ R N q ×1 . The block structure of the MIMO system is given in Fig. 4. In this block diagram, matrix G d is the dynamic gain matrix representing the network topological changes via switch sw [40]. Measured output signals,ỹ, can be either ambient or transient responses. Ambient data are considered as low-level random perturbations at e and u. The randomness of these perturbations is typically conceptualized as white Gaussian noise [41], resulting into a response coloured by the system dynamics. On the other hand, transient responses are caused by a step or pulse input at u (including also probing signals) or an excitation at sw, representing changes in network topology. In the transient response case, the m-th system response can be approximated by a sum of N damped sinusoids, superimposed with the underlying ambient noise e m (t) aŝ whereŷ m (t) is the estimated response; ω i , σ i is the mode angular frequency and damping factor and A i , ϕ i the associated amplitude and phase [41]. Assuming the signal is sampled at period T s , (3) is written in z-domain [18] where K is the total number of samples.

B. ARMAX MODELLING
The ambient or transient response of the m-th system output can be analyzed by using an ARMAX model, [41]- [43] with a , b and c indicating the autoregressive (AR), exogenous inputs and the moving average (MA) terms of order p a , p b and p c respectively. Note that, a refer to the coefficients of the characteristic equation related with the system modes. Alternatively, the ARMAX model in z-domain can be described in terms of the discrete transfer function as VOLUME 10, 2022 Here, B m (z)/A m (z) corresponds to the deterministic part of the model, describing the system response to a known input signal; C m (z)/A m (z) describes the stochastic part, that represents the impact of the non-measurable effects on the states of the deterministic part. Polynomials A m (z), B m (z), and C m (z) are defined as and their unknown parameters are calculated by means of a prediction error method [44]. The s-domain poles (eigenvalues), of the system, s i , are estimated as and the associated residues are derived by compiling (1). Once the residues have been determined the mode shape for the i-th mode can be calculated by where the reference channel can be a specific channel of interest, e.g., the channel with maximum magnitude for the i-th mode [24]. Additionally to mode damping and frequency information, mode shapes may be also useful to system operators for control action decisions [45]. Converting the ARMAX model into the s-domain, by using the bilinear transform [44], the transfer function of the system, G(s), can be represented by where α and b are the corresponding continuous-time coefficients.

C. SIGNAL PROCESSING
In case of measurements, a series of data processing tasks must be performed to improve the quality of the recordings. Signals recorded at a sampling rate of 10 samples/s (sps) or higher can be considered at the TN. However, for the DN, the sampling rate limit is at least 50 sps [37].
First, to facilitate parameter estimation and to capture the important system dynamics, the analysis window size must be fixed. For ambient data, the identification time window length varies from seconds to minutes and is updated at regular intervals. Considering transient responses, redundant postdisturbance response samples are excluded from the analysis, by applying a sliding window technique [37]. The end of the window length is obtained, when the mean value of the latest acquired dy/dt (y is the associated signal of the corresponding application, see Section III) samples compared to the mean value of the newest ones is less than a predefined threshold.
In the next step, signal preprocessing is applied to improve the measured signal quality. First, responses are normalized to adjust measured values to a notionally common scale [37]. Next, detrending is applied by subtracting the mean or a best-fit line (in the least-squares sense) from the data. Removing possible trends in the signal, e.g., due to load or generation variation in time, enables the analysis to focus on the fluctuations about the trend, facilitating the parameter estimation procedure. To eliminate the effect of high frequency measurement noise, transient signals are passed through a finite impulse response (FIR), zero-phase, low-pass filter (FIR) with cut-off frequency in the range of 5 -10 Hz [46]. Respectively, ambient data are downsampled; besides filtering, downsampling also improves the computational efficiency of the algorithms. Finally, in case of incomplete data sets, recovery is achieved by applying forward and backward AR modelling. Details on all above techniques are discussed in [37].

V. ONLINE APPLICATIONS
In this Section, the online applications supported by the proposed architecture are described. In particular, small-signal, transient and frequency (including inertia estimation) analysis along with the associated techniques are briefly discussed.

A. SMALL-SIGNAL ANALYSIS
Small-signal stability analysis involves the identification of power system oscillatory modes, namely, the inter-area, intraarea, local and intra-plant modes [47]. In case of multiple signals, the procedure of Fig. 5 is applied separately by each operator to its own grid signals, i.e., by the TSO and the DSO. Assuming M signals recorded from the different recording units of a specific grid, e.g., the TN, the following steps are performed: Step 1: ARMAX algorithm is applied in parallel to the M signals [48]. To determine the best-fit model order for each signal and exclude possible spurious modes additional to the dominant ones [46], mode estimation is performed iteratively, assuming an initial model order estimate. The procedure completes when the calculated R 2 = R 2 h − R 2 h−1 between two successive steps is lower than a predefined tolerance value; where,ȳ is the mean value of the measured signalỹ. It should be indicated that, not all signals participate in the same sense in the formation of the oscillation modes. Therefore, the number of modes contained in the measured data of the M signals may differ and consequently, the resulting order of the set of poles and residues. Step 2: The modes of oscillation are characterized by frequencies in the range of 0.05 Hz to 3 Hz [47]. Therefore, at this step, screening of the derived mode estimates is applied to exclude possible artificial components with frequencies higher than 3 Hz and damping factor lower than −10 s −1 .
Step 3: Multi-signal analysis is performed to obtain an optimum value of the dominant system modes. In particular, the derived mode estimates from the multiple signals (W in total) are automatically grouped via an unsupervised machine learning algorithm [49]. Clustering is used with input feature set θ i = f i , σ i ,Ē i , i.e., the mode frequency, damping factor and normalized energyĒ i . The normalized mode energy for the i-th (here, i = 1, . . . , W ) mode estimate is defined as The optimal number of clusters is determined by applying the heuristic elbow rule (or ''knee'' criterion) to the resulting within cluster sum of squares to between cluster variation (WCBCR) curve [50]. In this sense, the resulting frequency and damping factor centroid of each cluster corresponds to the representative parameters of each group.
Step 4: Representative mode estimates are classified to low and high frequency oscillatory modes considering frequencies below and above 1 Hz, respectively [47]. A further classification is performed to distinguish the mode estimates of each class to poorly-and well-damped [51]. Afterwards, the critical mode, λ cr = σ cr ± jω cr , i.e., the mode presenting the lowest damping factor is determined and is used as a direct indicator of the overall system small-signal stability [52]. The mode classification procedure is illustrated in Fig. 5. Besides transient responses, the system dominant modes can be also estimated by using ambient data [41]- [44] in the same sense.

B. TRANSIENT-STABILITY ANALYSIS
A two-step methodology is applied for the online transient stability assessment (TSA) of the power system in terms of both aperiodic/non-periodic and oscillatory transient instability [5].

1) APERIODIC INSTABILITY ANALYSIS (STEP 1)
when a large-disturbance event is detected, online prediction and evaluation analysis for the TSA of the system is performed.
Potential aperiodic instability is predicted online by means of the methodology proposed in [53] or other machine learning [54] and fuzzy-logic algorithms. Once a disturbance is detected, transient performance evaluation is applied, related to the calculation of wide-area measurement-based indices [55]. In this architecture three indices are adopted covering different facets of power system risk analysis. First, the well-known transient stability index (TSI) defined in (13) is used [23]. TSI evaluates the overall system stability by calculating the difference of the rotor angle of two generators following a disturbance.
Here, δ max is the maximum rotor angle deviation between two generators calculated at the same time instant. Assuming as stable operational boundary the 360 • [23], TSI results into negative values for unstable cases and into positive ones when the power system remains stable after the disturbance. However, in most cases, generator rotor angle signals are rarely time-synchronized. Additionally, using the abovementioned approach, measurements from ADNs cannot be utilized to evaluate the transient stability of the system. An appealing alternative to assess system stability is to use the bus voltage angle signal. In this sense, the transient voltage angle (TVA) index is defined as where θ max is the maximum voltage angle difference. Also, in this case the stable operation boundary is 360 • . Instead of analyzing the overall system performance by means of TSI and TVA, the bus specific normalized transient voltage magnitude (TVM) index is introduced in (15). TVM determines the maximum absolute deviation, U n , of the bus voltage after the fault clears with respect to the pre-disturbance voltage, U 0 , where n ∈ N buses refers to the network bus; N buses is the total number of the system buses. High TVM values, e.g., TVM > 0.1 indicate significant voltage deviations and consequently indication that the system is facing a large disturbance which might potentially lead to instability. On the other hand, TVM values close to zero indicate insignificant voltage deviations and thus minor disturbance where the system is likely to retain stable operation.

2) OSCILLATORY INSTABILITY ANALYSIS (STEP 2)
For cases where no aperiodic instability is detected, information on the damping of the resulting oscillations after the disturbance is obtained by using the matrix pencil (MP) method [56]. MP is applied on the measured signals as a complementary tool to investigate possible loss of system synchronism after a few cycles of oscillations by identifying the dominant modes contained in the recorded oscillatory responses.
To identify the poorly/negatively damped oscillatory modes, characterized by initial low energy, MP is applied to consecutive signal sliding windows. As the energy of the poorly/negatively damped mode increases with time, they can be eventually identified at a later time period. Akin to small-signal analysis, the final selection of poorly/negativelyand well-damped modes is achieved by applying k-medoids clustering [51]. Following the application of k-medoids, the generators participating in the identified unstable oscillatory mode can also be defined.

C. FREQUENCY RESPONSE ANALYSIS
Similarly to transient stability analysis, frequency stability is studied by applying a set of indices to the responses of events that cause a variation to the system frequency [5]. Performance indicators are used instead of responses to help operators to better interpret events and to provide a systematic framework for post-mortem analysis. As a general assessment of the frequency event severity the f nadir/zenith and the f post are used. f nadir/fzenith is defined in (16) and refers to the system frequency nadir/zenith during an event where f nadir,n/zenith,n is the frequency deviation at bus n (see Fig. 6) from the system nominal value, f n . f post , is used to quantify the frequency after the disturbance as follows  where, f post,n is the post-disturbance frequency deviation at bus n (see Fig. 7). High values of f nadir/fzenith and f post imply high frequency deviations. Frequency deviations can be also studied on the basis of RoCoF, calculated explicitly for each bus n as RoCoF n = df n dt (18) where df n dt is the power system frequency derivative at bus n. Frequency transient responses can be also processed to estimate online the power system inertia [21]. In case the system inertia is low, TSOs may decide to dispatch frequency response services or synchronous compensators to enhance frequency stability. A second order ARMAX model (n a = n b = n c = 2) is used, assuming that the frequency and the active power responses are the model input and output, respectively. Specifically, ARMAX modelling is applied to a short window length close to the inertial response period (see Fig. 7), to minimize the impact of additional system dynamics, e.g., caused by frequency control systems and inter-area oscillations [57]. The total window length includes pre-and post-disturbance data of length B 1 and B 2 , respectively, as shown in Fig. 7. Typically, B 1 is set to 2 samples and B 2 to a minimum of 6 samples [58].
In transient response analysis, the effect of the system low level perturbations can be neglected, thus only the deterministic part of the ARMAX model in (6) is considered.
To determine H , the resulting second order ARMAX model is reduced to a first-order transfer function, G(s), as shown in (19) Hence, H , is explicitly calculated by ARMAX modelling can be also used to estimate H from ambient data. In this case the impact of additional dynamics is more pronounced [21], [58]. Therefore, a high order ARMAX model ranging from 2 to 28 is required to fit all dynamics, assuming an analysis window of several seconds [21], [58]. Finally, H is determined by computing the impulse response of the identified ARMAX model. Indeed, the impulse response is equal to the inverse of the effective inertia [21].

VI. OFFLINE APPLICATIONS
Besides online analysis, the proposed architecture supports a number of offline applications presented in this Section.

A. DYNAMIC EQUIVALENCING
Dynamic equivalent models describe the real and reactive power response of a system bus following a perturbation that causes a voltage bus disturbance [60], [61]. In terms of measurement-based modelling, the real/reactive power model parameters (θ P /θ Q ) are determined by using measured voltage and power responses as input and output, respectively, to a curve fitting problem.
In the proposed architecture, the dynamic equivalent model structure proposed in [62] is used. Contrary to most conventional models found in the literature, this formulation is suitable to simulate and analyze the complex behavior of modern ADNs, e.g., reverse power flow phenomena [63]. The block diagram representation of this model is depicted in Fig. 8; the corresponding mathematical formulation is given by (21)-(26) where the polynomial coefficients yield λ 1 + λ 2 = 1 and κ 1 + κ 2 = 1. Function y(t) represents either real or reactive power responses and G(s) is a variable-order linear continuous transfer function. Polynomial functions y t (t) and y s (t) determine the transient and steady-state power variation, respectively, describing the nonlinear network behaviour when subjected to a disturbance. If there is no system event, dynamic network modelling is also feasible by means of ambient data analysis. In this case, y t (t) and y s (t) are neglected and the model is simplified to G (s) = Y (s) V (s); thus, G(s) directly relates the system input V (s) to the system output Y (s). It is worth mentioning that under the ambient data modelling premise, (26) can also result by rearranging (6). DN equivalent models derived by using ambient data are mostly applicable to studies involving the time-varying properties of the grid or smalldisturbances [64]. On the other hand, DN equivalent models derived by using transient responses can be used to investigate a wide range of transient events; however, large disturbance events are less frequent to happen. Therefore, both modelling approaches are of equal importance to cover all types of studies [65].
Considering both transient response and ambient data analysis, ARMAX modelling involves the estimation of the G(s) parameters. In order to determine the optimal G(s) order, Step 1 procedure of small-signal analysis is adopted accordingly by considering as y the ARMAX model output.

B. HISTORICAL DATA ANALYSIS
The identified model parameters, performance indicators and processed data are stored in a database for future analysis by system operators. Different types of analyses can be performed by applying simple visualization techniques, statistical analysis (boxplots, cumulative distribution functions and quantile-quantile plots) [23] or advanced artificial intelligence and machine learning algorithms (cluster analysis, artificial neural networks, etc.) [62]. Some of the most important offline functions include: • post-mortem analysis to evaluate an event and obtain concrete conclusions regarding the challenges and its impact. This way, system operators can plan the best course of actions to mitigate reliability impacts of disturbances [29], • identify tendencies within the system and/or over time considering the dynamic behavior of the system, • model revalidation and calibration of power system components, e.g., power plants, substations, etc. instead of performing expensive and less robust staged tests [62], • develop generalized dynamic equivalent models and derived robust model parameters to represent network dynamics over a wide range of operational conditions different to those originally developed, [62], [64]. VOLUME 10, 2022

VII. SIMULATION RESULTS
The distinct features of the proposed architecture are demonstrated by means of simulations performed on the DIgSI-LENT software [66], using the combined T&D network of Fig. 9. As shown, the examined power system consists of a high voltage (HV) TN, as well as a primary medium voltage (MV) and a secondary low voltage (LV) DN For the analysis, two scenarios are examined. In the first scenario, namely, the passive DN scenario, DRESs located at the DN are neglected and a conventional power system with passive DNs is simulated. In the second scenario, all DRESs are considered in order to emulate a modern power system topology with ADNs; this case is referred as the ADN scenario. For both scenarios, ambient and transient responses are generated by means of RMS simulations to derive equivalent models for DN analysis, to estimate oscillatory modes, to calculate inertia constants and evaluate the frequency stability. Additionally, for both scenarios transient stability indices are computed by processing system responses obtained during large disturbances. Ambient data are generated by varying the load (P and Q power) of the TN bus B6 by ±1% of the nominal power every 2 s, assuming Gaussian distribution to imitate the continuous random variations of actual operating conditions. The simulations are performed, considering a rate of 50 sps to simulate realistic PMU data streams. In addition, to replicate measurement errors, the simulated responses are distorted by means of white Gaussian noise, assuming different SNR levels, i.e., 40 dB, 30 dB, and 20 dB. The noisy signals are processed as described in Section IV.C; specific details if needed are provided in the corresponding subsections with regard to the examined case.

A. SMALL-SIGNAL STABILITY
The applicability of the small-signal stability analysis is tested considering a ringdown event, simulated by disconnecting the line between buses B2 and B3 at t = 0.1 s and by reconnecting it after 0.1 s. Dynamic responses from different buses of the TN and the DN are obtained. For the TN and the primary DN, the recorded signals at each bus include  Table 1.
At each signal SNR = 30 dB is assumed and the corresponding modal components are identified by using a 16 th order ARMAX model. Indicative results for signals recorded at the three grids (TN, primary and secondary DN) are depicted in Fig. 10 for the passive DN scenario. Note that, in the figure, the uppercase delta ( ) is used to denote DC component removal from the processed signals. The estimated responses are compared with the corresponding original (simulated) ones; as shown, very accurate approximates of the original signals are derived, resulting for all cases in R 2 values higher than 94.0%. In particular, in Fig. 10a the TN frequency signal at B5 is examined; it contains a low-and a high-frequency mode at 0.22 Hz and 1.05 Hz, respectively. The mode shapes are presented in Fig. 11a and 11b, respectively. The DN signals depicted in Fig. 10b and 10c, i.e., the active power signal at bus MV11 and the voltage signal at bus LV1, contain each a high frequency mode of oscillation at 1.05 Hz; the corresponding mode shapes are presented in Fig. 11c and 11d, respectively. Note that, the mode shapes have been calculated assuming bus B9 as reference for the TN signal and bus MV1 for the DN signals.
The performance of the ARMAX method is further evaluated via comparisons with two other well-known mode estimation algorithms, namely Prony [72] and MP [56]. The resulting mean R 2 for all signals of Table 1 is summarized in Table 2 as well as the mean computational time. Simulations were performed in an Intel CoreTM i7-930 CPU @ 2.8 GHz processor and 24 GB RAM. Note that, for consistency in the comparison a 16 th order is assumed for all methods. It is shown that, ARMAX and MP methods present very high accuracy, while Prony results to significantly lower R 2 values. Although MP presents slightly higher R 2 to ARMAX, major drawbacks of this method are the higher computational burden, being less efficient for online applications and the limited applicability to transient responses, being unable to analyze data under ambient conditions [73]. The multiple signals obtained from the different parts of each grid may result into varying mode estimates; this can be more intense due to the influence of noise. Therefore, it is very difficult for operators to determine the accurate values of the identified modes [49]. This motivates the use of the procedure of Fig. 5 to calculate representative mode estimates via clustering for the TN and the DN, separately. The signals of Table 1 are used and for each signal one or two  oscillatory modes are identified, resulting into 32 and 24 modal estimates for the TN and the DN, respectively. The mode estimates and their calculated normalized energy are used as inputs in the clustering algorithm to obtain representative mode parameters. Results for the passive DN and the ADN scenarios are summarized in Figs. 12 and 13, respectively. The k-means++, the k-medoids and the Fuzzy C-Means have been used. By applying the elbow rule the number of clusters has been determined to four. In Figs. 12 and 13, the associated true system modes are compared with the identified clusters and the corresponding centroids. The true modes have been calculated by applying the eigenvalue analysis tool of the DIgSILENT software using the linearized power system model [66]; the low frequency mode is at 0.075 Hz with damping factor −0.1343 s −1 (critical mode) and the high frequency mode is at 1.22 Hz with damping factor −0.913 s −1 .
By assessing the results, it is shown that two of the resulting clusters can be related to the two true system modes. Satisfactory classification results are provided by the k-medoids and the Fuzzy C-Means algorithms as also the corresponding Silhouette plots in Fig. 13 reveal; results refer to the ADN scenario though similar plots were obtained also for the passive DN. For the k-medoids and the Fuzzy C-Means algorithms most of the points in the clusters have a high silhouette value, VOLUME 10, 2022  indicating that these clusters are well-separated from the neighboring ones. Generally, clusters without negative values or wide fluctuations in the size of the plots are observed. On the contrary the k-means++ does not provide very consistent results. It should be indicated that, the Silhouette plot analysis also verifies the selection of the optimal number of clusters to four. By comparing the TN and the DN mode estimates, it can be inferred that the dominant system modes can be identified by using responses from all grid levels [48]. Note that, both modes are well-damped as also substantiated from Fig. 10 (damped oscillations).
Finally, the performance of the proposed k-medoids approach is compared with two other multi-signal analysis methods. In the first approach, the oscillatory modes identified by using the signals of Table 1 are classified to low-and high-frequency modes according to Step 4 of the proposed small-signal analysis procedure (see Section V-A). For each group of modes, the arithmetic mean of the corresponding mode estimates is computed [24]. In the second approach, instead of using simple arithmetic mean, weighted averaging is applied assuming mode energy as the weighting factor [48]. Indicatively, the resulting damping factor and mode frequency are compared in Table 3 for the passive DN scenario. Results are grouped according to the grid level and the low-and high-frequency classification. Comparing the mode estimates with those obtained from the eigenvalue analysis, it is shown that the proposed clustering approach and the arithmetic mean provide the more accurate calculations. However, the distinct advantage of the proposed method is that provides an automatic procedure to divide the available data without human intervention.

B. TRANSIENT-STABILITY ANALYSIS
In this Section, the transient stability of the power system in terms of aperiodic instability is evaluated via the TSI , TVA, and TVM indices, computed by using dynamic responses acquired from short circuit (SC) events. For the analysis,  three-phase to ground SCs are introduced sequentially at buses B1, B2, B3, B4, B5, B6, B7, B8, B9, and MV1 of the examined power system. The fault resistance is considered in all cases equal to 0 . Concerning the duration of the examined faults, two cases for the clearing time are considered. In the first case, all SCs clear 100 ms after the occurrence of the fault and in the second case, the clearing time is equal to 300 ms.
Results for the TSI are presented in Fig. 15 by means of boxplots. The TSI has been calculated by using signals, δ, of the three SGs. As shown, SCs that last 100 ms do not cause significant deviation of the rotor angles for this system. Indeed, in all these cases, the TSI is higher than zero, indicating that the power system remains stable for all applied faults. On the other hand, and as shown in Fig. 15, only one unstable case is observed (i.e., with negative TSI value) for the case with clearing time 300ms. This case is illustrated with a red cross in Fig. 15 and corresponds to a SC applied to bus B9.
Indicative results for the TVA index are also reported in Fig. 16; TVA has been calculated by using the responses at  the same bus where the faults have been applied. Note that in Fig. 16 results only for SCs that present a clearing time equal to 300 ms are summarized. As shown, the TVA index receives a negative value only for the SC that is applied to bus B9. As discussed in Section V.B, a negative value indicates system instability. For the rest of the examined SCs, the TVA is positive, indicating that the power system retains its stability. This analysis is in complete line with the TSI results, verifying that the TVA index can be also used to assess system stability by analyzing voltage angles, θ.
By using the same responses as for the TVA, results for the TVM index are summarized by means of a heatmap in Fig. 17. In particular, this figure presents the value that the TVM index receives at different system buses, assuming discrete SC locations. As already discussed, SCs applied at buses B1, B2, B3, B4, B5, B6, B7, B8, and MV1 do not affect system stability. For all these cases, the TVM index receives very low values, i.e., values well below the adopted threshold of 0.1 (see Section V-B). On the other hand, as verified by the TSI , when a SC is applied to bus B9, system instability occurs. For this specific case, the TVM index receives for all system buses values higher than the adopted threshold. These high values constitute a strong indicator regarding system instability.
Finally, the voltage angle responses of bus B2 and MV1, shown in Fig. 18a, considering a SC at B4with 300 ms duration are used to analyze the contained oscillatory modes after the fault has been cleared in terms of oscillatory instability analysis. The MP method is applied to consecutive sliding windows of duration 10 s; as time proceeds the window starts at the next sample. In the examined signals a single mode of oscillation is identified at 1.22 Hz by using a 10 th order MP model. In Fig. 18b the estimated mode damping factor is presented as the window slides up to 4 s. Note that, in the x-axis the time refers to the starting time of the window. It is shown that the damping factor varies against time, especially of the MV1 signal but remains negative for all cases; thus, no sustained oscillations or system instability is identified.

C. FREQUENCY RESPONSE ANALYSIS
The main objective of this analysis is to demonstrate inertia time constant calculation in close-to-real-time by using either ambient or transient response data.
Initially, ambient data are generated for 200 s. Frequency and active power responses of SGs G1, G2, and G3 are recorded. Subsequently, data preprocessing is applied to the noisy signals. For this purpose, both f and P responses are converted into per unit (pu) values. Considering active power responses, the rated apparent power of the SGs is used as the corresponding base. Finally, all responses are filtered by means of downsampling to 5 sps. Representative processed f and P responses of G1 are presented in Fig. 19, considering SNR equal to 40 dB.
To estimate inertia constants of individual SGs through ambient data, the following procedure is applied: f and P responses obtained during a 100 s interval are forwarded as inputs to an ARMAX model and a characteristic transfer function is estimated, as in (19). Afterwards, the impulse response of this transfer function is computed, and the inertia estimate is derived [21]. Nevertheless, inertia calculations  obtained from single estimations are uncertain [59]. Therefore, to reduce estimation error, the sliding window concept is utilized. The length of the employed window is 100 s and the refresh rate is 1 s. In total 100 inertia estimates are derived using the available ambient data. The identified inertia constant is the median value of all individual estimates.
To exemplify the procedure, Fig. 20 is used. In this figure, inertia estimates for SG1 are presented across 100 sequential sliding windows. For the analysis, a 6 th order ARMAX model is used. The starting time of the first sliding window is t = 0 s, whereas the ending time is t = 100 s. Starting and ending times for the last sliding window are t = 100 s and t = 200 s, respectively. As shown, individual estimates vary from 3.52 s to 2.26 s. The median value of all individual estimates is 2.71 s. The absolute percentage error between this median value and the actual inertia constant of SG1 is 3.04%.
The impact of the ARMAX model order on inertia estimates is also quantified by means of the absolute percentage error in Fig. 21. As shown, for the examined power system, a 6 th order ARMAX model provides the most accurate estimates for the total aggregated inertia constant H sys . Indeed, the corresponding error is lower than 8.5%. For H 1 and H 3 , errors lower than 4% are observed; for H 2 an error equal to 17.92% is reported. The impact of SNR on inertia estimates is evaluated on Table 4. In this Table, the absolute percentage error is presented, assuming the different SNR levels. Results indicate that the performance of the proposed approach is not influenced significantly by SNR.
To demonstrate further inertia estimation by using transient responses, a frequency event is simulated by increasing 30% the active power of the load connected to TN bus B6; the event occurs at t = 2 s. Frequency and active power responses of all SGs are recorded and representative sets of transient responses used for inertia estimation are presented in Fig. 22.  To estimate inertia constants through transient responses, the length of windows B 1 and B 2 should be specified. For this purpose, recommendations proposed in [58] are used and B 1 is set to 2, i.e., the pre-disturbance window contains 2 data samples, while B 2 is set to 16, i.e., post-disturbance window contains 16 data samples. Using these settings, inertia estimates of Table 5  The generated ambient signals are initially used to evaluate the accuracy of the dynamic equivalencing modelling procedure. The DN dynamic equivalent model is identified by using a 5-min data window (15000 samples). An indicative instance of the estimated real power response of the DN equivalent model is compared with the original response in   functions (CDFs) for the passive DN scenario are plotted in Fig. 24a and 24c for the real and reactive power models, respectively; CDF plots for the ADN scenario are presented in Fig. 24b and 24d. Results are presented for the optimal model order, i.e., 4 and 2 considering real and reactive power responses for the passive DN scenario; for the analysis of the real and reactive power responses of the ADN the optimal model order is 4 for both cases. It is shown that, the derived ARMAX models can accurately simulate the real and reactive power of both scenarios. This is more evident for measurement errors related to SNR levels higher than 20 dB. Moreover, by comparing Fig. 24a with Fig. 24b and Fig. 24c with Fig. 24d, it can be realized that higher R 2 values are generally observed for the passive DN than for the ADN scenario.
In addition, to demonstrate the ARMAX identification process by using transient responses, a voltage disturbance is caused by applying a 6% step-down at the on-load-tapchanger transformer of the DN. Indicative real and reactive power responses for SNR = 30 dB are presented in Fig. 25 for the passive DN and the ADN scenarios. In the same figure, the corresponding simulated responses by applying the exponential model are compared [19]. Note that, ∼70 % of most of the power system operators worldwide still use static load  models to represent DNs in stability studies; the exponential is one of the most widely adopted models of this type of equivalents. However, results reveal that static models cannot predict the transient behavior of passive DNs with motors as well as ADNs, but can only approximate their new steadystate. On the contrary the superiority of the proposed model is evident, showing R 2 higher than 94% for all signals.
The performance of the derived ARMAX models is statistically evaluated by conducting 100 MC simulations for different noise levels. The resulting R 2 CDF plots for all cases are summarized in Fig. 26 and the corresponding estimated poles are analyzed by means of boxplots or scatter plots in Fig. 27. For the real power modelling of the ADN, a pair of complex conjugate poles is required. On the other hand, one real pole is used for the rest of the examined responses, i.e., P and Q modelling of the passive DN and Q modelling of the ADN. The comparison of the CDF plots in Figs. 24 and 26 referring to the use of ambient data and transient responses, respectively, reveals that the accuracy of the derived models for the latter case is higher for all noise levels. In fact, dynamic equivalencing is feasible for SNR levels also lower than 20 dB. From Fig. 27, it can be noticed that the variability of the pole estimates increases as the SNR decreases. Nevertheless, even for low SNR levels, the model parameters vary in a relatively narrow range; this is mostly evident for reactive power modelling.

E. COMPUTATIONAL BURDEN ASSESSMENT
The computational time for all online applications is evaluated in Table 6 by using the transient responses of Section VII (see Table 1). Parallel computing techniques are used in CMUs, e.g., parallel processing for multi-signal analysis. To thoroughly evaluate the performance of the proposed techniques, discrete sampling rates of 10, 50 and 100 sps are examined. The computer characteristics are reported in Section VII.A. In particular, for small-signal analysis the computational time of the TN and DN multi-signal analysis of Fig. 5 is analyzed, considering 16 signals at each grid level (see Table 1). For transient stability analysis, the computational time regarding the TN refers to the calculation of TVA and TVM of 9 signals (Bus 1 -Bus 9) as well as TSI by using 3 signals; for the DN only TVA and TVM are calculated for bus MV1. Frequency stability calculations of f nadir , f post and RoCoF include signals from the three generators and respectively inertia estimation, H sys is calculated from H 1 , H 2 and H 3 . Note that, in all calculations the required signal processing time is also taken into account; DN results refer to the ADN scenario. As expected, the computational time increases with the sampling rate. Nevertheless, it remains considerably low for all cases, verifying the close-to-real time applicability of the proposed architecture. In fact, even lower simulation times can be achieved, by directly implementing the methods in programming language, instead of using commercial packages, e.g., MATLAB.

VIII. CONCLUSION
In this paper, a three-level network monitoring architecture has been introduced for the coordinated dynamic analysis and evaluation performance of the transmission, primary and secondary DNs, by exploiting ambient and transient response measurements. The proposed architecture supports several online and offline applications, including smallsignal, transient and frequency stability as well as system inertia estimation and dynamic equivalencing. All processed data can be stored in a database for future analysis by system operators.
The proposed architecture incorporates a set of measurement-based tools consisting of an event triggering and classification rule-based algorithm, a signal processing procedure and a unified dynamic analysis methodology based on ARMAX modelling. The latter is used for the estimation of the system dominant oscillation modes, mode shapes, system inertia and the development of dynamic equivalent models. Additionally, a set of indices has been developed, consisting of already known or newly introduced in this work, to evaluate the small-signal, transient and frequency stability of the power system. The performance of the proposed architecture has been demonstrated and tested by using simulation responses from a combined T&D network. Further conclusions taken from this research are: • The adopted ARMAX modelling approach can accurately estimate the dominant modes of oscillation and mode shapes by using measured data recorded both at the TN and also at the DN. A multi-channel cluster analysis method has been proposed to calculate representative modal parameters of low frequency and high frequency modes at each grid. Consistent results have been obtained by using the k-medoids and the Fuzzy C-Means algorithms. Compared to other single-signal or multi-signal approaches, the proposed machinelearning ARMAX-based method succeeds a trade-off among accuracy, computational burden and automatic processing.
• Unlike most research on transient stability analysis focusing only on the TN, the proposed architecture has extended the application of transient stability indices to analyze the dynamic behavior of ADNs, allowing the coordinated monitoring and operation of the overall power system.
• ARMAX models can be applied to both ambient and transient response data to estimate the inertia of synchronous generators and of the overall system. In general, more accurate inertia estimates are obtained by using transient responses regardless of the noise level. However, special care should be taken to determine the model order, especially when using ambient data, or the window length for transient responses. Transient frequency responses are also evaluated in terms of three well-known indexes. Using this information, statistical methods or machine learning techniques can be developed to provide insights concerning frequency stability of the power system, e.g., techniques to correlate system inertia with expected RoCoF values, etc.
• Similarly, to inertia estimation, ARMAX modelling can be used to derive dynamic equivalent models of DNs from both ambient data and transient responses; similar remarks concerning the modelling procedure are also deduced for this case. Accurate ARMAX models can be developed by using transient responses regardless the noise level; for ambient data it is also feasible, but for relatively high SNR levels, e.g., higher than 20 dB. A comparison of the results considering passive DN and ADN scenarios has shown that for the former case first-order models can be generally preferred. However, for the latter, higher order may be required, especially for the analysis of the real power responses. In line with this remark, comparisons with the widely adopted exponential model have clarified that static equivalent models should be avoided for dynamic studies of both passive DNs and ADNs.
• The computational burden of all applications has been tested verifying their efficiency for online analysis. In summary, the proposed architecture offers a holistic framework for online and offline analysis of the different aspects of dynamic phenomena that may occur at the different parts of modern power systems. This concept can be used to enhance the accuracy and reliability of dynamic analysis performed by TN and DN power system operators and extend their visibility and coordination. Future extensions of the three-level architecture may include voltage stability analysis, multi-stability boundaries assessment (in terms of small-signal, transient, frequency and voltage stability), investigations of interactions between TNs and DNs, the fine-tuning of control systems, e.g., power system stabilizers, the incorporation of protection algorithms and investigations on signal processing and communication requirements.