Wireless Multiple Access and Analog Oscillation Analysis in Microgrids: A Prony Approach

A key challenge in smart grid is to monitor the system modal characteristics of disturbance over a microgrid for system stability and diagnosis. The measurements at different phasor measurement units (PMUs) need to be reported to a processing center. Therefore, a multiple access scheme for wireless transmission is necessary to collect the measurements from PMUs. An analog oscillation analysis is carried out without demodulating the received signal. Both schemes of time-division multiple access (TDMA) and code-division multiple access (CDMA) are discussed. In both schemes, the Prony-based algorithms using the received signal to estimate the system modal characteristics, including the eigenvalues of the linear system and the amplitudes of different response modes, are proposed and analyzed. Two countermeasures, max-min ratio measure (MMRM) and variance measure (VM), are also proposed to prevent possible attacks of false reports. Numerical simulation results demonstrate the effectiveness of Prony-based algorithms and the two countermeasures.


Introduction
In recent years, smart grid has attracted many studies in the community of power system, communications, and networking.According to the Energy Independence and Security Act of 2007 (EISA), smart grid is defined as follows: "A modernization of the Nation's electricity transmission and distribution system to maintain a reliable and secure electricity infrastructure that can meet future demand growth" [1].An important component in smart grid is the microgrid [2,3], which is a subsystem containing multiple distributed energy generators (DEGs, such as microturbines, solar panels, or wind turbines) and multiple loads.A microgrid can work either together or independently with the power grid.In normal situations, the microgird is connected to the power grid and transfers energy either to or from the major power grid according to its own energy generation and consumption.When the power grid is unstable, for example, in a large area blackout, the microgird can be disconnected from the power grid and work in an island mode.At this time, its own DEGs will support the power load within the microgrid.Such a microgrid has been implemented in many testbeds, for example, the Electric Reliability Technology Solutions (CERTS) Microgrid supported by the US Department of Energy.
One problem impacting the reliability of power system is power system oscillations.Power system oscillations, caused by either network change or disturbance, appear as modes which are complex exponentially modulated sinusoids.Insufficiently damped modes may incur uncontrolled separation of the power system into islands, consequently resulting in blackouts.Therefore, for the reliability of a power system, it is important to estimate the power system state, detect and identify oscillation modes accurately and as early as possible.Oscillation exists not only in large power grids but also in microgrids.It may be even more severe in a microgrid than in a large power grid.However, the research of oscillations is very limited in microgrids.Therefore, it is important to study power system oscillation in microgrids for achieving a better reliability.
Currently, mode identification of power system oscillation is well studied in power grid.With the deployment of Phasor Measurement Unit (PMU), power oscillation detection and identification have been widely implemented in Wide-Area Monitor System (WAMS) [4,5].Typically, the system modal characteristics are estimated from the measurement of several synchronized PMUs.After taking the reports from the PMUs, a processing center will analyze the data and estimate the system modal characteristics, for example, the eigenvalues of the transfer function, as well as the amplitude of different components with different dampening ratios.If instability is detected, the system can take certain actions to stabilize the power grid.An important aspect of microgrid is the communication infrastructure for monitoring the oscillation due to the requirement of real-time processing and high reliability.In such a communication infrastructure, the observations at PMUs are sent to a processor.One standard related to the communications of PMUs is IEEE C37.118 [6].However the communication protocol defined in IEEE C37.118only specifies a few message types with few options and allows the usage of any possible communication approaches for the PMUs to report to the processing center, for example, wireless communications or optical fibers.In traditional microgrids, wired communications are employed.For example, an RS-485 (digital computer based) link or an ethernet is used in the Oak Ridge National Lab (ORNL) microgrid system [7].An alternative approach for the communication for DEG control is to use wireless communication systems due to its fast deployment.For wireless communication systems, the multiple access as shown in Figure 1 is a key challenge to the system monitoring since there exist several wireless PMUs reporting their observations.Since the data traffic is intensive, it is improper to use random access schemes like carrier sense multiple access (CSMA), which might incur too many collisions.Depending on whether collisions are allowed or not, we can consider orthogonal or nonorthogonal multiple access schemes.
(i) Orthogonal: for simplicity, we consider only the time-division multiple access (TDMA) scheme.The corresponding system estimation is straightforward.The challenge is how to deal with the time offset of reports from different PMUs due to the TDMA.
(ii) Nonorthogonal: we consider the code-division multiple access (CDMA), which allows collisions.One traditional approach is to first recover the original reports from the mixed received signal.However, this may not be the most efficient approach.We will study the approach of estimating the system state directly from the received signal.
In this paper, we will study both approaches.The power system will be modeled as a linear system with the assumption of small signal perturbations.The system modal characteristics will be estimated using Prony's algorithm [8,9], which will be made compatible with the TDMA and CDMA schemes.The Cramer-Rao-bound (CRB) covariance matrix for both TDMA-and CDMA-based Prony methods will be derived.Both TDMA and CDMA approaches will be compared using numerical simulations.The conclusions will provide insights for the design of communication infrastructure in microgrids and advance the signal processing techniques in the system modal characteristics estimation of microgrids.Moreover, we consider the case in which there exists an attacker and discuss how to detect the attack when attacking happens.We propose two variance-based countermeasures to identify the attack.
The remainder of this paper is organized as follows.Related work for modal identification of power system oscillation will be provided in Section 2. The system model will be introduced in Section 3. The multiple access schemes of TDMA and CDMA will be discussed in Section 4. The CRB covariance matrices for both TDMA-and CDMA-based Prony methods are derived in Section 5.The countermeasures for possible attacks are proposed in Section 6.Then, we will show the numerical simulation results in Section 7. Finally, the conclusion will be drawn in Section 8.

