Low Computational Complexity Model Reduction of Power Systems With Preservation of Physical Characteristics

A data-driven algorithm recently proposed to solve the problem of model reduction by moment matching is extended to multi-input, multi-output systems. The algorithm is exploited for the model reduction of large-scale interconnected power systems and it offers, simultaneously, a low computational complexity approximation of the moments and the possibility to easily enforce constraints on the reduced order model. This advantage is used to preserve selected slow and poorly damped modes. The preservation of these modes has been shown to be important from a physical point of view and in obtaining an overall good approximation. The problem of the choice of the so-called tangential directions is also analyzed. The algorithm and the resulting reduced order model are validated with the study of the dynamic response of the NETS-NYPS benchmark system (68-Bus, 16-Machine, 5-Area) to multiple fault scenarios.


Low Computational Complexity Model Reduction of Power Systems With Preservation of Physical Characteristics
Giordano Scarciotti, Student Member, IEEE Abstract-A data-driven algorithm recently proposed to solve the problem of model reduction by moment matching is extended to multi-input, multi-output systems.The algorithm is exploited for the model reduction of large-scale interconnected power systems and it offers, simultaneously, a low computational complexity approximation of the moments and the possibility to easily enforce constraints on the reduced order model.This advantage is used to preserve selected slow and poorly damped modes.The preservation of these modes has been shown to be important from a physical point of view and in obtaining an overall good approximation.The problem of the choice of the so-called tangential directions is also analyzed.The algorithm and the resulting reduced order model are validated with the study of the dynamic response of the NETS-NYPS benchmark system (68-Bus, 16-Machine, 5-Area) to multiple fault scenarios.Index Terms-Coherency, dynamic equivalents, model reduction from data, model reduction of power systems.

