Modeling and Quantifying the Impact of Wind Power Penetration on Power System Coherency

TL;DR

This paper analyzes how wind power integration affects power system coherency by altering the spectral properties of the system's Laplacian matrix, with a theoretical framework validated through simulations and measurement analysis.

Contribution

It introduces a novel theoretical framework to quantify the impact of wind penetration on power system coherency and validates it with numerical and measurement-based methods.

Findings

Wind penetration modifies the Laplacian matrix and affects generator coherency.

The model predicts changes in coherent areas due to wind farm placement and level.

Simulation results confirm the theoretical predictions.

Abstract

This paper presents a mathematical analysis of how wind generation impacts the coherency property of power systems. Coherency arises from time-scale separation in the dynamics of synchronous generators, where generator states inside a coherent area synchronize over a fast time-scale due to stronger coupling, while the areas themselves synchronize over a slower time-scale due to weaker coupling. This time-scale separation is reflected in the form of a spectral separation in the weighted Laplacian matrix describing the swing dynamics of the generators. However, when wind farms with doubly-fed induction generators (DFIG) are integrated in the system then this Laplacian matrix changes based on both the level of wind penetration and the location of the wind farms. The modified Laplacian changes the effective slow eigenspace of the generators. Depending on penetration level, this change may result in changing the identities of the coherent areas. We develop a theoretical framework to quantify this modification, and validate our results with numerical simulations of the IEEE 68-bus system with one and multiple wind farms. We compare our model based results on clustering with results using measurement-based principal component analysis to substantiate our derivations.

Tables

Table 1. TABLE I : Inter-area modes and areas of the nominal system

Slow modes	Frequency in Hz
$-0.1321\pm j1.877$	0.299
$-0.1449\pm j2.843$	0.453
$-0.1434\pm j3.763$	0.599
$-0.1772\pm j3.910$	0.622

Table 2. TABLE II : Inter-area modes and areas (bus 66, γ = 650 𝛾 650 \gamma=650 )

Slow modes	Frequency in Hz
$-0.129\pm j1.761$	0.280
$-0.145\pm j2.825$	0.449
$-0.142\pm j3.752$	0.597
$-0.1436\pm j3.844$	0.612

Table 3. TABLE III : Inter-area modes and areas (bus 37 37 37 , γ 𝛾 \gamma = 700 700 700 )

Slow modes	Frequency in Hz
$-0.1299\pm j1.82$	0.289
$-0.1447\pm j2.83$	0.450
$-0.1416\pm j3.756$	0.598
$-0.1698\pm j3.95$	0.629

Table 4. TABLE IV : Inter-area modes and areas with 3 3 3 wind plants

Slow modes	Frequency in Hz
$-0.1363\pm j2.244$	0.3571
$-0.1666\pm j2.65$	0.4219
$-0.1348\pm j3.118$	0.4964
$-0.1482\pm j3.986$	0.6342

Equations

$${\dot{\delta }}_i$$

$$P_{si}=$$

$$Q_{si}=$$

$$0=$$

$$\Delta \dot{\delta }=I\Delta \omega ,M\Delta \dot{\omega }=K_{11}\Delta \delta +K_{12}\Delta V+\Delta P_m,$$

$$\Delta P_s=-K_{11}\Delta \delta -K_{12}\Delta V,\Delta Q_s=-K_{21}\Delta \delta -K_{22}\Delta V.$$

$$\Delta \ddot{\delta }=M^{-1}\underset{{\mathcal{L}}_0}{\underbrace{(K_{11}-K_{12}{A}_3^{-1}{A}_1)}}\Delta \delta ,$$

$$A_3={\left[\begin{array}{cccc}(-K_1-K_{12}{)}^T& K_1^{{}^{\prime }T}& (-K_2-K_{22}{)}^T& K_2^{{}^{\prime }T}\end{array}\right]}^T$$

$$A_1={\left[-K_{11}^{T}0-K_{21}^{T}0\right]}^T.$$

$$q_s^{\alpha }=\frac{{\sum }_{i=1}^{n_{\alpha }}{M}_i^{\alpha }\Delta {\delta }_i^{\alpha }}{M^{\alpha }},q_{fj}^{\alpha }=\Delta {\delta }_j^{\alpha }-\Delta {\delta }_1^{\alpha },$$

$$\left[\begin{array}[]{c}q_{s}\\ q_{f}\end{array}\right]=\begin{bmatrix}\hat{M}^{-1}U^{T}M\\ G\end{bmatrix}\Delta\delta:=\begin{bmatrix}C\\ G\end{bmatrix}\Delta\delta,$$

$$\left[\begin{array}{c}{\ddot{q}}_s\\ {\ddot{q}}_f\end{array}\right]=\left[\begin{array}{cc}{T}_{11}& T_{12}\\ T_{21}& T_{22}\end{array}\right]\left[\begin{array}{c}{q}_s\\ q_f\end{array}\right]:=T\left[\begin{array}{c}{q}_s\\ q_f\end{array}\right],$$

$$T_{11}=ϵC{M}^{-1}{\mathcal{L}}_0^{E}U,T_{12}=ϵC{M}^{-1}{\mathcal{L}}_0^E{G}^{†},T_{21}=ϵG{M}^{-1}{\mathcal{L}}_0^{E}U,$$

$$T_{22}=G{M}^{-1}({\mathcal{L}}_0^I)G^{†}+ϵG{M}^{-1}{\mathcal{L}}_0^E{G}^{†},$$

