Development and Applicability of Online Passivity Enforced Wide-Band Multi-Port Equivalents For Hybrid Transient Simulation

TL;DR

This paper introduces an online passivity enforced method for developing frequency dependent network equivalents that enable real-time hybrid transient simulations of power systems with unknown parameters.

Contribution

It proposes a novel online recursive least squares algorithm for passivity enforced FDNE development applicable in real-time power system simulations.

Findings

Successfully verified on IEEE 39 and 68 bus power systems.

Achieves accurate, stable, and passive network equivalents in real-time.

Handles unknown network parameters effectively.

Abstract

This paper presents a method for developing single and multi-port frequency dependent network equivalent (FDNE) based on a passivity enforced online recursive least squares identification algorithm, which identifies the input admittance matrix in z-domain. Furthermore, with the proposed architecture, a real-time hybrid model of the reduced power system is developed that integrate transient stability analysis and FDNE. Main advantages of the proposed architecture are, it identifies the FDNE even with unknown network parameters in the frequency range of interest, and yet can be implemented directly due to discrete formulation while maintaining desired accuracy, stability, and passivity conditions. The accuracy and characteristics of the proposed method are verified by implementing on two-area, IEEE 39 and 68 bus power system models.

Tables

Table 1. TABLE I: Coherency grouping of generators for IEEE 39 bus system

Group	I	II	III	IV
Generators	4,5,6,7,9	1,10,8	3	2

Table 2. TABLE II: Comparison of Reduced (EMT+TSA) and Original (EMT) Models

	EMT+TSA(AGG)	EMT+TSA
Relative Speed (Fig. 4)	0.3589	0.2998
Active Power (Fig. 5)	0.0356	0.0351

Table 3. TABLE III: FDNE Formulation Comparison

Type	Proposed Method	Vector Fitting(VF)
FDNE Order	17	21
RMS Error	2.176e-6	2.0101e-6

Table 4. TABLE IV: Comparison of Reduced (EMT+FDNE) and Original (EMT) Models of Two-Area System

	EMT+FDNE(Proposed)	EMT+FDNE(VF)
Relative Speed (Fig. 10)	0.481261	0.4827015
Active Power (Fig. 11)	0.077046	0.0771309

Table 5. TABLE V: Comparison of Reduced Order (FDNE Only) Models

Case	Eigen Value	Frequency(Hz)	Damping(%)
EMT Based Model	-1.323$\pm$i7.544 -0.880$\pm$i5.157	0.6566 0.8327	21.92 16.82
EMT+TSA Based Model(AGG)	-0.906$\pm$i7.096 -1.381$\pm$i6.537	0.5494 1.0634	43.09 20.67
EMT+TSA Based Model	-1.306$\pm$i7.583 -0.598$\pm$i5.014	0.6501 0.8038	5.7543 11.8466
EMT+FDNE Based Model	-0.664$\pm$i6.993 -0.373$\pm$i2.507	0.3822 0.4034	61.26 61.26
EMT+FDNE+TSA Based Model(AGG)	-1.321$\pm$i7.652 -0.670$\pm$i4.984	0.6623 0.8804	4.1457 13.336
EMT+FDNE+TSA Based Model	-1.325$\pm$i7.601 -0.610$\pm$i5.002	0.6518 0.8021	5.2343 12.1216

Table 6. TABLE VI: Comparison of Reduced (EMT+FDNE+TSA) and Original (EMT) Models of Two-Area System

	EMT+FDNE+TSA (AGG)	EMT+FDNE+TSA
Relative Speed (Fig. 13)	0.2286	0.0923
Active Power (Fig. 14)	0.0163	0.0076

Table 7. TABLE VII: Comparison of Reduced (EMT+FDNE+TSA) and Original (EMT) Models of IEEE 39 Bus System

	EMT+FDNE+TSA (AGG)	EMT+FDNE+TSA
Active Power (Fig. 17)	0.1634	0.0489
Voltage (Fig. 18)	0.0112	0.0081
Relative Speed (Fig. 19)	0.0905	0.0433