I. INTRODUCTION
T HE problem of model reduction [1] consists in finding a simplified model of a given complex mathematical description maintaining at the same time some key properties of the original system.Since the mathematical models used to describe power systems can easily reach hundreds of states, inputs and outputs, the simulation of power systems for dynamic analysis, trajectory sensitivity analysis and control design is a computationally intensive task.As the computational power has advanced, the complexity of these mathematical descriptions and the complexity of the simulated scenarios have increased as well.This trend, which is not unique of power systems, has maintained the computational needs at the top or over the available possibilities [2,Section 5.7].Hence, the model reduction problem, sometimes also referred as the problem of dynamic equivalencing [3], is central to modern research on power systems.At the basis of the use of dynamic equivalencing in power systems there is the idea of distinguishing between a study area, the description of which is maintained in full detail, and an external area, consisting of the remaining part of the network, which is reduced.The study area can be precisely analyzed and controlled while still considering the interaction of the interconnection of such area with a less faithful representation of a larger external area.Thus, when doing model reduction of power systems the main aim is having a dynamic behavior of the interconnection between the reduced order model and the study area as close as possible to the one of the original nonlinear power system.In other words, we are interested in the faithful reproduction of the behavior of the system for a specific class of input signals, neglecting the behavior outside the operating conditions.
Historically, coherency-based methods have been used in model reduction of power systems, see, e.g., [4]- [12] and [3], [13] for a comparison of different coherency-based methods.These methods are based on the physical properties of the electrical machines connected to the network.The idea is to find coherent generators, i.e., machines which behave similarly when the same input is applied.Once coherent generators are identified a dynamic equivalent generator is used to replace them.In recent years, the power system community has started to be interested in reduction techniques based on mathematical properties instead of physical ones, see, e.g., [14]- [25].One of the reasons of this interest is the flexibility of having a reduction technique that is not based on the physics of the generators and, as a consequence, the possibility of reducing networks with renewable energy sources.Among these methods, balanced truncation and Krylov projectors have been successfully used in power systems reduction, see, e.g., [3], [26], [27].One of the drawbacks of the techniques based on Krylov projections, also called moment matching methods, is the difficulty in enforcing or preserving specific properties of the system to be reduced.Since the moment matching technique is completely decoupled from the physics of the system, the real applicability of the reduced order models is limited.For instance, one would expect that the preservation of some physical characteristics, such as some specific modes, is of paramount importance.Poorly damped modes, also called electromechanical modes [28], are important in the small-signal stability analysis of a power system since they are responsible for most of the oscillating behavior.Similarly slow modes characterize how fast the response of the system reaches the steady-state.While previous attempts to maintain these modes are essentially ad hoc, since with the classical Krylov methods it is very hard to preserve a certain set of eigenvalues, the problem of assigning a prescribed set of modes has been alleviated in [29].
In a preliminary version of this paper [30], the model reduction techniques presented in [29] have been applied for the first time to power systems.In particular, these methods have been used to assign arbitrary eigenvalues to the reduced order model, i.e., maintaining slow and poorly damped modes.Therein it has been shown that it is not necessary to increase the order of the reduced order model to improve the quality of the approximation of the system: to achieve this goal it may be sufficient to select a different set of eigenvalues to be preserved in the reduced order model.Moreover, in the paper it is also shown that increasing the order of the reduced model gives no guarantee of improving significantly the quality of the approximation if the wrong set of eigenvalues are preserved.The results of [30] give a physical justification to a mathematical deduction, since the resulting reduced order model is better both in terms of approximation error and physical properties.
In the present paper we shift the focus from demonstrating the reasons that justify a new approach to the problem of model reduction of power systems to the actual detailed description of the method.We begin giving a solid theoretical foundation.Power systems are not only described by a large number of states, but they have also a high number of inputs and outputs.Thus, the problem of model reduction of multi-input, multi-output (MIMO) systems is described in the framework of the so-called tangential interpolation problem.An algorithm to approximate the tangential directions is given and connections with previous published results are drawn.Since classical efficient Krylov algorithms suffer the drawback of not preserving specific modes, we present here a low computational complexity algorithm for the model reduction of MIMO systems which is able to preserve specific properties.Originally devised for single-input, singleoutput (SISO) systems for which the mathematical description is unknown [31], the algorithm is here extended to MIMO systems.The algorithm, exploiting the inputs and outputs of the system instead of operations on large matrices, offers a computationally efficient method to approximate reduced order models.Note also that this algorithm lays the foundations to achieve model reduction by moment matching directly from time-domain measurements.This would complement recent results obtained with coherency-based methods in which frequency-domain measurements are used to achieve model reduction, see, e.g., [32], [33] in which phasor measurements units are used to determine reduced order models.The simulation section illustrates the performance of the algorithm, the design of the reduced order model and the dynamic behavior of the interconnection of the study area with the approximated model.
The rest of the paper is organized as follows.In Section II we provide some preliminaries.In particular, in Section II-B we introduce the problem of the reduction of MIMO systems and the role of tangential interpolation.In Section III we derive the algorithm for the approximation of the moments of MIMO systems whereas in Section IV an heuristic to approximate the tangential directions is given.In Section V-A we describe the nonlinear model which represents the study area of the power system and the linear model which represents the external area.In Section V-B a complete algorithmic overview summarizes the method.In Section V-C the NETS-NYPS benchmark system is described.In Section V-D the reduced order model is designed.In Section V-E the use of the algorithm is illustrated and in Section V-F the fault response of the interconnection between the study area and the reduced order model is simulated in multiple scenarios.Notation: We use standard notation.C < 0 denotes the set of complex numbers with negative real part; C 0 denotes the set of complex numbers with zero real part.The symbol I denotes the identity matrix and σ(A) denotes the spectrum of the matrix A ∈ R n ×n .The symbol ||A|| indicates the induced Euclidean matrix norm.The vectorization operator is denoted by vec(A).Given a list of n elements a i , diag(a i ) indicates a diagonal matrix with diagonal elements equal to the a i 's.The superscript denotes the transposition operator and ι indicates the imaginary unit.