$${\mathcal{L}}_0={\mathcal{L}}_0^I+ϵ{\mathcal{L}}_0^E,$$

$$J_r{\dot{\omega }}_r(t)$$

$$J_g{\dot{\omega }}_g(t)$$

$${\dot{\theta }}_T(t)$$

$$v_{ds}$$

$$v_{dr}$$

$$v_{qr}$$

$${\psi }_{qs}=L_{ls}{i}_{qs}+L_m(i_{qs}+i_{qr}),{\psi }_{ds}=L_{ls}{i}_{ds}+L_m(i_{ds}+i_{dr}),$$

$${\psi }_{qr}=L_{lr}{i}_{qr}+L_m(i_{qs}+i_{qr}),{\psi }_{dr}=L_{lr}{i}_{dr}+L_m(i_{ds}+i_{dr}),$$

$$T_g(t)=\frac{3}{2}\frac{p}{2}{L}_m\left[i_{qs}(t)i_{dr}(t)-i_{ds}(t)i_{qr}(t)\right].$$

$$P_w=\gamma (v_{qs}{i}_{qs}+v_{ds}{i}_{ds}),Q_w=\gamma (-v_{ds}{i}_{qs}+v_{qs}{i}_{ds}).$$

$$\Delta \dot{\delta }=I\Delta \omega ,M\Delta \dot{\omega }={\bar{K}}_{11}\Delta \delta +{\bar{K}}_{12}\Delta V+\Delta P_m,$$

$$\Delta P_s=-{\bar{K}}_{11}\Delta \delta -{\bar{K}}_{12}\Delta V,\Delta Q_s=-{\bar{K}}_{21}\Delta \delta -{\bar{K}}_{22}\Delta V.$$

$$\Delta \dot{z}=A\Delta z+B\Delta V.$$

$$\displaystyle\hskip 36.98866pt\left[{\begin{array}[]{*{20}{c}}{{\Delta P_{w}}}\\ {{\Delta Q_{w}}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{{C_{1}}\Delta z}\\ {{C_{2}}\Delta z}\end{array}}\right]+\left[{\begin{array}[]{*{20}{c}}{{D_{1}}\Delta V}\\ {{D_{2}}\Delta V}\end{array}}\right],$$

$$\displaystyle C_{1}=\frac{\partial P_{w}}{\partial z}\Bigr{|}_{\hat{p}_{0}},C_{2}=\frac{\partial Q_{w}}{\partial z}\Bigr{|}_{\hat{p}_{0}},D_{1}=\frac{\partial P_{w}}{\partial V}\Bigr{|}_{\hat{p}_{0}},D_{2}=\frac{\partial Q_{w}}{\partial V}\Bigr{|}_{\hat{p}_{0}}.$$

$$0={\bar{A}}_1\Delta \delta +{\bar{A}}_2\Delta z+{\bar{A}}_3\Delta V,$$

$${\bar{A}}_1$$

$${\bar{A}}_3$$

$$(D_2-{\bar{K}}_4){\bar{K}}_6^T{]}^T.$$

$$\displaystyle\left[\begin{array}[]{c}\Delta\dot{\delta}\\ \Delta\dot{\omega}\\ \Delta\dot{z}\end{array}\right]=\begin{bmatrix}0&I&0\\ \mathcal{R}_{1}&0&\mathcal{R}_{2}\\ \mathcal{R}_{3}&0&\mathcal{R}_{4}\end{bmatrix}\left[\begin{array}[]{c}\Delta\delta\\ \Delta\omega\\ \Delta z\end{array}\right],$$

$${\bar{K}}_{11}=K_{11}+{\Delta }_{k11}(\gamma ),{\bar{K}}_{12}=K_{12}+{\Delta }_{k12}(\gamma ).$$

$$\displaystyle D_{1}=\bordermatrix{&&\underset{\downarrow}{(n+1)}&&&\underset{\downarrow}{(n+m+1)}&&\cr&\underbrace{0}_{\frac{\text{sync-gens}}{1\times n}}&\underbrace{\frac{\partial P_{w}}{\partial V_{(n+1)_{Re}}}\Bigr{|}_{\hat{p}_{0}}}_{\frac{\text{wind gen}}{(1\times 1)}}&\underbrace{0}_{\frac{\text{non-gens}}{1\times(m-n-1)}}&0&\frac{\partial P_{w}}{\partial V_{(n+1)_{Im}}}\Bigr{|}_{\hat{p}_{0}}&0},$$

Full text

I Introduction

Over the past two decades, significant amount of research has been done in studying various impacts of wind power penetration on power system dynamics, stability, and control. Results reported in [1, 2, 3], for example, demonstrate the impact of wind integration on transient stability and small-signal stability. Results in [4] show the low-inertial effects of wind integration resulting in degradation of frequency control. Results in [5] overview the effects of penetration level of wind farms on the damping characteristics of inter-area oscillations. Multiple studies have also been done on developing new control schemes for improving the steady-state and dynamic operation [6], and damping of inter-area oscillations in wind-integrated power systems [7, 8, 9]. In [10], an analytical framework has been developed to evaluate the effect of wind penetration on time-scale separation properties of power systems that lead to inter-area oscillations. The majority of these works describe how wind power impacts the eigenvalues of the small-signal model of a power system. However, in order to understand its holistic effect, it is equally important to evaluate how the eigenvectors of these models change as more DFIGs penetrate the grid.