Table 8. TABLE VIII: Time Required for Generation of Equivalents

Model	TSA(AGG)	TSA	FDNE
Two-Area System	10s	8s	40s
39 Bus System	14s	13s	91s
68 Bus System	19s	17s	97s

Equations

$$L_{index,i}=\sum _{k=1}^n(1-P_{ki}{)}^n$$

$$Y_{red(m\times m)}=\left[Y_{m\times m}-Y_{m\times n}{Y}_{n\times n}^{-1}{Y}_{n\times m}\right]$$

$$\displaystyle\left[\begin{array}[]{c}I_{b(j\times 1)}\\ I_{g(i\times 1)}\end{array}\right]=Y_{red(m\times m)}\left[\begin{array}[]{c}V_{b(j\times 1)}\\ V_{g(i\times 1)}\end{array}\right]$$

$$V_g\angle {\theta }_g=\left(I_g\angle {\delta }_g-Y_{gb}{V}_b\angle {\theta }_b\right)Y_{gg}^{-1}$$

$$I_b\angle {\delta }_b=Y_{bb}{V}_b\angle {\theta }_b+Y_{bg}{V}_g\angle {\theta }_g$$

$$I_{inj}\angle {\beta }_b=I_b\angle {\delta }_b$$

$$relativeerror=\frac{{\Vert y_{ref}-y_{act}\Vert }_2}{{\Vert y_{ref}\Vert }_2}$$

$$Y_f(z^{-1})=\frac{I_F(k)}{V_F(k)}=\frac{b_1{z}^{-1}+b_2{z}^{-2}+...+b_n{z}^{-n}}{1+a_1{z}^{-1}+a_2{z}^{-2}+...+a_n{z}^{-n}}$$

$$Y_{f}(z^{-1})_{m\times m}=\left[\begin{array}[]{cccc}Y_{f(1,p)}+Y_{f(1,1)}+\dots+Y_{f(1,m)}&-Y_{f(1,2)}&\dots&-Y_{f(1,m)}\\ Y_{f(2,1)}&.&\dots&-Y_{f(2,m)}\\ .&.&\dots&.\\ -Y_{f(m,1)}&.&\dots&Y_{f(m,p)}+Y_{f(m,1)}+\dots+Y_{f(m,m)}\\ \end{array}\right]$$

$$\frac{y(k)}{u(k)}=\frac{b_1{z}^{-1}+b_2{z}^{-2}+\cdots +b_n{z}^{-n}}{1+a_1{z}^{-1}+a_2{z}^{-2}+\cdots +a_n{z}^{-n}}$$

$$\displaystyle\left[\begin{array}[]{c}y(k)\\ y(k-1)\\ .\\ .\\ .\\ y(k-N+1)\end{array}\right]_{N\times 1}=\left[X_{N\times 2n}\right]\left[\begin{array}[]{c}a_{1}\\ .\\ .\\ a_{n}\\ b_{1}\\ .\\ .\\ b_{n}\end{array}\right]_{2n\times 1}$$

$${\Phi }_{model(N\times 1)}=X_{N\times 2n}{\Theta }_{2n\times 1}$$

$$ϵ={\Phi }_{measured}-{\Phi }_{model}$$

$$J={ϵ}^tϵ$$

$$\Theta ={\left[X^{t}X\right]}^{-1}{X}^t{\Phi }_{measured}$$

$$\Theta =S^{-1}{X}^t\Phi$$

$$\displaystyle\Theta(k)=S^{-1}\left[X(k)X^{t}(k-1)\right]\left[\begin{array}[]{c}\Phi(k)\\ \Phi(k-1)\end{array}\right]$$

$$\Theta (k)=S^{-1}\left[X(k)\Phi (k)+X^t(k-1)\Phi (k-1)\right]$$

$$\Theta (k)=S^{-1}\left[X(k)\Phi (k)+X^t(k-1)X(k-1)\Theta (k-1)\right]$$