A. SISO Systems
Consider a linear SISO continuous-time system described by the equations with state be the associated transfer function and assume that (1) is minimal, i.e., controllable and observable.1 Definition 1: [1, Chapter 11] Let s i ∈ C\σ(A).The 0moment of system (1) The model reduction technique by moment matching is based on the idea of interpolating a certain number of points η k (s i ) on the complex plane: a reduced order model is such that its transfer function (and derivatives of this) takes the same values of the transfer function (and derivatives of this) of system (1) at s i .In [34], [35], a characterization of the moments of system (1) has been given in terms of the solution of a Sylvester equation.The importance of this formulation is that it establishes, through the Sylvester equation, a relation between the moments and the steady-state response of the output of a specific system.These observations, which resulted in several developments in the area of model reduction by moment matching, see, e.g., [36]- [41], are summarized in the following result and illustrated in Fig. 1.
Theorem 1: [29] Consider system (1) and suppose s i ∈ C\ σ(A), for all i = 1, . . ., η.Consider any non-derogatory 2 Then there exists a one-to-one relation between the moments η 0 (s 1 ), . .., η k 1 −1 (s 1 ), . .., η 0 (s η ), . .., η k η −1 (s η ) and r the matrix CΠ, where Π ∈ R n ×ν is the unique solution of the Sylvester equation for any L ∈ R 1×ν such that the pair (L, S) is observable; r the steady-state response, provided σ(A) ⊂ C < 0 , of the output y of the interconnection of system (1) with the system for any L and ω(0) such that the triple (L, S, ω(0)) is minimal.Remark 1: The reduction technique based on this notion of moment consists in the interpolation of the steady-state response of the output of the system: a reduced order model is such that its steady-state response is equal to the steady-state response of the output of the interconnected system (provided it exists).
Remark 2: The method can be applied to systems which are not minimal, but in this case the resulting reduced order models describe only the controllable and observable subsystems.This is not at all a restriction but a sensible feature of any model reduction method.In fact, note that not only any Krylov method reduces the minimal subsystems (since these methods interpolate the transfer function) but unobservable and uncontrollable modes are the first candidates for elimination also for other methods.For instance, the balanced truncation method eliminates the less controllable and less observable modes.Since no controller or observer can be designed for such modes, when we want to economize on the order of the system it is natural to disregard these modes.See [1] for a thorough analysis of the relation between controllability, observability, the problem of model reduction and, more in general, the problem of "realization" (obtaining a state-space model).For a specific analysis of the problem of quantifying controllability and observability in power systems see [42]- [44].

B. MIMO Systems
Consider a linear MIMO continuous-time, system described by the equations C p×m be the associated transfer function and assume that (4) is minimal.If we characterize the moments of this MIMO system as in Definition 1, we have mpν interpolation conditions (the mp components of W at ν points), whereas the dimensions of CΠ is only pν.While a possible solution to this problem is to inflate the order of the reduced order model, this substantially increases the dimension of the model and the reduction may even lose sense if mν ≥ n, which is a likely situation in the analysis of power systems.A workaround to this drawback con-sists in "merging" some of the conditions together.However, the merging of the interpolation conditions has to be done in a specific way to preserve the information of the original single moments.This is the idea behind the tangential interpolation approach, initially proposed in [45].The tangential interpolation approach can be embedded in the framework we have presented with a redefinition of moment.Definition 2: Let s i ∈ C\σ(A) and l i ∈ R m ×1 .The kmoment of system (4) at s i along l i is a vector of p complex numbers defined as Hence, the moment matching conditions become [46], [47] with i = 1, . . ., ν, where Ŵ (s) is the transfer function of the reduced order model.