In this paper we address this problem by studying the impact of wind penetration on a specific eigen-property of power systems that requires the use of eigenvectors - namely, coherency [11],[12]. Coherency is a fundamental property of power systems that arise from the separation of time-scales in the dynamics of synchronous machines. Machines that are strongly coupled tend to swing together, and synchronize over a fast time-scale, thereby forming coherent groups or clusters, while the groups themselves swing against each other and synchronize over a slow time-scale due to their weaker coupling. A large literature exists on coherency theory, for example see [13, 14, 15, 16], which to date is still highly useful for transmission planning and operations including dynamic equivalencing [17], controlled islanding [18], and oscillation damping control [19]. Recent papers such as [20, 21] have used data-driven techniques to evaluate how the traditional notions of coherency are affected by non-synchronous power generation from wind, but no theoretical reasonings have been presented. The results derived in this paper compensate for this gap by inferring that the intrusion of DFIGs can be viewed as an addition of heterogeneity to the homogeneous dynamics of synchronous machines. This heterogeneity changes the effective dynamic coupling between the synchronous machines, and thereby perturbs the eigenvalues and eigenvectors of the swing dynamics such that the coherent groups change. The change depends on the amount of wind power injected, and the bus locations of the wind farms.

The main contributions of the work can be summarized as follows. We compare electro-mechanical models of multi-machine power systems both with and without wind injection, and quantify the perturbation caused by this injection in the weighted Laplacian matrix associated with the swing dynamics of the synchronous generators. The quantification is done in terms of both penetration level and location of wind plants. Using singular perturbation theory, an equivalent Laplacian matrix is derived to capture the modified interaction between the synchronous generators in presence of wind. Thereafter, a coherency grouping algorithm is stated in terms of the eigenvectors of this modified Laplacian matrix. A motivational example for the study is drawn from the model of the New York state power grid with large-scale wind penetration. The theoretical analyses and algorithms are all verified using simulations of the IEEE benchmark 16-machine, 68-bus power system model with multiple wind plants.

II Recapitulation of Coherency

We first recall the fundamental theory of coherency in synchronous machines, as detailed in [11]. Consider a power system with $m$ buses and $n$ synchronous generators. Considering classical model of synchronous generators [22], the dynamics of the $i^{th}$ generator can be written as,

[TABLE]

where $\delta_{i}$, $\omega_{i}$, $M_{i}$, $x^{\prime}_{di}$, $E_{i}$, $P_{mi}$ are respectively the phase angle, machine speed deviation from nominal speed ($120\pi$ rad/s), inertia, direct-axis transient reactance, internal machine voltage, and the mechanical power input to generator $i\;(i=1,\dots,n)$. $V_{i_{Re}}$ and ${V_{i_{Im}}}$ are the real and imaginary parts of the bus voltage phasor. The active and reactive power outputs of the $i^{th}$ generator can be written as,

[TABLE]

The active and reactive power flow balance at any bus $j,j=1,\dots,m$ can be written as,

[TABLE]

where $G_{j}$ and $B_{j}$ are the conductance and the susceptance of the shunt load at bus $j$ with line charging. Assuming the transmission lines to be lossless, $B_{jk}$ denotes the susceptance of the tie-line connecting bus $j$ and bus $k$. Linearizing (1) and (2) about a stable operating point $p_{0}=\{\delta_{0},0,V_{Re_{0}},V_{Im_{0}}\}$ governed by the power flow solution we get,

[TABLE]

Here, $M$ is the diagonal matrix of the machine inertias, $K_{ij},i,j=\{1,2\}$ are the Jacobian matrices of appropriate dimensions following from (1)-(3), and $I$ is the identity matrix. $\Delta\delta,\Delta\omega,\Delta P_{m},\Delta P_{s},\Delta Q_{s}$ and $\Delta V$ are the vectors constructed by stacking the state, the input and the algebraic variables. Specifically, $\Delta V=[\Delta V_{1_{Re}}\;\dots\;\Delta V_{m_{Re}}\;\Delta V_{1_{Im}}\;\dots\;\Delta V_{m_{Im}}]^{T}$. Following Kron-reduction, and the unforced small-signal model (4) can be written as

[TABLE]

[TABLE]

The Jacobian matrices $K_{1},K_{2}$ follow from the linearization of (3) for the generator buses, while $K^{\prime}_{1},K^{\prime}_{2}$ are those for non-generator buses.

We impose a two-time scale behavior on (5) by assuming the network to be divided into $r$ distinct and non-overlapping coherent areas [11]. The two time-scale model can be derived as follows. Let there be $n^{\alpha}$ generators in area $\alpha,\alpha=1,\dots,r$. Let $\Delta\delta_{i}^{\alpha}$ and $M_{i}^{\alpha}$ be the small-signal phase angle and inertia of the $i^{th}$ machine in area $\alpha$. Define two variables $q_{s}^{\alpha}\in\mathbb{R},q_{f}^{\alpha}=\mbox{col}(q_{fj}^{\alpha})\in\mathbb{R}^{n_{\alpha}-1}$ for $\alpha=1,\dots,r,j=2,\dots,n_{\alpha}$ as,

[TABLE]

where $M^{\alpha}=\sum_{i=1}^{n_{\alpha}}M_{i}^{\alpha}$. Stacking $q_{s}^{\alpha}$ and $q_{f}^{\alpha}$ into vectors for $\alpha=1,\dots,r$, one can write,

[TABLE]