$$\Theta (k)=S^{-1}\left[X(k)\Phi (k)+S(k-1)\Theta (k-1)\right]$$

$$S(k)=S(k-1)+X(k)X^t(k)$$

$$\Theta (k)=S^{-1}\left[X(k)\Phi (k)+\{S(k)-X(k)X^t(k)\}\Theta (k-1)\right]$$

$$\Theta (k)=\Theta (k-1)+[S(k-1)+X(k)X^t(k){]}^{-1}X(k)$$

$$[\Phi (k)-X^t(k)\Theta (k-1)]$$

$$P(k)=P(k-1)\left[I-\frac{X(k)X^t(k)}{1+X^t(k)P(k-1)X(k)}\right]P(k-1)$$

$$K(k)=\frac{X(k)}{1+X^t(k)P(k-1)X(k)}$$

$$P(k)=\left[I-K(k)X^t(k)\right]P(k-1)$$

$$\Theta (k)=\Theta (k-1)+K(k)\left[\Phi (k)-X^t(k)\Theta (k-1)\right]$$

$$K(k)=\frac{P(k-1)X(k)}{\gamma +X^t(k)P(k-1)X(k)}$$

$$P(k)=\frac{\left[I-K(k)X^t(k)\right]P(k-1)}{\gamma }$$

$${\left[Y_d\right]}_{m\times m\times k}={\left[Y_f(z^{-1})\right]}_{m\times m}$$

Full text

Nomenclature

EMT Model

The study and external areas are modeled as EMT type (Original model).

EMT+TSA

The external area is modeled as TSA type equivalent with only network aggregation.

EMT+TSA(AGG)

The external area is modeled as TSA type with both network and generator aggregation.

EMT+FDNE

The external area is modeled as FDNE type generated using the proposed algorithm.

EMT+FDNE (VF)

The external area is modeled as FDNE using the Vector Fitting (VF) from literature.

EMT+FDNE+TSA

The external area is modeled as a combination of FDNE and TSA type with only network aggregation.

EMT+FDNE+TSA (AGG)

The external area is modeled as a combination of FDNE and TSA type with both network and generator aggregation.

I Introduction

Real-time EMT simulation requires detailed modeling of transmission systems to understand the effect of transients and harmonics arising due to varying operating conditions and disturbances in power grid. Effect of power electronic components associated with renewable energy sources on the power grid and performance of different controllers can be analyzed using EMT simulations [1]. Typically, integration time step of EMT simulation is in microseconds ($\mu s$). This makes modeling of a large transmission system for EMT studies impractical as detail modeling increases complexity and computational burden. One solution is to model the transmission system as TSA type with larger integration time as TSA simulations can run faster than EMT. However, in TSA type, due to large integration time, high-frequency transients following a disturbance in the system are not preserved making this approach not very accurate [2].

Another approach is to model large transmission network as frequency dependent reduced order systems that can represent the power grid under any operating condition. One way to reduce large power grid is to model part of the grid which is of interest (“study area”) in detail and the remainder of the system (“external area”) by an efficient equivalent such as FDNE [3]. In this process, initially, the network is divided into study and external area based on the coherency grouping of generators such that all coherent generators are present in the external area [4]. The boundary between the study and external area is divided considering the fact that interconnecting points should have the least minimum number of ports. Generally, in TSA type equivalent, the network admittance is evaluated only at the fundamental frequency, hence this representation ignores high-frequency oscillations. The high-frequency behavior of the external area can be preserved by using FDNE [5], however, FDNE ignores electromechanical oscillations. In order to cover both electromechanical and high-frequency behavior, the external area should be modeled as a combination of TSA equivalent and FDNE.