Remark 3:
The pν tangential interpolation conditions (5) are weaker than the original pmν matching conditions.Since this approach consists in replacing m equations with one, the resulting model is not in general as good as a model which interpolates all the m conditions.This is the price to pay to maintain the order of the reduced order model independent from the number of inputs.Moreover, this suggests that the selection of the vectors l i is fundamental to obtain a reliable reduced order model.
Remark 4: In [26] the authors report the block Krylov approach of [48].This approach, predating the development of the tangential interpolation theory in [45], inflates the dimension of the reduced order model to satisfy the pmν matching conditions.However, the authors recognize this drawback and do not use this method in the application part of the paper.In fact, claiming "in order to limit the size of the Krylov subspaces, we consider that the matrix B [...] is the sum of the input matrices" they are applying a single input algorithm on an approximated system.They also recognize that "in general, such kind of heuristics for economy in the size of the base [...] does not work well".Interestingly, a theoretical explanation for the poor general performance can be given revisiting their approximated approach in the tangential framework.It can be easily shown that adding the columns of the matrix B corresponds to using moment matching with tangential directions l i = [1 . . .1] .Thus, the performance can be improved selecting different directions, i.e., applying the method to a linear combination of the columns of B instead of the simple addition.
Exploiting Definition 2, Theorem 1 can be adapted to the MIMO case.
ν, be such that the pair (L, S) is observable.Then the moments of the system are in one-to-one relation with CΠ, with Π ∈ R n ×ν the unique solution of the Sylvester equation (2).Finally, the family of systems [29] with S − GL ∈ R ν ×ν , G ∈ R ν ×m and CΠ ∈ R p×ν contains all the models of dimension ν interpolating the moments of system (4) at S if G is such that σ(S) ∩ σ(S − GL) = ∅.Hence, we say that system ( 6) is a model of (4) at S. System ( 6) is a reduced order model of system (4) at S if ν < n.Remark 5: All the models that can be obtained using Krylov projectors are encoded in the family of systems (6) [29].Thus the models obtained with the two approaches are equivalent.The advantage of this formulation is that the family of systems ( 6) is parametrized in G, which allows to set with ease several properties of the reduced order model, as shown in [29].For instance, setting the eigenvalues of the reduced order model is an easy task, whereas with the classic Krylov method this is rather difficult.
Remark 6: The family of systems ( 6) represents reduced order models of system (4) for any matrix L, i.e., for any tangential direction.However, as already remarked, the quality of the approximation depends on the choice of L. We also note that the tangential directions depend upon the reduced order model and the reduced order model depends upon the tangential directions.Thus, the determination of the directions is a difficult problem.
Remark 7: In the interpolation framework, algorithms have been proposed to iteratively approximate the vectors l i , see, e.g., [24].However, these solutions do not allow one to chose the interpolation points nor to preserve a prescribed set of eigenvalues.While these algorithms work well when specific physical properties are not important, they are somewhat limited in the case of power systems for which the preservation of slow and poorly damped modes is fundamental for giving a physical meaning to the reduced order models [30].

III. A LOW COMPLEXITY ALGORITHM FOR THE COMPUTATION OF THE MOMENTS
In this section we extended to MIMO systems an algorithm recently proposed for the computation of the moments of a SISO system from input/output data [31], [39].Although the algorithm has been primarily devised to compute the moments when the matrices A, B, C are not available, it has also the advantage of being a computationally fast method for the approximation of the moments of a system.
Remark 8: We do not claim that the algorithm we are presenting here can be used to approximate the moments from actual input/output data of power systems.In fact, for a power system, which is a nonlinear system, we do not have any guarantee on the approximation given by the algorithm when measurements generated by such nonlinear system are used.However, in addition to propose this algorithm for the fast computation of the moments, this paper lays also the foundations to solve the problem of input/output reduction from real time-domain data in two ways: the input/output data of a power system can be filtered by means of a linear filter; the nonlinear version of the algorithm proposed in [40] can be extended to MIMO systems and used with the measured data.However, the use of these methods for the reduction of power systems is not trivial.For the aims of this paper, we limit ourself to use the algorithm just as a mean to approximate the moments without solving equation (2), which is computationally expensive.
In the next theorem 3 we make some technical assumptions for which we clarify the meaning here.Assuming σ(A) ⊂ C < 0 implies that the system is stable.Assuming that the triple (L, S, ω(0)) is minimal implies that it is always possible to find a sequence {t k } of sampling times t i , with lim k →∞ t k = ∞, such that rank ω(t k −ν +1 ) . . .ω(t k ) = ν, see [49].
Theorem 2: Consider the interconnection of system (4) with the signal generator (3).Assume σ(A) ⊂ C < 0 , σ(S) ⊂ C 0 and that the triple (L, S, ω(0)) is minimal.Define the time-snapshots R k ∈ R gν ×ν and Υ j k ∈ R gν as and , where y j (t i ) is the j-th row of y(t i ).If g = 1, the matrix R k has full rank and is an approximation of C j Π, with C j the j-th row of C, i.e., there exists a sequence {t k } such that the equation Then a MIMO version of the algorithm proposed in [31] can be formulated.