where, $\hat{M}=diag(M^{1},M^{2},..,M^{r})$, $U=blockdiag(U_{1},U_{2},..,U_{r})$ with $U_{\alpha}\in\mathbb{R}^{n_{\alpha}}$ being the vector of all ones, $G=blockdiag(G_{1},G_{2},..,G_{r})$ where definition of $G_{\alpha}$ can be found in [11]. The transformation (7) is invertible with the inverse given by $[U\ \ G^{T}(GG^{T})^{-1}]$. Using (7), one can rewrite (5) in the time-scale separated form,

[TABLE]

[TABLE]

where $G^{\dagger}=G^{T}(GG^{T})^{-1}$, the matrices $\mathcal{L}_{0}^{I},\mathcal{L}_{0}^{E}$ follow from partitioning $\mathcal{L}_{0}$ as

[TABLE]

and $0<\epsilon\ll 1$ is a singular perturbation parameter arising from the worst-case ratio of the tie-line reactances internal and external to the coherent areas. For precise definition of $\epsilon$, please see [11]. The model (8) is in the singularly perturbed form where $q_{s},q_{f}$ are the slow and fast variables. For small values of $\epsilon$, (5) will exhibit a two time-scale behavior, reflected through one DC mode, $n-r$ fast oscillation modes, and $r-1$ slow oscillation modes. The $r$ coherent groups can be identified from the eigen-analysis of $M^{-1}\mathcal{L}_{0}$. Algorithm 1 as in [11] recalls the steps for this identification.

III Motivating Example

We next provide a motivating example from the New York State (NYS) power grid to show how wind penetration can change coherency. The utility-scale model of the NYS grid is simulated using the PSS/E. The model consists of over 70,000 buses, and thousands of dynamic elements. The grid is divided into eleven zones following the NYISO zonal separation based on the similarity of frequency responses of generators in any zone following contingencies. In each zone, a representative bus is chosen corresponding to the largest generating unit. The system is excited with different NYISO-specified contingencies. Time responses of the frequencies at these buses are recorded in a matrix $\mathbb{M}\in\mathbb{R}^{\rho\times s}$. For this example $\rho=11$ is the number of representative buses, and $s=3571$ is the number of data samples. The simulations were run for $15$ seconds following the contingencies, and a sampling time of $0.0042$ second was used.

Next, Principal Component Analysis (PCA) [23] is applied on $\mathbb{M}$ with the objective of expressing this matrix as $\mathbb{M}=\mathcal{K}\mathcal{V}^{T}$ where columns of $\mathcal{V}$ are orthonormal basis vectors. To achieve this the following steps are applied.

Step 1: The singular value decomposition of $\mathbb{M}$ is computed as $\mathbb{M}=U\Sigma V^{T}$ where $U\in\mathbb{R}^{\rho\times\rho}$ and $V\in\mathbb{R}^{s\times s}$ are respectively the matrices of left and right singular vectors, and $\Sigma\in\mathbb{R}^{\rho\times s}$ is the diagonal matrix of singular values.

Step 2: Since the columns of $V$ are the normalized right eigenvectors of $\mathbb{M}^{T}\mathbb{M}$, one can write $\mathcal{K}=U\Sigma,\mathcal{V}^{T}=V^{T}$.

Step 3: The columns of $\mathcal{K}\in\mathbb{R}^{\rho\times s}$ represent the weighting for each principal component. Since $s\gg\rho$, $\mathbb{M}$ can be expected to be a low rank matrix. The weightings for $c$ dominant principal axes of $\mathbb{M}$ (for any chosen $c$) can be found by identifying the $c$ columns of $\mathcal{K}$ that have highest variance. These $c$ columns of $\mathcal{K}$ are finally plotted in a $c$-dimensional plot.

For our example we consider two scenarios A and B, where the wind penetration in the NYS model is respectively $10\%$ and $8\%$ of the nameplate wind capacity. Most of the wind generation is located in zone $1$ (western NY) and zones $2$-$5$ (central-northern NY). PCA results for these scenarios are shown in Figs. 1-2 for the same contingency, considering $c=2$ principal axes. From the two figures it can be seen that subsystem $1$ in scenario A moves towards the group formed by subsystems $2$ through $5$ in scenario B. On the other hand, subsystem $10$ departs from its own coherent group in scenario A, and forms a separate cluster in the scenario B.

This example shows that depending on operating conditions wind penetration can result in changes in coherent clusters. In the next sections we derive conditions that quantify this movement depending on the amount and the location of wind penetration for the power system model (5).

IV Quantification of perturbation

IV-A Wind farm (WF) model

Consider the $m$-bus power system model introduced in Section II. Without loss of generality we assume that a wind farm is connected at the $(n+1)^{th}$ bus. Following [10], the farm is assumed to consist of parallel combinations of $\gamma$ individual wind turbine and DFIG units. The turbine model is considered as

[TABLE]

where ${T_{a}(t)}=\frac{{\rho{A_{s}}\nu^{3}(t){C_{p}}}}{{2{\omega_{r}}(t)}}$ is the aerodynamic torque input. The physical meanings of all variables can be found in [10]. The DFIG dynamics can be expressed in a power-invariant synchronously rotating d-q reference frame as

[TABLE]

where expressions for flux linkages and the generated electrical torque are given by

[TABLE]

Here, $\omega_{ge}:=\frac{p_{e}}{2}\,\omega_{g}$ is the electrical speed of the rotor of the DFIG, $p_{e}$ is the number of electrical poles. $L_{ls}$, $L_{lr}$, $L_{m}$ are the stator and rotor leakage inductances and the magnetizing inductance, respectively. Standard meanings of the voltage and flux variables can be found in [10] and are skipped here for brevity. The active and reactive power output of the wind farm can be written as

