Grid Topology Identification using Electricity Prices
Abstract
The potential of recovering the topology of a grid using solely publicly available market data is explored here. In contemporary whole-sale electricity markets, real-time prices are typically determined by solving the network-constrained economic dispatch problem. Under a linear DC model, locational marginal prices (LMPs) correspond to the Lagrange multipliers of the linear program involved. The interesting observation here is that the matrix of spatiotemporally varying LMPs exhibits the following property: Once premultiplied by the weighted grid Laplacian, it yields a low-rank and sparse matrix. Leveraging this rich structure, a regularized maximum likelihood estimator (MLE) is developed to recover the grid Laplacian from the LMPs. The convex optimization problem formulated includes low rank- and sparsity-promoting regularizers, and it is solved using a scalable algorithm. Numerical tests on prices generated for the IEEE 14-bus benchmark provide encouraging topology recovery results.
uclear norm regularization; compressed sensing; alternating direction method of multipliers; economic dispatch; locational marginal prices; graph Laplacian.
1 Introduction
Data analytics is a major driver towards the smart grid transformation: from load and wind forecasts, to consumer preference learning, and cyber-physical attack detection. Among grid data mining tasks, topology inference is one of vital importance. Currently, system operators and utilities maintain a detailed physical model for the underlying transmission and distribution grids. The model is updated regularly by the network topology processor, and constitutes the foundation for several monitoring and market-related tasks; see e.g., [1].
On top of these conventional cases, identifying the grid topology can be useful for various purposes. A well-designed cyber attack on the state estimator requires grid topology information. Knowing the power network structure could assist in informed bidding strategies, or even lead to market manipulation; see e.g., [2]. In addition, the Laplacian matrix of the graph corresponding to the grid provides a meaningful notion of inter-bus electrical distances. Such distances could surrogate the spatial correlation among pricing nodes [3]; or adopted by clustering techniques to reveal influential nodes.
Albeit the extensive research on data attacks, grid topology recovery using readily accessible data has been overlooked. The possibility of arbitrarily perturbing power system state estimates by counterfeiting a few measurements is recognized in [4]. The strength of such data attacks and their impact on market outcomes is characterized in [5], [6]; yet designing stealth attacks oftentimes assumes that the attacker knows the topology of the grid [7]. From the system operator’s point of view, detecting topological changes in the form of transmission line outages has been studied in [8], [9], [10]. Given hourly topology updates and nodal voltage phases, the sparse overcomplete representation employed in [10] can reveal outages of even geographically distant lines at affordable complexity. In a microgrid scenario, the timing difference of power line communication signals to travel between electrical nodes is utilized to recover the grid structure in [11]. Grid topology recovery is cast as a blind identification problem in [12]. Although [12] exploits the positive semidefiniteness and the sparsity of the associated weighted grid Laplacian, it presumes that bus voltage phases are linearly independent and nodal injections are available across the entire interconnection. Under a Markov random field (MRF) assumption, transmission network faults are localized using the sample covariance matrix of nodal voltage phases in [13]. Likewise, if nodal voltage magnitudes constitute an MRF, the topology of a distribution grid is pursued by means of the sample covariance matrix of voltage magnitudes [14].
Distinct from the previous setups, this work considers grid topology recovery using publicly available market data. Real-time locational marginal prices (LMPs) are calculated based on the Lagrange multipliers of the network-constrained economic dispatch (Section 2). The fresh idea here is that LMPs could unveil grid topology information. Specifically, let be the matrix collecting the real-time LMPs announced over a network of nodes over successive market clearing intervals. Interestingly, satisfies the model , where is a positive definite and sparse matrix having non-positive off-diagonal entries, and is a sparse and low-rank matrix (Section 3). Recognizing such a rich structure is the first contribution of this paper. Inspired by recent advances in learning using sparse and low-rank models, a novel recovery scheme is devised, and constitutes the second contribution. As a third contribution, an efficient algorithm based on the alternating direction method of multipliers (ADMM) is developed (Section 4). Numerical tests on market data generated on the IEEE 14-bus benchmark corroborate the validity of our findings (Section 5).
Notation. Lower- (upper-) case boldface letters denote column vectors (matrices); while denotes the all-ones (all-zeros) vector, and the identity matrix. Prime stands for matrix transposition; is the matrix inner product; is the entrywise (Hadamard) matrix product; is the diagonal matrix having scalars on its diagonal; and is the matrix determinant. Symbol () is reserved for the set of real symmetric (positive semidefinite) matrices. Regarding matrix norms, is the nuclear norm (sum of matrix singular values); is the Frobenius norm; and is the -norm.
2 Data Modeling and Problem Statement
Consider a power grid represented by the graph , where the set of nodes corresponds to buses, and the edges in to transmission lines. Grid connectivity is captured via the branch-bus incidence matrix : if its th row corresponds to edge , then , , and zero otherwise [15]. For a connected grid, the nullity of is one; and by construction .
To describe real-time market clearing, the linear DC power flow model is summarized next. Let be the voltage phase and the active power injected at bus . The active power flowing from bus to bus over line is approximated as , where is the line reactance. Collecting and in and , respectively, it holds that , where . Power conservation dictates that with defining the weighted Laplacian of . As such, is positive semidefinite, and for a connected grid, its zero eigenvalue has multiplicity one with being the corresponding eigenvector, i.e., ; see e.g., [15].
Real-time electricity prices are determined by solving the network-constrained economic dispatch problem [16]. In a simple, yet sufficiently representative form, the latter is typically formulated as the linear program
(1a) | ||||
s.to | (1b) | |||
(1c) | ||||
(1d) |
which is solved regularly to determine the (incremental) power schedules for the five-minute interval ahead.
If the -th bus bidder is a generator, is non-negative and bounded by (1b) in accordance to generation limits and selling offers. Moreover, in (1a) is the lowest price in dollars per MWh the generator is willing to be paid to inject power in bus . If the -th bus bidder is a consumer, is non-positive and it is either fixed to the value predicted for a fixed load, or, it is bounded by the buying offer for elastic loads. When the load is elastic, the corresponding is the highest price the consumer is willing to pay for withdrawing electricity at bus ; or zero for fixed loads. Zero-injection buses are modeled by simply setting in (1b). Having multiple bidders and/or multi-block bids per bus does not harm the generality of the ensuing results.
The outcome of (1) depends not only on participant bids , but also on the underlying physical system limitations. Indeed, bus injections are related to nodal voltage phases via (1c), while power flows cannot exceed the transmission capacity limits imposed on both flow directions via (1d).
The optimization in (1) can be simplified after recognizing that due to the nullspace of , the vector pairs are minimizers of (1) too. This phase shift ambiguity is resolved by fixing , in which case bus 1 serves the role of a reference bus. Consider the partition and ; and let be the full-column rank matrix obtained after maintaining all but the first column of . Define also the reduced grid Laplacian , which is strictly positive definite.
The set described by (1c), (1d), and , can be represented by , , and . Upon eliminating , the problem in (1) can be rewritten as
(2a) | ||||
s.to | (2b) | |||
(2c) | ||||
(2d) |
To describe the pricing mechanism, let denote the optimal Lagrange multiplier associated with constraint (2c); and and be the optimal Lagrange multipliers corresponding to the lower and upper limits of (2d), respectively. The vector of ex-ante real-time LMPs is
(3) |
where , and is a relatively small term accounting for the heat losses on transmission lines. Oftentimes, is approximated as the product of times the gradient of the system-wide loss evaluated at .
Evidenced by (3), LMPs consist of three components:
(C1) the marginal energy component ;
(C2) the marginal congestion component; and,
(C3) the marginal loss component .
Given that price components are usually announced separately, our focus will be on the topology-related information that can be extracted from (C2) and (C3). Even if LMPs are announced as a sum, (C1) can be easily subtracted from since , and (C3) can be modeled as low-variance noise.
As mentioned earlier, the problem in (1) or (2) is solved every five minutes. The grid topology captured by , , and , will be assumed approximately constant over a period of such intervals. However, the triplets are time-varying: generator offers can be altered hourly and load demands fluctuate on a 5-min basis. Let vector denote the LMPs derived for a specific 5-min interval indexed by . With a slight abuse of notation, comprises only the (C2) plus (C3) LMP components assuming that (C1) has been removed. Due to (3) and for , it holds that
(4) |
where and captures unmodeled physical system variations and the loss component. Upon collecting prices and Lagrange multiplier differences over time periods in and , model (4) reads
(5) |
where and . Given in (5), our goal is to identify the topology (Laplacian) matrix .
3 Market Data Factorization
Finding from (5) given only the observed prices in constitutes a blind recovery problem, yet with a very rich structure as recognized next. According to (3), the -th column of is . Due to complementary slackness in (2), the -th entry of () is strictly positive only when line is congested, i.e., it has reached its lower (upper) capacity limit. Since typically only a few transmission lines are congested over a period of intervals, is expected to be sparse.
The sparsity of endows with two important properties. Note that the -th column of can be expressed as . Hence, the -th entry of is non-zero only if at least one of the lines adjacent to node is congested at time . Assuming that the entry is nonzero only for a few across all , vectors are expected to be sparse too. Moreover, the columns of are expressed by linearly combining only a few of . Thus, is not only sparse, but low-rank as well.
By definition of the original graph Laplacian , its -th off-diagonal entry is if , and 0 otherwise. Granted that power grids are sparingly connected, most of the off-diagonal entries of are zero, while the rest have negative values. These properties readily carry over to .
A meaningful matrix factorization exploiting the aforementioned properties can be found by solving the problem
(P1) | ||||
s.to |
where are positive tuning parameters; and matrix is introduced to select the off-diagonal entries of . The cost in ((P1)) involves a least-squares data fitting term and four regularization components: and promote sparse and minimizers; encourages a low-rank solution for ; while ensures and excludes the uninteresting case of having zero minimizers.
Problem ((P1)) can be interpreted as a regularized maximum likelihood estimator (MLE) of . Suppose the noise term in (4) is drawn from a multivariate Gaussian distribution with zero mean and covariance matrix for some . The negative log-likelihood of reads . Assuming independence across , the MLE of given the observations is found as the solution of
(6) |
The extra penalties and the constraint in ((P1)) regularize the MLE to further promote the structural properties of .
4 Algorithms
Optimization problem ((P1)) is convex and it can be actually expressed as a semidefinite program (SDP). However, high-dimensional market data ( and in the order of thousands), exclude the possibility of using standard interior point-based SDP solvers. Instead, an efficient algorithm based on the alternating direction method of multipliers (ADMM) is developed next; see e.g., [17] for a review on ADMM. First, ((P1)) is equivalently reformulated as
(P2) | ||||
s.to | (P2a) | |||
(P2b) |
To efficiently handle the non-differentiable cost and the constraints in ((P1)), the original variables and are replaced by and in ((P2)). Consensus among these variable duplicates is guaranteed by the constraints in ((P2b)).
If , , and are the Lagrange multipliers associated respectively with the three constraints in ((P2b)), the augmented Lagrangian function of ((P2)) is given by
(7) | ||||
where is a predefined constant and denotes the indicator function for set . Each ADMM iteration consists of two primal and one dual update steps. During a primal update step, is minimized over a subset of variables, while the remaining variables are fixed to their most recent values in a Gauss-Seidel fashion. In the dual update step, Lagrange multipliers are updated via gradient ascent. Letting denote the iteration index, our ADMM algorithm cycles through the three steps detailed next.
In the first step, are updated as the minimizers of , where the superscript indicates the variable value at the end of the -th iteration. The minimization is separable over and . Regarding the update of , upon completing the squares and ignoring constant terms, turns out to be the solution of
(8) |
that is provided in closed form as
Optimization with respect to entails solving
(9) |
whose minimizer becomes available in closed form as [17]
(10) |
and is the so termed soft thresholding operator applied entrywise to matrix .
In the second step, are updated as the minimizers of . The optimization is again separable over the three variables. Specifically, the update is found as the solution of
(11) |
Problem (11) decouples over the entries of and after using the constraint to simplify the cost, it can be shown that
Next, variable is updated as the minimizer of
(12) |
By Lemma 1 (proved in the Appendix), the minimizer of (12) is neatly expressed as
(13) |
where with the eigenvalue decomposition . Finally, upon completing the squares, is found via
(14) |
whose solution can be expressed as [18, Th. 2.1]
(15) |
where , with the singular value decomposition . The closed-form update in (15) is essentially a soft thresholding operator on the singular values of its matrix argument.
In the third step, the Lagrange multipliers are updated via gradient ascent simply as
(16) | ||||
(17) | ||||
(18) |
Bus number | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Type (Generator/Load/Zero) | G(R) | G | G | L | L | G | Z | G | L | L | L | L | L | L |
Unmodulated mean loads [MW] | - | 21.7 | 94.2 | 47.8 | 7.6 | 11.2 | - | - | 29.5 | 9.0 | 3.5 | 6.1 | 13.5 | 14.9 |
Mean generation bid [$/MWh] | 18 | 31 | 30 | - | - | 15 | - | 22 | - | - | - | - | - | - |
Upper generation bounds [MW] | 200 | 140 | 100 | - | - | 100 | - | 100 | - | - | - | - | - | - |
5 Numerical Tests
The novel topology recovery approach was evaluated using market data generated for the IEEE 14-bus grid. Ex-ante real-time prices were simulated over a one-day period by solving (2) for 5-min intervals using YALMIP and SDPT3 solvers [19], [20]. Matrices , , and were obtained from MATPOWER [21]. Among the 14 buses, buses are generators (bus 1 is the reference), bus 7 is a zero-injection bus, and the rest are inelastic loads. Lower generation bounds were set to zero and upper generation bounds were selected as listed in Table 1. Following MATPOWER’s line ordering, the flow limits on the transmission lines were set to 70, 90, 50, 70, 50, 20, 50, 70, 90, 90, 20, 70, 50, 70, 20, 50, 90, 50, 50, 70MW.
Generation bids were modeled as uniform random variables having the means shown in Table 1, and standard deviation 2.88$/MWh. Loads were drawn as independent Gaussian random variables with standard deviation MW, and means simulated as the product of the loads shown in Table 1 times a factor of 0.77, 0.74, 0.73, 0.74, 0.77, 0.83, 0.90, 0.95, 0.97, 0.99, 0.99, 1.00, 0.99, 0.99, 0.97, 0.96, 0.94, 0.92, 0.91, 0.90, 0.91, 0.88, 0.81, 0.75 depending on the hour of the day. These factors were the normalized hourly loads over the Midwest Independent System Operator market in June 1st, 2012. After solving (2) for all 288 intervals, the price matrix was constructed as according to (4)-(5). Among the 288 intervals, congestion occured only in intervals; and these non-zero ’s comprised the columns of the nominal price matrix .
Before inferring from ((P1)), the parameter vector needs to be tuned. Since cross-validation cannot be applied here [22], the following heuristic was used instead. A grid of values was chosen for all ’s. For each candidate , 10% of the entries of chosen uniformly at random were considered unknown, and set to zero to yield the observed price matrix . Matrices were then found as the minimizers of ((P1)) using . The entries of assumed unknown were recovered as the corresponding entries of . The process was repeated 10 times for each . The squared reconstruction error on the missing entries of was averaged over all 10 runs. The configuration attaining the lowest aggregate error was .
Problem ((P1)) was solved using the algorithm of Section 4. Parameter was set to , while the stopping criteria suggested in [17, Sec. 3.3] were employed. The ADMM solver terminated in less than 1 min on a 2.4 GHz Intel Core i7 processor with 4 GB RAM, whereas the SDPT3 solver could not run. Topology recovery results are shown in Figs. 1-2. The results are quite encouraging given that the novel scheme is based solely on the observed prices . Collecting prices over a longer observation period and over more diverse market conditions are expected to offer enhanced recovery.
6 Conclusions
Recovering grid topology from LMPs was the theme of this work. It was first recognized that the price matrix admits an interesting bilinear decomposition. A convex optimization problem involving sparse and low-rank regularizers was formulated to reveal the constituent matrix factors. The problem was solved by an iterative ADMM-based algorithm entailing only closed-form updates. The novel scheme yielded encouraging topology recovery results on market data generated using the IEEE 14-bus grid. Performing identifiability analysis and experimenting with real market are challenging future directions.
Lemma 1.
Given and , the convex problem where admits the minimizer
(19) |
with the eigen-decomposition .
Consider the decomposition into the symmetric and the anti-symmetric . Recall that the inner product between a symmetric and an anti-symmetric matrix is zero. For to be the minimizer, it suffices to show that for all . The gradient of at is . Since for all , what remains to be shown is for all . Observing that , , and share the same eigenvectors , it is easy to show that , which completes the proof.
Footnotes
- footnotetext: Work in this paper was supported by the Inst. of Renewable Energy and the Environment (IREE) under grant no. RL-0010-13, Univ. of Minnesota, and NSF Grant ECCS-1202135.
References
- A. Abur and A. Gómez-Expósito, Power System State Estimation: Theory and Implementation. New York, NY: Marcel Dekker, 2004.
- Y.-Y. Lee, J. Hur, R. Baldick, and S. Pineda, “New indices of market power in transmission-constrained electricity markets,” IEEE Trans. Power Syst., vol. 26, no. 2, pp. 681–689, May 2011.
- V. Kekatos, Y. Zhang, and G. B. Giannakis, “Electricity market forecasting via low-rank multi-kernel learning,” IEEE J. Sel. Topics Signal Process., 2014 (accepted). [Online]. Available: http://arxiv.org/abs/1310.0865
- Y. Liu, P. Ning, and M. K. Reiter, “False data injection attacks against state estimation in electric power grids,” ACM Trans. Info. and System Security, vol. 14, no. 1, pp. 13:1–13:33, May 2011.
- O. Kosut, L. Jia, J. Thomas, and L. Tong, “Malicious data attacks on the smart grid,” vol. 2, no. 4, pp. 645–658, Dec. 2011.
- L. Jia, J. Kim, R. J. Thomas, and L. Tong, “Impact of data quality on real-time locational marginal price,” IEEE Trans. Power Syst., 2014, early access.
- L. Xie, Y. Mo, and B. Sinopoli, “Integrity data attacks in power market operations,” vol. 2, no. 4, pp. 659–666, Dec. 2011.
- R. Emami and A. Abur, “Tracking changes in the external network model,” in Proc. North American Power Symposium, Arlington, TX, Sep. 2010, pp. 1–6.
- J. E. Tate and T. J. Overbye, “Double line outage detection using phasor angle measurements,” in Proc. of IEEE Power & Energy Society General Meeting, Calgary, Alberta, Canada, Jul. 2009, pp. 1–5.
- H. Zhu and G. B. Giannakis, “Sparse overcomplete representations for efficient identification of power line outages,” IEEE Trans. Power Syst., vol. 27, no. 4, pp. 2215–2224, Nov. 2012.
- T. Erseghe, S. Tomasin, and A. Vigato, “Topology estimation for smart micro grids via powerline communications,” IEEE Trans. Signal Process., vol. 61, no. 13, pp. 3368–3377, Jul. 2013.
- X. Li, V. Poor, and A. Scaglione, “Blind topology identification for power systems,” in Proc. IEEE Smart Grid Communications Conf., Vancouver, BC, Canada, Oct. 2013.
- M. He and J. Zhang, “A dependency graph approach for fault detection and localization towards secure smart grid,” vol. 2, no. 2, pp. 342–351, Jul. 2011.
- S. Bolognani and L. Schenato, “Identification of power distribution network topology via voltage correlation analysis,” in Proc. IEEE Conf. on Decision and Control, Florence, Italy, Dec. 2013.
- G. B. Giannakis, V. Kekatos, N. Gatsis, S.-J. Kim, H. Zhu, and B. Wollenberg, “Monitoring and optimization for power grids: A signal processing perspective,” IEEE Signal Process. Mag., vol. 30, no. 5, pp. 107–128, Sep. 2013.
- D. Kirschen and G. Strbac, Power System Economics. West Sussex, England: Wiley, 2010.
- S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Found. Trends Mach Learning, vol. 3, pp. 1–122, 2010.
- J.-F. Cai, E. J. Candes, and Z. Shen, “A singular value thresholding algorithm for matrix completion,” SIAM J. on Optimization, vol. 20, no. 4, pp. 1956–1982, 2008.
- J. Lofberg, “A toolbox for modeling and optimization in MATLAB,” in Proc. of the CACSD Conf., 2004. [Online]. Available: http://users.isy.liu.se/johanl/yalmip/
- R. H. Tutuncu, K. C. Toh, and M. Todd, “Solving semidefinite-quadratic-linear programs using SDPT3,” Mathematical Programming Ser. B, vol. 95, pp. 189–217, 2003.
- R. D. Zimmerman, C. E. Murillo-Sanchez, and R. J. Thomas, “MATPOWER: steady-state operations, planning and analysis tools for power systems research and education,” IEEE Trans. Power Syst., vol. 26, no. 1, pp. 12–19, Feb. 2011.
- T. Hastie, R. Tibshirani, and J. Friedman, The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer Series in Statistics, 2009.