FDNEs are generally formulated as frequency-dependent black-box terminal equivalents based on rational functions. In [6], FDNE is formulated using modal vector fitting (VF) method and in [7] passivity of the rational function is enforced using mode revealing algorithm. The authors in [8], used VF method to generate FDNE combined with TSA component for real-time simulation and, in [9] this method was improved by including generator coherency. Further, the applicability of these models for large power system real-time simulation is studied in [10]. Formulation of FDNE using time domain VF and simulated time domain response evaluations are performed in [3], and in [11], authors compared different rational approximation methods for simulated time domain responses. A rational approximation approach for formulating FDNE based on matrix pencil method is proposed in [12] where as [13] proposed a method for electromechanical and electromagnetic transient analysis using FDNE.

However, state-of-the-art algorithms formulate FDNE with continuous domain transfer function which makes it very complex to implement these algorithms in the real-time simulator (EMT based). Also, in many of these methods FDNE formulation is dependent on the availability of admittance data over a wide range of frequency which is not always possible. Further, off-line calculation of admittance over a wide frequency range, and storing and retrieving of the data for curve fitting is tedious and time consuming. For example, if the frequency range of interest is from 0 to 5000 Hz with a step size of 0.1 Hz, one should construct 50001 admittance matrices at each frequency sample which is extremely complex to compute. This complexity increases with sample size, number of ports, and size of network under consideration.

This paper introduces a novel method for formulating FDNEs based on online RLS identification. In this method, the external area is energized with constant voltage source at the boundary buses after all the voltage sources and current sources in the network are short-circuited and open circuited respectively. Subsequently, by tracking input voltage and output current, FDNE is formulated. The proposed method simplifies FDNE formulation as it is independent of the availability of network parameters over a wide frequency range and its formulation is in discrete domain directly. Also the computational burden is reduced by using Kron’s node elimination method in the external area network. The advantages of the proposed architecture are:

• •

FDNE formulation is independent of the availability of network parameters over a wide frequency range.

• •

The architecture formulates FDNE in discrete domain, which reduces complexity in interfacing FDNE with the real-time simulator.

• •

The architecture formulates coherency based network equivalents of complex networks with less computational burden and desired accuracy.

• •

The methodology enforces stability and passivity conditions to ensure stable EMT simulations.

• •

The methodology can be directly implemented for real-time control [14].

The rest of the paper is organized as follows: In section II overall architecture is discussed. In section III implementation test on the interconnected power grid is discussed and Section IV concludes the paper.

II Proposed Methodology

In large power grids, it has been frequently observed that, after disturbances, the generators swing together in groups, meaning units near a disturbance respond faster and together, whereas distant machines show relatively damped oscillatory behavior. This physical property is known as coherency and group of machines with similar responses are termed coherent generators. In our work, the generators are coherently grouped based on a localness index. Further, power system model order reduction is performed by dividing the original system into a study area and an external area. The proposed method further divides the external area into two parts. The first one is a low-frequency equivalent (TSA) and the second one is high-frequency equivalent (FDNE). The reduced order modeling of power system involves the following steps:

II-A Aggregation of External Area for TSA type modeling and Real-time Integration

For retaining the electromechanical behavior of the system under consideration, aggregated TSA model is used. Fig. 1 shows the flowchart for TSA type modeling.

II-A1 Coherency Based Classification of Study and External Area

For coherency grouping of the power system under study, initially the small signal stability study of the power grid model is performed and the generator electromechanical modes of oscillation are evaluated. Further, based on participation factors of all the generators, a localness index is calculated as follows:

[TABLE]

where $n$ is a number of synchronous generators connected in the system and, $P_{ki}$ is the normalized participation factor of the $k^{th}$ machine in the $i^{th}$ mode. For example, Table. I shows the coherency grouping of generators for IEEE 39 bus system based on the localness index. More details of the localness index are discussed in [15].

II-A2 Network Aggregation

For aggregating the external area network, the admittance matrix $(Y_{n\times n})$ of the external area is formulated using the bus and line data (at 60Hz). For example, if there are $n$ number of buses in the external area and we want to retain $m$ buses (i.e. $m=i+j$ , where $j$ be the number of boundary buses and $i$ is the number of generator buses) and eliminate the remaining $n-m$ buses in the external area, then using Kron node elimination method [16], reduced admittance matrix $(Y_{red})$ can be obtained as:

[TABLE]

[TABLE]

where subscript $b$ and $g$ represents the boundary and generator buses respectively.

II-A3 Generator and Associated Controller Aggregation

After network aggregation and generators are left intact, the reduced admittance matrix is of the size $m\times m$ ($m=i+j$, i.e. there are $j$ boundary buses and $i$ generator buses). To further reduce computational burden and to reduce complexity in modeling, generators and associated controllers can be aggregated. With generator aggregation, the reduced admittance matrix is of the size $(m-i+1)\times(m-i+1)$ ($m=i+j$, i.e. $i$ coherent generators can be aggregated into one generator). In the case where generators are not aggregated the accuracy of coherency grouping has no significant effect on the TSA, but in the case of generator aggregation, all the generators which are to be aggregated must be in the same coherent group. For scenarios where multiple operating condition changes, the grouping can be performed using online coherency grouping [17]. Method of generator and controller aggregation is discussed in [18]. Thus additional details are not explained in the paper.

II-A4 Interfacing TSA type modeling with real-time simulator

In this step, voltages at the boundary buses are the input to the TSA block, where output currents from the TSA block are injected back to the boundary buses. Here, the generators are modeled in detail to observe the electromechanical behavior. Conversion of boundary bus voltage from time domain to phasor domain is then performed using discrete sequence analyzer [19]. The result gives magnitude $|V_{b}|$ and, phase angle $\angle V_{b}$. Also, phase angle of $V_{b}$ with reference to $I_{b}$ can be determined as $\theta_{b}=\angle\delta_{b}-\angle V_{b}$, where $\angle\delta_{b}$ is the angle of $I_{b}$. Then, using (8) the generator bus voltage can be calculated as

[TABLE]

where $I_{g}$ is the generator current injection, $\angle\delta_{g}$ is the angle of generator current, $V_{b}$ is the boundary bus voltage, $V_{g}$ is the generator bus voltage and, $\angle\theta_{g}$ is the angle of generator bus voltage. Using $I_{g}$ (for the first iteration the initial value is obtained from the power flow solution) and $V_{b}$, the generator voltage $V_{g}$ is calculated. The generator bus is energized with $V_{g}$ and subsequently $I_{g}$ is obtained from generator phasor model. From calculated $V_{g}$ and boundary bus voltage $V_{b}$, boundary bus current injection $I_{inj}$ is calculated as shown in (9)-(10). A graphical description is shown in Fig. 2.

[TABLE]

[TABLE]

After calculating boundary bus current $I_{inj}\angle\beta_{b}$ in phasor form, it is converted into the time domain and injected into the boundary bus. The overall implementation is shown in Fig. 1.

II-B Study and analysis of TSA type equivalent modeling

For study and preliminary analysis, the TSA type equivalent is implemented on two area power system model as shown in Fig. 3. In this test system, area-1 consists of generators G1, G2, and area-2 consists of generators G3, G4 [20]. For analysis purpose, area-2 is considered as external area with boundary bus as bus 10. Two cases are considered here. In the first case EMT+TSA Based Model and in the second case EMT+TSA Based Model (AGG) are analyzed. For bench-marking, both the test cases are compared with (EMT Based Model).

For validating the proposed approach a 3-ph fault is created on Bus-8 at 1s for a duration of 0.1s. Fig. 4 shows the comparison of the relative speed of Gen. 2 with respect to (w.r.t) gen. 1, and Fig. 5 shows the comparison of active power flow from bus 10 to bus 9. In general, after a disturbance, high-frequency transients occur for a short duration whereas electro-mechanical oscillation will be for longer duration. From Fig. 4, it can be observed that due to large simulation time step compared to EMT type, TSA type equivalent can only preserve electro-mechanical oscillations, but high-frequency oscillations are not preserved as it can be seen from Fig. 5.