[TABLE]

Aligning $v_{qs}$ with the wind bus voltage phasor, we have ${v_{qs}}=|V_{n+1}|$ and ${v_{ds}=0}$. The DFIG is equipped with active and reactive power control loops with PI controllers whose setpoints are computed using Maximum Power Point Tracking (MPPT) and power flow calculations, respectively.

IV-B Linearized wind-integrated model

Since the stator of the DFIG is directly connected to the wind bus through a step-up transformer, the swing states of the synchronous generators will now be dynamically coupled with the wind farm states $z$ which are considered to be the average of the individual unit states [10]. Let the operating point of the power system after wind integration be denoted as $\hat{p}_{0}=\{\hat{\delta}_{0},0,z_{0},\hat{V}_{Re0},\hat{V}_{Im0}\}$. Equations (1) and (2) are linearized about $\hat{p}_{0}$ as

[TABLE]

The Jacobians $\bar{K}_{ij}$ are different than $K_{ij}$ in (4) due to the shift in operating point from the wind injection. The linearized wind farm model from (12)-(14) around $\hat{p}_{0}$ is written as

[TABLE]

Expressions for $A$ and $B$ are skipped for brevity. Note that $B$ is a zero padded matrix where zeros correspond to the non-wind buses. The linearised power output equations are written as

[TABLE]

The zero-padded matrices $D_{1}$ and $D_{2}$ depend on wind penetration level $\gamma$ and on the location of the wind farm. More detailed structures of $D_{1}$ and $D_{2}$ will be shown shortly. The nonlinear power flow equations (3) are linearized at $\hat{p}_{0}$ with Jacobian matrices $\bar{K}_{1},\bar{K}_{2}$ for active and reactive power flows for synchronous generator buses, $\bar{K}_{3},\bar{K}_{4}$ for wind generator bus, and $\bar{K}_{5},\bar{K}_{6}$ for the non-generator buses as

[TABLE]

[TABLE]

Using (16a), (17) and (25), considering $\Delta P_{m}=0$ we write the final set of state-space equations for the unforced wind-integrated power system model:

[TABLE]

where, $\mathcal{R}_{1}=M^{-1}(\bar{K}_{11}-\bar{K}_{12}\bar{A}_{3}^{-1}\bar{A}_{1})$, $\mathcal{R}_{2}=M^{-1}(-\bar{K}_{12}\bar{A}_{3}^{-1}\bar{A}_{2})$,$\mathcal{R}_{3}=-B\bar{A}_{3}^{-1}\bar{A}_{1}$, $\mathcal{R}_{4}=A-B\bar{A}_{3}^{-1}\bar{A}_{2}$.

We next quantify the perturbation of the wind-integrated model (32) from the wind-less model (5). Note that $\mathcal{R}_{1}=M^{-1}\mathcal{L}$ where $\mathcal{L}=(\bar{K}_{11}-\bar{K}_{12}\bar{A}_{3}^{-1}\bar{A}_{1})$ mimics the role of $\mathcal{L}_{0}$ in (5). However, unlike $\mathcal{L}_{0}$, which is a weighted Laplacian matrix (by definition, a symmetric matrix $\mathcal{L}_{0}=\mathcal{L}_{0}^{T}$ is Laplacian if the diagonal entry of each row equals to the negative sum of the other entries of that row), $\mathcal{L}$ is not a Laplacian matrix. To quantify the difference between $\mathcal{L}_{0}$ and $\mathcal{L}$, we must compare each of their constituent matrices separately. This is shown as follows.

IV-C Quantification of perturbation in $\mathcal{L}$

• •

Perturbation in $K_{11}$ and $K_{12}$:

The nominal Jacobian matrices $K_{11},K_{12}$ are perturbed because of the shift in operating point from $p_{0}$ to $\hat{p}_{0}$. We write this perturbation as

[TABLE]

The unstructured perturbations $\Delta_{k11}(\gamma),\Delta_{k12}(\gamma)$ are implicit functions of $\gamma$ and the location of wind farm, and can be determined numerically for a particular wind integration scenario.

• •

Perturbation in $A_{3}$:

The matrices $K^{\prime}_{1}$ and $K^{\prime}_{2}$ in (5) are rewritten as $K^{\prime}_{1}=[-K_{3}\;K_{5}],K^{\prime}_{2}=[-K_{4}\;K_{6}]$ where $K_{3},K_{4}$ are the Jacobians after linearizing (3) for the nominal model with respect to the $(n+1)^{th}$ bus voltage, and $K_{5},K_{6}$ are that with respect to the voltages of non-generator buses. To compare $A_{3}$ and $\bar{A}_{3}$ we need to consider the structures of $D_{1}$ and $D_{2}$. The entries of $D_{1},D_{2}$ can be partitioned in terms of synchronous generator, wind farm and non-generator buses as follows:

[TABLE]