Algorithm 1:
Let k be a sufficiently large integer.Select η j > 0, with j = 1, . . ., p sufficiently small.Select g ≥ 1 (recall that gν is the number of rows of R k ).
1: Construct the matrices R k and Υ j k .2: If rank R k = ν then compute C j Π k solving equation (7).Else increase g.If k − gν < 0 then restart the algorithm selecting a larger initial k.
Remark 9: If g = 1 the matrix R k is full rank.However, in practice numerical approximations in the computation of R k may cause loss of rank.For this reason we formulate the algorithm with the possibility of increasing the number of measurements.As a consequence Algorithm III gives a least square approximation of the moments.According to our experience, good values of g can usually be selected as 1 ≤ g ≤ 10 depending on the system.In the simulations presented in this paper g = 4.
Remark 10: For SISO systems, the algorithm has a computational complexity of4 O(gν 2.373 ), whereas the Arnoldi or Lanczos procedures for the approximation of the moments have a computational complexity of O(νn 2 ), with n ν [1, Section 14.1] (both complexities are multiplied by p in the MIMO case).

IV. APPROXIMATION OF L
In this section we present a heuristic that can be used to construct the matrix L of the signal generator (3).
Algorithm 2: Let r = 1 and j = 1.Consider C j the j-th row of C, B r the r-th column of B and l r i the r-th element of l i .
1) Consider the system (A, B = B 1 , C 1 ) and design a reduced order model with l 1 i = 1 for i = 1, . . ., ν, i.e., select the interpolation points and the desired eigenvalues of the reduced order model to achieve the desired approximation of the resulting SISO system.
2) If r < m then r = r 1.Using the same interpolation points and desired eigenvalues compute a reduced order model of the multi-input system

Repeat.
Else with the obtained matrix L, the same interpolation points and desired eigenvalues compute a reduced order model of (A, B, C).

3) Stop.
This heuristic is justified by the observation that L is involved in combining m elements.Thus, we may expect that we 5 can ignore the fact that the system is multi-output in the approximation of L. In addition, one can expect that the set of interpolation points and prescribed eigenvalues which have been selected to obtain a desired approximation on the truncated SISO system can be used on the MIMO system.This is justified by the observation that we still want to maintain an approximation on the truncated system as close as possible to the one we have designed.As a consequence of these observations, the approximation of the mν elements of L has been reduced to the problem of determining m − 1 independent scalars l r i .Remark 11: There is no guarantee that this heuristic works in general.However, our simulations have shown that the resulting L is at least locally optimal with respect the H 2 norm.A small variation of any of the obtained elements of L causes a rapid increase of the error norm.