For quantitative analysis, relative error between the two cases and full EMT type model is calculated using (11). The results are tabulated in Table. II.

[TABLE]

where $y_{ref}$ represents the output from full EMT model (EMT Based Model) and $y_{act}$ represents the output obtained in each case (EMT+TSA Based Model (AGG) and (EMT+TSA Based Model) respectively)

It can be seen that the aggregated model has a large error. This motivates the use of FDNE representations.

II-C FDNE Formulation and Real-time simulator Integration

For retaining the high-frequency behavior of the system under consideration, FDNE is formulated. Fig. 6 shows the FDNE modeling flowchart. In this method, the external area is energized with a constant voltage source with varying frequency in steps of 0.01 (from few Hz to kHz) after short-circuiting all voltage sources, and open circuiting all current sources. Since, the goal here is to identify the frequency dependent admittance $(Y_{f})$, the boundary bus voltage and current are the required signals. FDNE is then formulated in $z$-domain using RLS [21, 22] by tracking input voltage and output current. The basic principle is as follows. If $V_{F}$ is the voltage input to the boundary bus and $I_{F}$ is the current output from the boundary bus, then $Y_{f}$ can be written as

[TABLE]

where $k$ is the number of samples. For a $m$-port network (that means $m$ boundary buses), $Y_{f}$ can be represented as in (13), where $Y_{f(m,m)}$ and $Y_{f(m,p)}$ in (13) are the self and mutual admittance respectively. FDNE model formulation and validation involves the following steps.

II-C1 Recursive Least Square Estimation

Identification of a dynamic process is performed using the process input $u(k)$ and the process output $y(k)$ at every sample $k$ . Considering the $z$-domain model of an $n^{th}$ order process, this can be represented as

[TABLE]

where $a^{\prime}s$ and $b^{\prime}s$ are the transfer function denominator and numerator coefficients respectively. For $N$ observation window length, (14) can be rewritten as

[TABLE]

Equation (29) can be written in the generic form as

[TABLE]

where $X$ is a matrix of past inputs and outputs, $\Phi$ is a matrix of past and present outputs, and $\Theta$ is the coefficient matrix of the transfer function.

Assume that the model identified is different from measurements, then

[TABLE]

where $\epsilon$ is the error between the performance of the system measurement (subscript measured) and the model (subscript model). For reducing this error, a criteria $J$ can be defined as

[TABLE]

By letting $dJ/d\Theta=0$, we get

[TABLE]

From (33) it can be seen that, in order to calculate the measured variable the inverse of the state matrix should be determined. This can drastically slow down the process and some time may not be achievable. To circumvent this issue, an RLS based computational algorithm that eliminates the matrix inversion is designed. Let $S=X^{t}X$, then (33) can be written as

[TABLE]

where $\Phi=\Phi_{measured}$ Then,

[TABLE]

[TABLE]

Using (30), (38) can be written as

[TABLE]

[TABLE]

[TABLE]

Substituting (41) in (40)

[TABLE]

[TABLE]

Let $P(k)=S^{-1}(k)$. Then by using matrix inversion lemma, $P(k)$ can be represented as

[TABLE]

Let

[TABLE]

Then, $P(k)$ can be re-written as

[TABLE]

Substituting (46) in (43), (43) can be represented as

[TABLE]

With weighted least square, (45) and (46) can be presented as

[TABLE]

[TABLE]

where $\gamma$ is the weighting factor.

Thus, with a given process input $V_{F}(k)$ and process output $I_{F}(k)$, $Y_{f}$ can be computed using RLS estimation. The validity of the proposed algorithm is verified by implementing on different test systems. In the first case, the proposed algorithm is implemented on two area test system with 1-port network. Fig. 7(a) and Fig. 7(b) shows the magnitude and angle of admittance of the external area and Table. III shows the comparison of FDNE formulation for two area system. It can be observed that even though this approach uses lower order transfer function (meaning less computational burden), it gives similar error compared to higher order VF method.

II-C2 Passivity Enforcement