Simple calculations show that the partial derivatives in the above matrices can be written as \frac{\partial P_{w}}{\partial V_{(n+1)_{Re}}}\Bigr{|}_{\hat{p}_{0}}:=\gamma\zeta_{1},\;\frac{\partial P_{w}}{\partial V_{(n+1)_{Im}}}\Bigr{|}_{\hat{p}_{0}}:=\gamma\zeta_{2},\frac{\partial Q_{w}}{\partial V_{(n+1)_{Re}}}\Bigr{|}_{\hat{p}_{0}}:=\gamma\zeta_{3},\;\frac{\partial Q_{w}}{\partial V_{(n+1)_{Im}}}\Bigr{|}_{\hat{p}_{0}}:=\gamma\zeta_{4}, where $\zeta_{i}$’s are linearization constants depending on the steady-state stator voltage and stator currents of the DFIG. Using the structures of $D_{1}$ and $D_{2}$, we can write:

[TABLE]

where the perturbation term $A^{\prime}_{3}$ has the following sparse structure,

[TABLE]

and $\Delta_{A3}(\gamma)$ captures the change in the nominal Jacobians $K_{12},K_{22}$ and $K_{i},i=1,\dots,6$ due to the operating point shift.

• •

Perturbation in $A_{1}$:

To compare $A_{1}$ in (5) with $\bar{A}_{1}$ in (25), we use (33). We write the perturbation as

[TABLE]

where $\Delta_{A1}(\gamma)=[\Delta_{k11}^{T}\;0\;0\;\Delta_{k12}^{T}\;0\;0]^{T}$ follows from (21).

• •

Perturbation in $\mathcal{L}_{0}$:

We recall from (20) and (5) that $\mathcal{L}=\bar{K}_{11}-\bar{K}_{12}\bar{A}_{3}^{-1}\bar{A}_{1}$ and $\mathcal{L}_{0}=K_{11}-K_{12}A_{3}^{-1}A_{1}$. Using the matrix inversion lemma we get,

[TABLE]

where $X=-(I+A_{3}^{-1}(\gamma A^{\prime}_{3}+\Delta_{A3}(\gamma)))^{-1}A_{3}^{-1}(\gamma A^{\prime}_{3}+\Delta_{A3}(\gamma))A_{3}^{-1}$. Using (33)-(35) we get,

[TABLE]

Here $\kappa_{\mathcal{L}}(\gamma)=\Delta_{k11}-K_{12}\Delta_{A1}-\Delta_{k12}(A_{3}^{-1}+X)(A_{1}+\Delta_{A1}))$, contains all the terms due to change in operating point, while the perturbation term $-K_{12}XA_{1}$ contains explicit information about penetration level and location of the wind farm.

IV-D Extension to multiple wind farms

Equations (21)-(25) provide the expressions for matrix perturbations when there is only one wind farm in the system at the $(n+1)^{th}$ bus. The approach can be easily extended to when the system has $p$ wind farms, say at $(n+1,n+2,\dots,n+p)^{th}$ buses. Let their penetration levels be $\gamma_{1},\gamma_{2},\dots,\gamma_{p}$, respectively. We consider the matrix $D_{1}=[D_{1}^{Re},D_{1}^{Im}]$, where $D_{1}^{Re}\in\mathbb{R}^{p\times m}$ is given by

[TABLE]

The definition of $D_{1}^{Im}$ is the same as (26) but with $V_{i_{Re}}$ replaced by $V_{i_{Im}},i=n+1,\dots,n+p$. Similarly, $D_{2}=[D_{2}^{Re},D_{2}^{Im}]$ can be defined in the same way by replacing $P_{w}$ by $Q_{w}$ in (26). As before, the partial derivatives are written as $\frac{\partial P_{w}}{\partial V_{i_{Re}}}=\gamma_{i}\zeta_{1i},\frac{\partial P_{w}}{\partial V_{i_{Im}}}=\gamma_{i}\zeta_{2i},\frac{\partial Q_{w}}{\partial V_{i_{Re}}}=\gamma_{i}\zeta_{3i},\frac{\partial Q_{w}}{\partial V_{i_{Im}}}=\gamma_{i}\zeta_{4i},i=n+1,\dots,n+p$. Accordingly, the admittance matrix $\bar{A}_{3}$ is expressed as

[TABLE]

where $A_{3}^{i}$ has the same structure as $A^{\prime}_{3}$ in (22) but with $\zeta_{1},\zeta_{2},\zeta_{3},\zeta_{4}$ replaced by $\zeta_{1i},\zeta_{2i},\zeta_{3i},\zeta_{4i}$, respectively. The comparison between $\mathcal{L}$ and $\mathcal{L}_{0}$ follows thereafter in the same way as in the foregoing subsection.

V Perturbation analysis for two time-scale property

We recall from (11), that $\mathcal{L}_{0}=\mathcal{L}_{0}^{I}+\epsilon\mathcal{L}_{0}^{E}$. We next analyze how the perturbation in $\mathcal{L}_{0}$ as given by (37) extends to its the internal and external connection components $\mathcal{L}_{0}^{I}$ and $\mathcal{L}_{0}^{E}$, respectively.

V-A *Perturbation in $\mathcal{L}_{0}^{I}$ and $\mathcal{L}_{0}^{E}$ *

Due to the existence of $r$ clusters one can separate $A_{3}$ into internal and external admittance matrices as, $A_{3}=A_{3}^{I}+\epsilon A_{3}^{E}$. Using exponential expansion we can write,

[TABLE]

where $X_{1\epsilon}=[(-(A_{3}^{I})^{-1}(A_{3}^{E})+\epsilon((A_{3}^{I})^{-1}(A_{3}^{E}))^{2}-\dots)(A_{3}^{I})^{-1}]$. Considering a single wind farm scenario we recall (37) as

[TABLE]

The matrix $X$ can be written as

[TABLE]