Related Work
Parametric modal identifications are widely used in power systems.One main method is based on eigenvalue analysis.However, the disadvantage of eigenvalue analysis is that it requires the knowledge of the system matrix.The difficulty of monitoring the whole state space system matrix prevents it from online tracking of power oscillation modes.In contrast, low-order linear estimations are suitable for online tracking of power oscillation modes.The advantage of these methods is that they do not require the linearization of the system and the associated large-scale system matrix.The low-order linear estimation methods include Steiglitz-McBride Algorithm (SMA) [10], Eigensystem Realization Algorithm (ERS) [11], and Prony method [12].SMA uses an iterative procedure to adjust the coefficients of a transfer function; ERA is based on singular value decomposition of the Hankel matrix associated with the impulse response of the system; the Prony method is based on a linear prediction model, which estimates parameters to fit a finite number of samples of a signal.In [13] the performances of these three algorithms are well compared.Numerical results show that Prony method has more accurate outputs for various inputs.Therefore, in this paper we focus on the Prony method for the mode identification of power system oscillation in microgrids.
The Prony method has been widely applied in power systems for mode identifications [12,[14][15][16][17][18][19][20].Currently the online Prony analysis has been deployed for operation in many places, for example, Jiangsu in China [4] and the Western Electricity Coordinate Council [5].The first application of Prony method for power system oscillation is [12,14] which used field measurement data.The Prony method proposed in [12,14] is based on the least square error criterion and processes each signal independently.However, Prony method is sensitive to noise.The Prony analysis of individual signals may have conflicting frequencies and damping coefficients.References [8,21] have extended the Prony method to analyze multiple signals simultaneously in order to improve the accuracy.In our paper, we also combine multiple signals to improve the accuracy of Prony analysis.The contribution of our paper is that the system state is directly estimated from the received signal without recovering the original information bits, which is coined analog oscillation analysis.References [22,23] derived the CRB covariance matrix of Prony method with single and multiple signals, respectively.[24,25] derived the covariance matrix for singular-value-decomposition-(SVD-) and totalleast-square (TLS) based Prony method, respectively.In our paper, the CRB covariance matrices for TDMA-and CDMAbased Prony methods are derived.

System Model
We model the power grid as a linear system, whose dynamics are given by [9] ẋ = Ax + bu, where x is an N-dimensional system state vector, Y = [y 1 , y 2 , . . ., y M ] is an observation vector for M PMUs, u is a scalar control signal, and n is the noise added to the observation.Each dimension of the observation is sensed by a PMU and sent to a processing center via the communication infrastructure.Note that power system is usually highly nonlinear.However, the linear model in (1) approximately holds for the small-signal case; that is, the perturbation to the signal is small.
Consider the mth output y m , m = 1, . . ., M, and the input u.The corresponding transfer function is given by where λ n is the nth eigenvalue of matrix A, which represents an oscillation mode, and R m n is the amplitude of the nth oscillation mode.
The observation at PMU m is then given by Suppose that each PMU samples the continuous-time observation with time step Δt.Then, the samples at PMU m are given by where z n = exp(λ n Δt).Note that the eigenvalues {λ n } n and the amplitudes {R m n } n,m are unknown; the task of system modal characteristics monitoring is to estimate the eigenvalues and amplitudes from the observations {y m (k)} m,k .
We suppose that the communication system for the monitoring system has a bandwidth W. The time is divided into time chips; each has a time interval T c approximately equaling 1/W.We suppose that each PMU can transmit one observation during one time chip.For simplicity, we ignore the quantization error.At the tth time chip, the received signal at the processing center in the base band is given by where a m is the channel amplitude gain, which is assumed to be time-invariant, h m (t) is the signal transmitted by PMU m at time chip t, and o(t) is noise, the variance of which is σ.

Multiple Access Schemes
In this section, we discuss two types of multiple access schemes, namely, the TDMA and CDMA.We assume that the reports from different PMUs are honest.Usually, the observations at the PMUs have a very high signal-to-noise ratio (SNR).We also assume that the communication channel quality is very good and has a high SNR.Therefore, we ignore the noise in the discussion in this section.The impact of noise will be evaluated in the numerical simulations where the algorithms with the assumption of no noise will be applied in the noisy case.

TDMA.
In TDMA, we simply let the PMUs report their measurements in a fixed order.We can set Δt = MT c ; that is, one sample per M time chips.In the time chip scheduled for PMU m, it transmits the newest measurement in this time chip.In other time chips, PMU m does not transmit.Therefore, the transmit symbol h m (t) at time chip t, in the general form (5), is given by Upon receiving MW time chips of reports, the received signal at the processing center can be written as where the observation matrix Y is given by and the matrix Z T is given by where z n = z M n , for all n = 1, . . ., N, and the amplitude matrix R T is given by We can consider the whole quantity of R m n z m−1 n as the unknown.Once we obtain R m n z m−1 n and z n from the observations, we can easily obtain R m n , thus removing the effect of time offset.
To solve the equation for the unknowns {R m n z m−1 n } n,m and { z n } n , we adopt the Prony's method [8,9].Here, we provide a brief introduction to the Prony's method, which will be modified for the CDMA scheme.The quantities { z n } n=1,..., N are necessarily the roots of an N-order polynomial with coefficients {a} n (unknown to the processing center), which is given by Then multiplying vectors in the form of (−a N , . .., −a 1 , 1, 0, . . ., 0) on the left side of (7), we have Then, the Prony's method includes the following four steps.
(1) The coefficient vector {a n } n is estimated by solving the previous overdetermined least square problem based on the assumption that (W − N)M is larger than 2N.Note that there are very well-known solutions to the least square problem [26].Many solutions are applied to the Prony method to reduce the impact of noise [4,24,25].(2) Solving the polynomial equation ( 11) for { z n } n .Note that ( 11) is a Nth order polynomial equation after the coefficient vector {a n } n is estimated in the above step.(3) Calculating eigenvalues {λ n } n from { z n } n ; that is: (4) Solving (7) for the amplitudes {R m n } n,m .(5) Order the eigenvalues { z n } n based on estimated {R m n } n,m to get the main eigenvalues.A widely used criterion to search for the optimal subset of eigenvalues that best represent the system is the Akaike information criterion (AIC) [27].

CDMA.
In the CDMA scheme, collisions among the PMUs are allowed.We define a symbol period as T s time chips during which each PMU transmits one observation.During the kth symbol period, PMU m transmits the T svector y m [(k − 1)T s + 1][s m (1), . . ., s m (T s )] over the T s time chips.The sequences {s m (t)} m,t are called signature waveforms, which are known to the PMUs and the processing center in advance.Note that we fix the signature waveform for all symbol periods, which facilitates the future processing.Then, the received signal, a T s -vector, during the kth symbol period is given by where Upon receiving W symbol periods and putting the received signal together, the received signal is given by where D = [r(1), . . ., r(W)] and Y = [y(1), . . ., y(W)].
According to (4), we have where and Z C has the same definition as Z T in ( 9) except z n = z Ts n .Combining the previously mentioned definitions and expressions, we have One approach is to use a CDMA demodulator to compute Y first and then apply the Prony's method for the decomposed signal from each PMU, as illustrated in Figure 2(a).
However, by observing the similarity between ( 7) and (19), we can propose a new algorithm as shown in Figure 2(b) to apply the Prony's method again as follows.
(1) Consider term RS T in (19) as matrix R T in (7).Apply the first step in the Prony's method in Section 4.1 to obtain the coefficients (a 1 , . . ., a N ).
(3) Then, we use (19) to compute the amplitude components in R. To that end, we apply the equation, which holds for any matrices A, B, and C with suitable sizes: where the operator vec means stacking the columns of a matrix to obtain a vector and ⊗ is the Kronecker produce [28].Then, we have Finally, the amplitudes in matrix R are obtained by solving the above linear equations.

Cramer-Rao-Bound Covariance Matrix for Multiple Access Scheme
In this section we derive the CRB covariance matrix for estimating {R m n } and {z n } in TDMA-and CDMA-based Prony method.CRB is a fundamental lower bound to the variance of any unbiased parameter estimator.The derived CRB can not only help estimate the impact of noise on the estimation of {R m n } and {z n } but also help evaluate the accuracies of TDMA-and CDMA-based Prony methods while considering the impact of noise.

TDMA.
In TDMA systems, the received signal from the mth PMU can also be written as where Then, the likelihood and log-likelihood of data from all M PMUs are given by L(r(1), . . ., r(WM)) where In order to calculate the CRB covariance matrix, the following parameters are defined.The real and imaginary parts of R m n are denoted by R m n and R m n , respectively.The real and imaginary parts of z n are denoted by z n and z n , respectively: Then the Fisher information matrix I θ can be calculated by taking partial derivatives of (25) with respect to θ which is shown in (28), where The proof of ( 28) can be found in Appendix A.
International Journal of Distributed Sensor Networks 7 Then the CRB Covariance matrix CRB θ is calculated by taking the inverse of the Fisher information matrix I θ shown in (31), where The corresponding proof can be found in Appendix A.

CDMA.
In CDMA systems, the received signal from all M PMUs for time slots (w − 1)T s , . .., wT s can also be written as where Then the likelihood and log-likelihood of data from all M PMUs are given by and where Then the Fisher information matrix I θ can be calculated by taking partial derivatives of (36) with respect to θ which is shown in (38), where The corresponding proof can be found in Appendix B.  Then the CRB covariance matrix CRB θ is calculated by taking the inverse of Fisher information matrix I θ shown in (41), where The proof is the same as that for (31).

Countermeasure for Attacks
In the previous procedure of system modal characteristics estimation, we implicitly assume that all PMUs are honest and report their observations faithfully.However, in practice, a PMU could be compromised by a malicious party.It is also possible that the sensor at the PMU malfunctions and gives wrong readings.The false report could result in a significant error in system estimation or could even incur power system instability as a result of the system taking an incorrect control action due to the wrong state estimation.In this section, we study the countermeasure for detecting possible attacks.

Attack Model.
Attacks are defined as all events which generate data including false or faked mode characteristics.We notice that the amplitudes {R m n } have no redundancies in the received signal since it is completely determined by the reports from PMU m.Therefore, it is impossible to detect the attack on the amplitudes unless we introduce more redundancies, which is beyond the scope of this paper.Therefore, we focus on only the attack on the eigenvalues, which is more important than the amplitudes since the eigenvalues directly determine the stability of the system.Meanwhile, we assume that all the PMUs share the same set of eigenvalues, thus resulting in the redundancy and facilitating the detection of attacks.In addition, we assume that attackers work independently; that is, there is no collaborative attack.In other words, the false or faked eigenvalues generated by attacks are different from each other.For simplicity, in this paper we consider the case with one attacker.

Attack Detection.
The procedure of detecting the attack and attacker is shown in Figure 3.As shown in the procedure of the Prony analysis in Section 4, the normal PMU and compromised PMU share the same set of estimated eigenvalues z i .The difference between normal PMU and compromised PMU is the estimated amplitudes R m n .For instance, suppose that PMU k and PMU l are PMUs without attack and with attack, respectively.z i is one of the eigenvalues of normal PMU k; however, it is not the eigenvalue of compromised PMU l.Then, the corresponding estimated R l i for attacked PMU l may be small with a high probability.In other words, the estimated R k i of normal PMU is greatly larger than that of the attacked PMU l with a high probability.For paths without attack, although R i i and R j i of PMU i and j for eigenvalue z i are different, the difference between R i i and R j i may not be very large.Based on this observation, we propose two countermeasures to detect the attack, namely, max-min ratio measure (MMRM) and variance measure (VM).

MMRM Approach.
MMRM is based on the ratio of maximum and minimum absolute amplitudes of R(i, m) in row i corresponding to eigenvalue z i ; that is: where N f is the number of eigenvalues which will be calculated in the measure and R(i, j) is the estimated coefficient of ith mode for PMU j.Note that we assume that the estimated eigenvalues are already sorted based on the AIC criteria.Figure 4 demonstrates the MMRM of normal data and attacked data for TDMA when SNR is equal to 30 dB.For both attacked data and normal data, there are three PMUs.The first two PMUs are normal, and the third PMU is attacked.The detailed expression of signals will be provided in Section 7.2.As shown in the Figure 4, MMRM for attacked data is significantly larger than that of normal data, and there is a gap between MMRM for attacked data and normal data.Therefore, MMRM can be used as a measure to detect the possible attack.

VM Approach.
VM is based on the variance of each row of R: where R(i, :) is the ith row of R and Var(X) is the variance of X. Figure 5 illustrates the VMs of attacked data and normal data for TDMA when the SNR is equal to 30 dB.The signals are the same as those for MMRM.As shown in Figure 5 there is a gap between the VMs for normal data and attacked data.Therefore, VM can also be used as a measure to detect the attack.

Attacker Detection.
Once the processing center finds that an attack has been launched, the next job is to detect the attacker from the received signal.One simple approach is to separate the observations from different PMUs and then use the Prony method to analyze the observation of each PMU separately and estimate eigenvalues of each PMU.Then, it is easy to detect the attacker as the estimated eigenvalues of attacked PMU are different from those of normal PMUs.Note that this detection procedure needs to execute the Prony's method for every PMU which is very time consuming.Another possible solution is to use MMRM or VM based cluster analysis to classify the signals into different clusters.
Considering the fact that in most cases the MMRM or VM among normal data is smaller than that between normal data and attacked data, the MMRM or VM based cluster analysis method puts the attacked data in the cluster that does not include normal data.Therefore, when there is only one attacked PMU, the signal in the cluster with only one element is identified as the attacked data.There are many papers (including those we have surveyed [29][30][31]) that discuss the cluster analysis.Due to the limited length of the paper, the detailed procedure for cluster analysis is not provided.

Numerical Simulations
In this section, the simulation results of proposed Prony analysis for TDMA and CDMA are provided.Then, the results of proposed MMRM and VM for attack detection for TDMA are shown.

Performance of Prony Analysis for TDMA and CDMA.
In the simulation, we assume that there are two PMUs.The signals are given by International Journal of Distributed Sensor Networks  The sampling interval T s is 0.1s, the total number of samples for each PMU is 100, and the spreading matrix 6 and 7 demonstrate the mode distribution of TDMA and CDMA, respectively, when the SNR is equal to 20 dB.There are three modes both PMUs: −0.3265 + 1.9604i, −0.4143 + 3.2924i, −0.5287 + 7.6215i.For convenience, these three modes are labeled as modes 1, 2, 3, respectively.As shown in the figures, for both PMU 1 and 2, the estimated modes for TDMA and CDMA are closer to the real modes than that of single PMU.
In order to further compare the performance of Prony analysis for TDMA, CDMA, and single PMU, the mean square error (MSE) of Prony analysis is illustrated in Figure 8.Here, MSE is defined as the average square error between the estimated modes and the real modes.As shown in the figure, the performance of both TDMA and CDMA is better than that of either PMU 1 or 2. Therefore, the Prony analysis for TDMA and CDMA has performance gain from multiple signals and has more capability to combat the noise, compared with Prony analysis for single PMU.

Countermeasure for Attacks.
As the performance of MMRM or VM of CDMA is similar with the result for TDMA, only simulation results for TDMA are provided in this part.Besides signals f 1 and f 2 , another correct signal f 3 and an attacked signal h 3 are added in the simulation.f 3 has the same eigenvalues as signal f 1 and f 2 .h 3 has different eigenvalues: (48) Figure 9 shows the performance of MMRM with different N f when the SNR is 20 dB.When N f = 3, the detection probability increases rapidly to 1 when false alarm probability is 0.05.The performance degrades gradually as N f increases.Figure 10 compares the performance of different SNRs when N f = 3.When the SNR is equal to 25 dB, the performance is very good.The detection probability is almost 1 with different false alarm probabilities.When the SNR is equal to 20 dB, the detection probability rises quickly to 1 when the false alarm probability is larger than 0.05.
Figure 11 shows the performance of VM with different N f when the SNR is 20 dB.As illustrated in the figure, except for N f = 3, the performance of other N f is good; the detection probability is almost 1 when the false alarm probability is greater than 0.1.Figure 12 compares the  performance of different SNRs for N f = 3.As shown in the figure, when SNR is equal to 25 dB, the performance is very good.The detection probability is very close to 1 with a low false alarm probability.When the SNR is equal to 20 dB, the detection probability is near 0.9 when false alarm probability is 0.1.

Conclusions
In this paper, TDMA and CDMA are proposed to transmit measurements of different PMUs to a processing center.In both schemes, the Prony algorithms using the received signal to estimate the system state, including the eigenvalues of the linear system and the amplitudes of different response modes, are proposed and analyzed.The simulation results have demonstrated that the performance of Prony analysis The derivatives of ln L in (25) with respect to z n and z n are given by Then, the derivatives of ln L with respect to z and z can be calculated, which are given by For the simplification of the expression for Fisher information matrix I θ , we assume E e(m)e T (s) = 0, for ∀m, s. (A.8) In order to obtain the final expression for I θ , we also need to use the following property equations which are proved in [22]: Then, the elements of I θ can be obtained via This concludes the proof of (28).
For the proof of (31), we need to use the following lemma which was proven in [22,23].Based on Lemma 1 and the property of partitioned matrix, the inverse of I θ can be directly calculated and (31) is proved.

B. Derivation of CRB Covariance Matrix for CDMA
The derivative of ln( Ľ) in (36) with respect to σ is given by ∂ ln Ľ In order to attain the final expression for Ǐθ , besides using (A.8), (A.9), and (A.10), we also need to use the orthogonal property of S; that is: where I is an identity matrix.
Then, the elements of Ǐθ can be obtained:

Figure 1 :
Figure 1: An illustration of the multiple access of PMUs.

Figure 2 :
Figure 2: Two approaches for the system estimation with CDMA.

Figure 3 :
Figure 3: An illustration of the detection of attack and attacker.

Figure 4 :
Figure 4: An illustration of the max-min ratio measure for Prony analysis for TDMA, SNR = 30 dB.

Figure 5 :
Figure 5: An illustration of variance measure for Prony analysis for TDMA, SNR = 30 dB.

T
CZ T T z (D z ) m−1 T R (m)

Lemma 1 .
Assuming that one is given a block matrix with the following form:B r = B − B B B , (A.22)where B is a b × b block matrix and B and B are the real and imaginary parts of B, then the inverse of B r is equal to(B r )where B = (B) −1 , and B and B are the real and imaginary parts of B .