For a stable EMT simulation, the admittance matrix should be passive. Presence of negative resistance (i.e angle of admittance $>$ $\pm 1.57\ rad$ in phasor form) in the frequency range of interest due to inherent approximations in identification violates passivity. For instance, from Fig. 7(b) it can be seen that for the previous test case, passivity is violated. For enforcing passivity, an algorithm is developed as shown in Algorithm 1. The details are as follows.

From (13), let $Y_{f}(z^{-1})$ is the fitted admittance transfer function matrix which of the size $m\times m$ (where $m$ is the number of ports). By substituting $z=e^{i2\pi fT_{s}}$ in (13), the fitted admittance matrix data for $k$ frequency samples can be represented as

[TABLE]

where $f\in 1:f_{max}$ ($f_{max}$ being the maximum frequency under consideration) and, $T_{s}$ is the sampling time. Let $G_{d}$ is the real part of the admittance matrix $(Y_{d})$. Then

[TABLE]

For a function to be passive

[TABLE]

This implies that, if the admittance transfer function matrix is positive definite then it is also passive. Considering this, if the fitted function $Y_{f}$ violates (52), then a new corrected transfer function matrix $Y_{c}$ is formulated. The conductance transfer function $(G_{f})$ and susceptance transfer function ($B_{f}$) are calculated as follows.

[TABLE]

[TABLE]

[TABLE]

where $*$ stands for complex conjugate.

Since passivity is related to real part of admittance matrix $(G_{f})$, correcting $G_{f}$ without affecting imaginary part $(B_{f})$ is sufficient. Then $Y_{c}$ can be represented as follows

[TABLE]

where

[TABLE]

The objective here is to calculate $\Delta G$. A real, symmetric matrix $G_{b}$ is said to be positive definite if $x^{T}G_{b}x$ $>$ [math] $\forall$ $x$ $\neq$ [math]. Thus $x^{T}G_{b}x$ can be written as

[TABLE]

This shows $G_{b}$ is positive definite if and only if $G_{b}+G_{b}^{T}$ is positive definite. This can be achieved by minimizing an objective function through optimization as

[TABLE]

[TABLE]

where $F$ stands for Forbenius norm of a matrix. We propose a convex optimization formulation to find $G_{b}$ using semi definite programming (SDP) [23]. The optimization solution is then used to calculate $\Delta G$ using (61).

[TABLE]

The validity of the proposed algorithm is implemented for enforcing passivity of 1-port network formulated previously. From Fig. 8 it can be seen that the passivity is enforced.

II-C3 Interfacing FDNE with real time simulator

FDNE can be directly implemented since it is computed in $z$-domain. The implementation process is as follows. With boundary bus voltage $(V_{b})$ as input to FDNE, and with $n^{th}$ order estimation, (12) can be written as

[TABLE]

where $I_{F}$ is current output from FDNE. For observing high frequency transients only FDNE part is required. To maintain boundary bus parameters at initial steady state, a constant current source is injected into the boundary bus as calculated from (63) [8, 9]. This can be represented as,

[TABLE]

[TABLE]

where $P_{b}$ and $Q_{b}$ are the active and reactive power flow respectively from the boundary bus, $V_{b}$ and $\theta_{b}$ are the voltage and angle respectively of the boundary bus. Since admittance at a fundamental component of frequency $(Y_{c}(60Hz))$ is included either in (9) for EMT+FDNE+TSA Based Model or in (63) for EMT+FDNE Based Model, the fundamental frequency component must be eliminated from FDNE. This is performed by subtracting $Y(60Hz)V_{b}\angle\theta_{b}$ term in (64) before injecting boundary bus current to remove fundamental frequency component from FDNE.

II-D Study and analysis of FDNE type equivalent modeling