where $x=(\gamma A^{\prime}_{3})+\Delta_{A3}$. Expanding $(I+A_{3}^{-1}x)^{-1}$ we get,

[TABLE]

where $P_{1a}=I+(A_{3}^{I})^{-1}x$, and $X_{2\epsilon}$ has the same structure as $X_{1\epsilon}$ with $A_{3}^{I},A_{3}^{E}$ replaced by $P_{1a}$ and $X_{1\epsilon}x$, respectively. In that case, we have

[TABLE]

Here, $P_{a}=P_{1a}^{-1}(A_{3}^{I})^{-1}x(A_{3}^{I})^{-1}$ and $P_{b}=X_{2\epsilon}((A_{3}^{I})^{-1}+\epsilon X_{1\epsilon})x((A_{3}^{I})^{-1}+\epsilon X_{1\epsilon})+P_{1a}^{-1}X_{1\epsilon}x((A_{3}^{I})^{-1}+\epsilon X_{1\epsilon})$. Then finally (25) can be rewritten as

[TABLE]

where $\Delta_{\mathcal{L}^{I}}=K_{12}P_{a}A_{1}+\kappa_{\mathcal{L}}(\gamma),\Delta_{\mathcal{L}^{E}}=K_{12}P_{b}A_{1}$. These perturbations, which are functions of $\gamma$ and the wind farm location, quantify the changes in the internal and external components of $\mathcal{L}_{0}$.

V-B Extraction of the equivalent Laplacian

We next substitute the expression (46) in $\mathcal{R}_{1}=M^{-1}\mathcal{L}$ in the wind-integrated model (32). From this model we can write

[TABLE]

Our intent is to capture the interactions between the synchronous machines in the wind-integrated system. For that we apply the transformation (7) on (47), which results in the following transformed unforced dynamics

[TABLE]

where,

[TABLE]

Here $\tilde{q}_{s},\tilde{q}_{f}$ are the slow and fast states corresponding to synchronous-only motions. Due to space constraints we skip the derivations of (49)-(52). Comparing the nominal coherency dynamics in (8) with that of the perturbed dynamics in (48), it is clear that if the matrix $\tilde{T}$ in (48) still has to reflect the coherency between synchronous generators in the wind-integrated model, then we must redefine $\mathcal{L}$ such that $CM^{-1}\Delta_{\mathcal{L}^{I}}G^{\dagger}=\hat{M}U^{T}\Delta_{\mathcal{L}^{I}}G^{\dagger}=0$ and $\Delta_{\mathcal{L}^{I}}U=0$. This can be done simply by defining an equivalent Laplacian matrix in the following way:

[TABLE]

This definition allows us to write

[TABLE]

where now both perturbations $\Delta_{\mathcal{L}_{eq}^{I}}$ and $\Delta_{\mathcal{L}_{eq}^{E}}$ are Laplacian matrices since $\mathcal{L}_{0}^{I}$ and $\mathcal{L}_{0}^{E}$ are Laplacian matrices by default. This satisfies $\hat{M}U^{T}\Delta\mathcal{L}_{eq}^{I}G^{\dagger}=0$ and $\Delta\mathcal{L}_{eq}^{I}U=0$, making $\tilde{T}$ structurally consistent with $T$.

Once constructed, the equivalent Laplacian matrix $\mathcal{L}_{eq}$ is next used to identify the coherent groups of synchronous generators in the system. This can be done via Algorithm $1$ with $\mathcal{L}$ replaced by $\mathcal{L}_{eq}$. Let $\tilde{V}$ consists of eigenvectors corresponding to smallest eigenvalues of $M^{-1}\mathcal{L}_{eq}$. The matrix $\tilde{V}_{r1}$ is constructed by aggregating the rows corresponding to reference machines obtained using Gaussian elimination. Permuting the rows of $\tilde{V}$ we will have $\tilde{V}_{r}=[\tilde{V}_{r1}^{T}\;\;\tilde{V}_{r2}^{T}]^{T}.$ Then the row vectors of $\tilde{V}_{r1}$ are used as unit coordinate vectors in the new coordinate system using the transformation

[TABLE]

Depending on the difference between $\mathcal{L}_{0}$ and $\mathcal{L}_{eq}$, the matrix $\tilde{L}$ may now produce a different clustering structure than the wind-less system. Note that even if the change in the individual entries of the $\mathcal{L}_{eq}$ from $\mathcal{L}_{0}$ may not be significant, its slow eigenspace can still be notably perturbed so that the coherent groupings before and after wind penetration are different, shown in the following simple example. Let $\alpha_{i},\tilde{\alpha}_{i},(i=1,\dots,n)$ be the row vectors of the nominal and perturbed eigen-spaces $V_{L}=[I,L^{T}]^{T}$ and $\tilde{V}_{L}$ following from $\mathcal{L}_{0}$ and $\mathcal{L}_{eq}$, respectively. From [12] it follows that these perturbed row vectors will also lie on the hyperplane $\sum_{j=1}^{r}\tilde{\alpha}_{ij}=1,\forall i=1,..,n,$ just like the nominal row vectors. Fig. 3 shows the visualization of this hyperplane for a $2-$area $4$-machine system in $2$-dimensional coordinates where $\alpha_{i}$ and $\tilde{\alpha}_{i}$ are plotted for $i=1,2$. The tip of both row vectors will remain on the same hyperplane as shown in the figure, but their positions may shift depending on $\gamma$ and the wind bus location, which, in turn, can change the coherent grouping.