A. Power System Model
We describe a power system composed of n m -machines and n b -bus with the classical model, see [52], [53], which is normally used in the literature of model reduction of power systems, see, e.g., [3], [26], [27].The model is described by the differential 5 See [51, Chapter 2] for the definition of ) with i = 1, . . ., n m , where δ i and ω i are the rotor angle and angular velocity, respectively, of the i-th machine, ω s is the reference angular velocity, H i and D i are the inertia and damping coefficients, respectively, of the i-th machine, E i is the internal voltage of the i-th machine, Y ij = G ij + ιB ij is the admittance between the machines i and j, G ii is the self-conductance of the i-th machine and T M i is the mechanical input power of the i-th machine.
In the literature on model reduction of power systems the study area and the external area are sometimes modeled as two separate entities interconnected each other with n p -tie-lines, see, e.g., [3], [26].However, this is a somewhat strong approximation.In fact, note that if the two power systems, study area and external area, are interconnected then we have a unique large power system and the power flow analysis which defines the parameters of system (8) has to be updated.Using the tie-lines only to exchange the input and output of the two systems gives the considerable simplification that the number of inputs and outputs corresponds to the arbitrary number of tie-lines.
On the contrary, in this paper the division in study area and external area is a pure exercise of labeling.In fact, the whole power system is described by equations ( 8) and the division in study and external area can be done over every region of the power system using the actual buses as interconnecting lines between the two areas.This approach has the advantage of improving the fidelity of the simulation of the power system.The drawback is that the number of inputs of the external area corresponds to the number of machines of the study area and the number of outputs of the external area corresponds to the number of machines of the external area, and vice-versa.However, since with the considered method the dimension of the reduced order model does not depend upon the number of inputs and outputs, a large number of inputs and outputs is not an issue for the technique we are presenting.
Thus, consider an external areas composed of e m -machines and a study areas composed of s m -machines, with s m + e m = n m .The study area is described by system (8) for i = e m + 1, . . ., s m .The input of the study area (output of the external area) is δ j , with j = 1, . . ., e m .The external area is described by the linearization of system (8) around an equilibrium point, namely where ) is an equilibrium point and the remaining matrices are defined as for i = 1, . . ., e m , j = 1, . . ., e m and j = i, for i = 1, . . ., e m , and for i = 1, . . ., e m and k = e m + 1, . . ., s m .

B. Algorithmic Overview
The steps of the method are summarized in the flowchart in Fig. 2. We start determining a linearized model of the power system.Note that the linear model is (9) if we use the classical model ( 8), but it can be different if we use a more complex original model.We then inspect the frequency response W (ιω) (e.g., the Bode plot) and we order in a list Σ the frequencies s i at which preeminent features, e.g., frequency peaks, appear.We inspect the spectrum of the matrix A and we determine the slowest and poorly damped modes.We order these modes in a list Λ.We generate the data y(t) and ω(t), and we determine CΠ η using Algorithm 1 applied to (A, B = B 1 , C 1 ).We determine the reduced order model ( 6) assigning the eigenvalues from the list Λ, e.g., using the function "pole" of MATLAB.If the SISO reduced order model obtained is not satisfactory, e.g., by comparing the Bode plot, we increase the order of the reduced order model and repeat the previous steps.Otherwise we generate a MIMO reduced order model using Algorithm 2. If the MIMO reduced order model obtained is not satisfactory we increase the order of the reduced order model and repeat the previous steps.Otherwise we stop.

C. New England Test System (NETS)-New York Test System (NYPS)
The theory presented in Sections II and III is validated on the interconnected NETS-NYPS 68-bus, 16-machine, 5-area power system, see [28].The study area, composed of the machines 14, 15 and 16, is interconnected with the bus-lines 18-50, 18-49 and 41-40 to the external area, composed of the machines from 1 to 13.The separation in study and external area follows the geographical distribution of the power system and the tie-lines are actual bus-lines of the power system.The system to be reduced has n = 26, m = 3, p = 13.The parameters of system (8) have been computed in MATLAB as follows.With the script Init_MultiMachine.m, which can be downloaded from [28], the line data, the power flow results, the generator direct axis transient reactance and the inertia values H i of the machines have been generated.Following the theory presented in [52], see also [53], we have computed the reduced bus admittances Y ij , the equilibrium voltages E i and the equilibrium angles δ i .The mechanical input powers T M i have been computed from these quantities and equations (8) written at the equilibrium point.The damping coefficients D i have been generated with the function rand.Finally, the quantities of the linearized system have been computed directly from (10)- (12).