For study and preliminary analysis, the proposed FDNE algorithm is implemented on two area test system. Two cases are considered for study. First case is the proposed method EMT+FDNE Based Model (Proposed) and the second case is an existing architecture (EMT+FDNE Based Model(VF)). For validation, both the test cases are compared with EMT Based Model. For the same type of disturbance as discussed in Section II-B, Fig. 10 shows the comparison of the relative speed of Gen. 2 w.r.t gen. 1, and Fig. 11 shows the comparison of active power flow from bus 10 to bus 9.

Table IV shows the comparisons between proposed FDNE and an offline VF based algorithm. Both algorithms gives similar results, proving that FDNE can be formulated online with less computational effort and lower order of transfer function when compared to offline algorithms. From Fig. 10 and Fig. 11 it can be seen that with this approach, high frequency oscillations are preserved whereas electromechanical oscillations are not preserved, proving the need to have combined FDNE models with TSA equivalents.

II-E Implementing TSA and FDNE on two area power system

Fig. 12 shows the implementation approach for combined TSA and FDNE type equivalents. Table. V shows the comparison of frequency and damping factor for original model and various reduced order model using eigen value realization algorithm for generator-3 speed data. It can be seen that the proposed approach provides very close result compared to full EMT based model. Fig. 13 shows the comparison of the relative speed of Gen. 2 w.r.t gen. 1, and Fig. 14 shows the comparison of active power flow from bus 10 to bus 9. From Fig. 13 and Fig. 14, it can be seen that both high-frequency and electromechanical oscillations are well preserved. Table. VI shows the error comparison between orignial and reduced model using (11).

III Implementation Test on Interconnected Power Grid

To prove scalability and implementation using a multi-port networks, IEEE 39 and 68 bus power system models are considered. Plesae refer [24] and [25] for IEEE 39 and 68 bus system details and one-line diagrams respectively.

III-A IEEE 39 Bus System Implementation Test Results

In this system first, based on the coherency grouping of the generators the test system is divided into study and external area. To assess the performance of proposed EMT+FDNE+TSA based reduced order model, Group-I which consists of generators 4, 5, 6, 7, 9 is considered as external area and the rest of the power system as study area. The study and external area is divided at bus 16 (Port-1), bus 17 (Port-2), and bus 26 (Port-3). Fig. 15(a) and Fig. 15(b) shows the admittance magnitude and angle of port-3. Fig. 16 shows the passivity enforcement for 3-port network.

For analysis, a three-phase fault is initialized for a duration of 0.1s and the simulation results are compared with the original model. Fig. 17 to Fig. 19 shows the validation results of the proposed algorithm. Table VII shows relative error (11) comparison of reduced and original models. It can be seen that the reduced model closely resembles the full EMT model.

III-B IEEE 68 Bus System Implementation Test Results

Further, in order to validate the proposed algorithm on a larger system, IEEE 68 bus system which consists of 16 generators consists of 5-areas is considered. The area with generators 1 to 9 is considered as $external$ area. The $study$ and $external$ areas are divided at bus 54 (Port-2), 27 (Port-3), and 60 (Port-1). Fig. 20(a) and Fig. 20(b) shows the admittance magnitude and angle of port-2. For analysis a 3-ph fault for a duration of 0.2s is simulated at bus-49 starting at at 1s. The responses of reduced order model are compared with the original model. Fig. 21 shows the Bus-60 voltage and Fig. 22 shows the relative speed of generator-10 w.r.t generator-16.

From the above comparisons it can be seen that the dynamics of the proposed reduced order model is similar to that of the original model. The computational time required to generate equivalents is shown in Table. VIII.

IV Conclusion

In this paper a novel real-time frequency based reduced order modeling of the large power system for EMT simulation is proposed. In the proposed architecture, external area is modeled as a combination of FDNE and TSA. FDNE is formulated using an online discrete RLS which can preserve high-frequency behavior, where as TSA preservers the electromechanical (low frequency) behavior of the system under consideration. The approach also enforces passivity conditions. Implementation results in two area, IEEE 39 and 68 bus power system models shows that the proposed reduced order model represents the aggregated power grid accurately and at the same time can capture the oscillations as close as the full detail model for EMT simulations.