VI Simulation Results

We perform simulations on a $16$-machine, $68$-bus IEEE benchmark power system model - first with a single wind farm and then with three wind farms. This is a simplified model of interconnected New York and New England power systems. Simulations have been performed in Matlab $2016a$ platform with intel(R) core(TM)-$i7-2600$ CPU$@3.40$ GHz processor. The total active load is $17.6$ GW. Each wind turbine-generator unit is rated at $1.76$ MW and wind parameters are taken from [10]. The internal PI controllers in the active and reactive power loops of the DFIGs are tuned to achieve stable operation.

VI-A Nominal system without wind plant

The slow oscillation modes of the system without any wind plant and the corresponding coherency structure using Algorithm $1$ are shown in Table-I. The areas are marked in the system diagram in Fig. 6(a). The orientation of row vectors of the slow eigenspace are shown in Fig. 4. The coherency algorithm assigns machines $5$ and $13$ as reference machines for Area $1$ and Area $2$, respectively.

VI-B Single wind plant connected to the grid

Several subcases are considered with the wind plant connected to buses $66,37,32$, and $38$. Bus location $66$ is considered as it is away from loads and generation zones in Area $1$. Bus $37$ is tested because of its close proximity to the loads in this area. Total connected load in Area $1$ is more than the total generation. Bus locations $32$ and $38$ belong to Area $2$, where total connected load is less than the generation. Here we increase wind penetration to higher values in order to test the change in coherency behavior.

We first present the case where the wind plant is located at bus $66$ with $\gamma=650$ ($1144$ MW). The blue and yellow stems in Fig. 5 denote the values of the entries of $\mathcal{L}_{0}$ and $\mathcal{L}_{eq}$, showing how the row entries of the Laplacian corresponding to generator $1$ change before and after wind penetration. The wind-integrated model has four slow modes as shown in Table-II. We construct $\mathcal{L}_{eq}$ and apply algorithm 1 to identify the five slow coherent areas arising from these four slow modes. These areas and the indices of the synchronous generators in each area are also shown in Table-II. Comparing tables I and II, it can be clearly seen that the wind injection forces generators $1$ and $8$ to move from Area 1 to Area 2. The other generators remain in their respective areas. Fig. 6 shows the system diagram with modified clustering before and after wind injection. Fig. 7 shows the row vectors of the perturbed slow eigen-space $\tilde{V}$ while Fig. 8 shows the compass plots for the selected row vectors of the transformed eigen-space $V_{L}$ of the nominal system (Fig. 8a) versus $\tilde{V}_{L}$ of the perturbed system (Fig. 8b). Both figures testify to the tendency of generators $1$ and $8$ to move from Area 1 to Area 2.

Next we consider the case when the wind plant is located at bus $37$ with $\gamma=700$ ($1232$ MW). Table-III shows the slow modes and the corresponding coherent clusters. Compared to the previous case the generator indices in the respective areas do not change; however, the reference machine of Area $1$ now changes from $6$ to $5$.

We further validate our results for this case using the model-free PCA technique. The comparison is done only for the validation of our derivations, and not intended for comparison between model-driven and data-driven methods. PCA is carried out on the rotor angle measurements of all $16$ generators for $100$ seconds with a sampling time of $0.01$ s, following a unit step change in the mechanical power of machine $1$. Fig. 9 shows the first three principal components with the maximum variance in the data set. When the wind plant is connected to bus $37$ with $\gamma=700$ the plot clearly shows that the positions of generators $1$ and $8$ move more towards generators $11$ and $12$. The same is predicted by the coherency algorithm using $\mathcal{L}_{eq}$.

The third case considers the wind plant connected at bus $32$ in Area $2$. For this case even when $\gamma=700$, the coherency structure was not found to change from the nominal. Similar observation is made when the wind plant is located at bus $38$. This indicates that Area $1$ is more prone to perturbation in coherency than Area $2$. This can be a useful message for transmission planners for deciding the location of wind installations without disturbing the coherency of their grid.

VI-C Wind plants at buses $32,66$ and $57$

We next consider the IEEE $68$-bus system with three wind farms located at buses $32,66$ and $57$ with penetration levels $\gamma_{1}=200$, $\gamma_{2}=250$, $\gamma_{3}=200$, respectively. The corresponding $\mathcal{L}_{eq}$ is constructed and Algorithm 1 is applied. The resulting clustering structure is shown in Table-IV. The table shows that with wind installed at these three locations, Area 1 now shrinks to only two generators - namely, generators $2$ and $3$. Area $2$, on the other hand, now expands to a much bigger geographical area covering a total of $11$ generators. Areas $3$ through $5$, however, remain unaffected. The row vectors of $\tilde{V}$ are plotted in Fig. 11 indicating the same result. The result is also validated using PCA, as shown in Fig. 11, where generators $2$ and $3$ depart from their nominal dynamic signatures forming an area between just the two of them.

VII Conclusion

A mathematical analysis of the perturbation in coherency of synchronous generators due to wind integration is presented in this paper. The dynamic interaction between the generators in a wind-integrated system is captured by an equivalent Laplacian matrix. Depending on the amount of wind injection and placement of wind plants the slow eigenspace of the equivalent Laplacian matrix may change, thereby changing the coherent groupings. Results are validated using the IEEE $68$-bus system with single and multiple wind farms. The results can be useful for transmission planners in deciding potential locations of wind installations, and also for readjusting wide-area control gains in case the wind injection changes the coherent groupings.