D. Design of the Reduced Order Model
The selection of the interpolation points and of the eigenvalues of the reduced order model has been done exploiting the analysis in [30].Therein it has been shown that good time-domain performance can be obtained with a re-Fig.3. Eigenvalues of the linear system (9) (crosses) and of the reduced order model ( 6) (squares).The dash-dotted lines represent the 10% damping ratio.duced order model of order between 6 and 12 if the correct slow and poorly damped modes are preserved.System (9) has been simulated with the input generated by the signal generator (3).The matrix of the generator has been selected as S = diag(S 0 , 2.12 S, 2.79 S, 3.42 S, 5 S), with S 0 = 0, 1; 0, 0 and S = 0, 1; −1 0 , which corresponds to choosing the interpolation points at 0 (zero and first order moment), 0.3374, 0.444, 0.5443 and 0.7958 Hz (all zero order moments).The reduced order model is computed with L approximated using Algorithm 2 and CΠ approximated with Algorithm 1.The resulting reduced order model has order ν = 10.We determine the six least damped and the four slowest eigenvalues of system (9) and we assign them to the reduced order model (6).In Fig. 3 the eigenvalues of system (9) are represented with crosses, whereas the eigenvalues of the reduced order model are depicted with squares.In the figure the modes in the area between the two dash-dotted lines are well damped (more than 10% damping ratio), whereas the others are considered poorly damped [28].

E. Approximation of the Moments
In this section we illustrate the performance of Algorithm 1 by showing how the approximated reduced order model improves as k in Algorithm 1 increases.The surfaces in Fig. 4 represent the magnitude (left graphs) and phase (right graphs) of the elements W 1,1 (top), W 10,2 (middle), W 11,3 (bottom) of the transfer matrix of the reduced order model as a function of t k , with 3.4448 ≤ t k ≤ 48.4299 s.For comparison, the solid/black lines indicate the magnitude and phase of the respective elements of the transfer matrix of the reduced order model for the exact moments CΠ, i.e., computed solving (2).The figures show how the approximated magnitude and phase of the reduced order model ( 6) evolve over time and approach the respective quantities of the exact reduced order model as t k → ∞.Finally note that we have chosen to show these three particular components of the transfer matrix because the rest of the components are very similar to these.These are the components that present the most distinctive graphical features.Fig. 5 shows the mag-Fig.4. Magnitude of the elements W 1 , 1 (top), W 10, 2 (middle) and W 11, 3 (bottom) of the transfer matrix of system (9) (solid lines) and of the reduced order model (dotted lines).Fig. 5. First scenario -Large: angular velocities of the study area when this is connected to the nonlinear system describing the external area (solid lines) and when it is connected to the reduced order model (dotted lines) with the eigenvalues shown in Fig. 3. Insert: absolute errors between the time histories.nitude of the elements W 1,1 (top), W 10,2 (middle) and W 11,3 (bottom) of the transfer matrix of system (9) (solid lines) and of the reduced order model (dotted lines).We note that the curves of the reduced order model are close to the curves of the system along all the frequencies.As expected, the approximation is not uniformly good for all the elements of W .In fact, the curve in the top graph does not approximate the given system as well as the one of the other two graphs.This is caused by the use of the tangential interpolation, i.e., we are trying to capture more information maintaining the same number of parameters (the order of the reduced order model).

F. Fault Behavior
Dynamic simulations of the power system have been performed.We have considered three fault scenarios.7. Second scenario -Large: angular velocities of the study area when this is connected to the nonlinear system describing the external area (solid lines) and when it is connected to the reduced order model (dotted lines) with the eigenvalues shown in Fig. 3. Insert: absolute errors between the time histories.
1) A self-clearing fault at bus 14 of the study area occurring at t = 1 s and cleared at 1.15 s (clearing time of 150 ms).2) A self-clearing fault at bus 42 of the study area occurring at t = 1 s and cleared at 2.50 s (clearing time of 1500 ms). 3) A self-clearing fault at bus 15 of the study area occurring at t = 1 s and cleared at 1.85 s followed by a self-clearing fault at bus 16 of the study area occurring at t = 3 s and cleared at 3.30 s (clearing time of 850 and 300 ms, respectively).In the first scenario, Fig. 6 shows the angular velocities (large) and respective absolute errors (insert) of the study area when this is connected to the nonlinear system describing the external area (solid lines) and when it is connected to the reduced order model (dotted lines).We see that the two time histories are almost indistinguishable as confirmed by the small value of absolute error.
Remark 12: The dynamic simulation of the fault has been implemented with the function "ode45" of MATLAB, with the option "odeset('RelTol',1e-9,'AbsTol',1e-9)".The elapsed time of the simulation of the nonlinear interconnected power system is 0.550235 s.The elapsed time of the simulation of the reduced order model interconnected with the nonlinear model of the study area is 0.300701 s.Thus, the simulation time of the reduced order model is 54.65% of the original time.It is expected that as the complexity of the original system increases the difference in simulation time will grows (in fact note that the benchmark system is a simple example for which only 13 machines have been reduced).
Since a clearing time of 150 ms may be a too simple scenario to verify that the reduced order model is a good approximation of the power system, we consider now the second and third scenarios.In the second, Fig. 7 shows, with the same color coding of Fig. 6, the angular velocities and absolute errors for a fault cleared in 1500 ms.We notice that, although the angular veloci-Fig.8. Third scenario -Large: angular velocities of the study area when this is connected to the nonlinear system describing the external area (solid lines) and when it is connected to the reduced order model (dotted lines) with the eigenvalues shown in Fig. 3. Insert: absolute errors between the time histories.ties undergo a variation ten times greater than the one in Fig. 6, the system is able to recover and the transient and steady-state behaviors are well approximated by the reduced order model.The third scenario is shown in Fig. 8, with the same color coding of Fig. 6.The second fault is occurring during the transient behavior caused by the first fault and it can be seen as a simulation of a domino effect.Moreover, the simulation targets the 16th generator that in the previous two scenarios underwent a small variation (with respect to the 14th and 15th).Also this more complex scenario is well approximated by the reduced order model.Thus, the analysis and control of the nonlinear system describing the study area can be performed using the interconnection of such area with the reduced order model instead of the full nonlinear description of the external area: in fact, the dynamic response to faults is almost identical but the number of the equations is reduced.

VI. CONCLUSION
We have presented a low complexity algorithm for the fast estimation of the moments of MIMO systems.The estimated moments have been exploited for the model reduction of largescale interconnected power systems.The technique that we have demonstrated offers, simultaneously, a low computational complexity approximation of the moments and the possibility to easily enforce constraints on the reduced order model.This possibility has been used to preserve selected slow and poorly damped modes which are important both from a mathematical and physical point of view.The problem of the choice of the so-called tangential directions has also been studied and an heuristic for their approximation has been given.The techniques have been validated with the study of the dynamic response of the NETS-NYPS benchmark system.

APPENDIX
Proof of Theorem: 2 The matrix Π defined in equation ( 2) is such that ẋ(t)| t k = ΠSω(t k ) + Ae At k (x(0) − Π k ω(0)) (13) hold.Consider the first equation of system (4) computed at t k , namely Substituting (13) in equation ( 14) yields ΠSω(t k ) + Ae At k (x(0) − Π k ω(0)) = Ax(t k ) + BLω(t k ) and multiplying on the left-hand side by C j A −1 we obtain The matrix C j Π k defined in equation ( 7) is such that which yields and

Fig. 1 .
Fig. 1.Diagrammatic illustration of Theorem 1.The first term on the righthand side (solid circled/red) describes the steady-state response of the interconnected system, whereas the second term (dashed circled/green) the transient response.

Fig. 6 .
Fig. 6.The colored mesh represents the magnitude (left graphs) and phase (right graphs) of the elements W 1 , 1 (top), W 10, 2 (middle), W 11, 3 (bottom) of the transfer matrix of the reduced order model as a function of t k , with 3.4448 ≤ t k ≤ 48.4299 s.The solid/black lines indicate the magnitude and phase of the respective elements of the transfer matrix of the reduced order model for the exact moments C Π, i.e., computed solving (2).

Fig.
Fig.7.Second scenario -Large: angular velocities of the study area when this is connected to the nonlinear system describing the external area (solid lines) and when it is connected to the reduced order model (dotted lines) with the eigenvalues shown in Fig.3.Insert: absolute errors between the time histories.