Online Energy Price Matrix Factorizationfor Power Grid Topology Tracking

Online Energy Price Matrix Factorization for Power Grid Topology Tracking


Grid security and open markets are two major smart grid goals. Transparency of market data facilitates a competitive and efficient energy environment, yet it may also reveal critical physical system information. Recovering the grid topology based solely on publicly available market data is explored here. Real-time energy prices are calculated as the Lagrange multipliers of network-constrained economic dispatch; that is, via a linear program (LP) typically solved every 5 minutes. Granted the grid Laplacian is a parameter of this LP, one could infer such a topology-revealing matrix upon observing successive LP dual outcomes. The matrix of spatio-temporal prices is first shown to factor as the product of the inverse Laplacian times a sparse matrix. Leveraging results from sparse matrix decompositions, topology recovery schemes with complementary strengths are subsequently formulated. Solvers scalable to high-dimensional and streaming market data are devised. Numerical validation using real load data on the IEEE 30-bus grid provide useful input for current and future market designs.



nline convex optimization; compressive sensing; alternating direction method of multipliers; economic dispatch; locational marginal prices; graph Laplacian.

1 Introduction

An independent system operator collects energy offers and bids, and dispatches power by maximizing the social welfare while meeting physical grid limitations. To guarantee competitive market operation, multiple data are communicated to market participants or are openly publicized, either in real-time or with certain delay [1]. Such market data may involve energy prices, bids and offers, congestion information, demand, and renewable generation. Looking forward, the smart grid vision calls for energy markets reaching the distribution level to promote participation, accounting for increased stochasticity at a finer time resolution [2]. New reliable market designs are hence to be developed.

From state estimation to load prediction, inference using data has been a major grid operation component. Facing smart grid challenges and opportunities, grid analytics are now extending to price and renewable forecasting, consumer preference learning, and cyber-physical attack detection [3]. Among grid learning tasks, topology monitoring is critical for security, market clearing, and billing. Although operators monitor grid topology via the generalized state estimator [4], topology tracking could be used for other purposes. Data attacks on the state estimator require precise physical network information [5]. Knowing congested transmission lines could assist in informed bidding or in market manipulation [6]. Line reactances could be used as a measure of electrical distance to cluster buses, or reveal influential nodes. Moreover, the Laplacian matrix of the graph representing a grid could capture the correlation across pricing nodes [7], or characterize the performance of decentralized algorithms.

Although there has been a long line of research regarding attacks on the state estimator, grid topology recovery using readily accessible market data has been overlooked. Stealth data attacks to the power system state estimator were first recognized in [8]. Their impact on state estimation and market outcomes was characterized in [9][10]. The possibility of data framing attacks deceiving the bad data processor were explored in [11]. Attacks and countermeasures on power system controllers have been studied in [12]. Procedures for detecting and identifying cyber-physical attacks have been also reported; see e.g., [13]. Designing cyber-physical attacks generally presumes the grid topology to be known [5], [14].

Detecting topological changes from the operator’s perspective has been studied in [15], [16]; while transmission line outages can be efficiently revealed via the sparse overcomplete representation of [17]. Grid topology recovery is cast as a blind factorization on the matrix of spatio-temporal nodal injections in [18]: Even though building on the sparsity and positive semidefiniteness of the grid Laplacian, [18] relies on linear independence across voltage phases. Considering a power line communication network, time delays of communication signals are leveraged to unveil the microgrid structure in [19]. By postulating a Gaussian Markov random field over nodal voltage phases, transmission network faults could also be localized [20]. Likewise for distribution grids, the topology recovery scheme of [21] exploits the sample covariance matrix of nodal voltage magnitudes.

All in all, existing topology recovery schemes presume access to a physical system quantity (power injections, voltage phases or magnitudes, communication delays) that is actually measured over all buses. In contrast, our previous work introduced the possibility of topology tracking using readily available cyber-system data [22]. The idea is that real-time energy prices are found by the system operator as the Lagrange multipliers of the network-constrained economic dispatch, which is a linear program (LP) typically solved every 5 minutes. Dispatch decisions are the primal variables of this LP, while grid topology and electricity offers/bids are its parameters. Observing the dual variables (prices) related to multiple offer/bids instances, the crux is to recover the quasi-stationary topology underlying this LP (Section 2). Our first contribution is recognizing that properties of the Laplacian matrix and sparsity in congested lines could be exploited: The matrix of spatio-temporal prices can be factorized as the product of a doubly positive matrix with sparse inverse times a sparse matrix (Section 3). Novel blind recovery schemes of complementary strengths constitute the second contribution (Section 4): Different from [22], the low-rank property of one of the matrix factors is not regularized here, thus significantly simplifying the problem. As our third contribution, algorithms handling big market data are developed based on the alternating direction method of multipliers and its online version (Section 5). Advancing tools from online optimization, an algorithm handling streaming market data is devised. Distinct from [22], such an online approach is pertinent to future smart grid market designs. Experiments with market data obtained using real load data over the IEEE 30-bus benchmark corroborate the validity of our findings (Section 6).

Notation. Lower- (upper-) case boldface letters denote column vectors (matrices); and denote the all-ones and all-zeros vectors. Symbols , , and , stand for matrix transposition, trace, and determinant, respectively. Symbol () is 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 .

2 Energy Price Data Model

Before delineating our price data model, this section reviews linear DC power flows and real-time energy markets.

2.1 Power Grid Modeling

Consider a power grid represented by the graph , where the set of nodes corresponds to buses, and the edges in to transmission lines. The grid topology is captured via the branch-bus incidence matrix  [3]. For a connected grid, the nullity of is one; and by definition, . If is the reactance of line and an diagonal matrix with , the bus reactance matrix can be defined as . Given that is the weighted Laplacian of , it is positive semidefinite, and is an eigenvector corresponding to ’s zero eigenvalue.

The DC power flow model can now be expressed in matrix-vector form. The active power flow from bus to bus over line can be approximated as , where is the voltage phase at bus ; while the power injection at bus is . By stacking and in and , respectively; it follows that and . By eliminating , the flows can be linearly expressed in terms of ; yet is non-invertible.

To resolve the singularity of , partition into the first and the rest of its columns as . For a connected , the reduced branch-bus incidence matrix has full column-rank. Thus, the reduced bus reactance matrix , is strictly positive definite. Setting , it readily follows


where .

2.2 Modeling of Real-Time Energy Markets

Building on this model, let us now review real-time energy markets. Energy markets determine the price for electricity by matching supply and demand. Due to time-varying demand and transmission grid limitations, the electricity cost varies across time and space (buses), giving rise to locational marginal prices (LMPs) [1]. Real-time markets are spot markets where hourly power schedules determined over the previous day are adjusted every five minutes to accommodate real-time deviations. Specifically, real-time LMPs are found via the network-constrained economic dispatch, typically formulated as the following LP [23]

(2a) (2b)

Problem (2) determines the incremental power injections for the upcoming 5-min interval indexed by . The optimum dispatch is found by minimizing the electricity cost in (2a); while satisfying the power limits in (2b), achieving the supply-demand balance via (2c), and confining line flows approximated by (1) to lie within a secure range [cf. (2d)]. The power injection bounds in (2b) model bid blocks.

By solving (2), the operator not only determines , but also calculates the LMPs as follows. Let be the optimal Lagrange multiplier associated with the supply-demand equality in (2c); and be the optimal Lagrange multipliers related to the lower and upper flow limits in (2d). By duality and upon defining , problem (2) can be equivalently expressed as


If is selected as the vector of nodal electricity prices at time and assuming are the actual marginal costs, then (3) reveals that maximizes the sum of the individual profits. In practice, to account for transmission line losses ignored by the DC model, LMPs are calculated as


where is a relatively small loss correction [23].

The LMPs in (4) consist of three summands: the marginal energy component (MEC) ; the marginal congestion component (MCC) ; and the marginal loss component (MLC) . According to (3), the MEC is the energy price at the reference bus (without loss of generality selected here as bus 1). When the power flow on line reaches the upper or lower limit at time , that is or , then line is deemed congested. Complementary slackness implies that if line is not congested at time , the -th entry of is zero. Apparently, if there are no congested lines and losses were ignored, all nodes would enjoy the same energy price .

2.3 Problem Statement

Depending on the market, the three LMP components are announced either separately or collectively as a sum. In the former case, the MCC is readily available. In the latter, the effect of MEC can be isolated by subtracting the first entry of from all entries of . It can be argued that subtracting the first entry does not harm the generality of this preprocessing step, even if the reference bus is not bus 1. Either way, collect all but the first bus prices in , for which we postulate:


where and captures the MLC. Slightly abusing terminology, will be henceforth termed the LMPs.

Market clearing occurs every five minutes, and market bids change partially over time, adapting to demand and generation fluctuations. Consider collecting the LMPs of (5) over the horizon of consecutive market intervals, and suppose grid topology remains invariant over . Upon stacking as columns of the matrices , , and , respectively, it follows from (5):


Model (6) asserts that if is ignored, the price matrix can be factorized as the product of the inverse grid Laplacian times matrix . With (6), topology recovery can be now formulated as the problem of finding given .

Remark 1.

Having multi-block offers and several bidders per bus does not harm generality of (4)-(6). Specifically, electricity offers and bids oftentimes consist of multiple blocks: For example, a generator may offer to produce the first 20MWh for at least 20MWh and the next 5MWh for at least 23MWh at the same bus. In this case, the corresponding should be decomposed as the sum of two extra optimization variables as with and ; while the summand in (2) is replaced by . Having multiple generators and/or consumers at the same bus is handled similarly. Either way, constraints (2c)-(2d) apparently remain unaltered. Hence, even though simplifying, (2) is sufficiently representative. Actually, Section 6 involves tests with multi-block offers.

3 Topology Recovery Approaches

If the MCCs are announced separately, matrix satisfies the noiseless counterpart of (6), namely


Decomposing into constitutes a blind matrix factorization problem. To uniquely recover the factors, their rich structure delineated next should be properly exploited.

Recall that is positive definite. Once has been recovered, the original grid Laplacian can be trivially found in light of the property . Note further that the -th entry of equals , if there is a line between buses and ; and zero otherwise. Granted power grids are sparingly connected, is not only sparse, but its off-diagonal entries are non-positive. Having positive eigenvalues and non-positive off-diagonal entries implies is an invertible M-matrix [24, Sec. 2.5]. Hence, has positive entries, i.e., .

As far as is concerned, its columns can be expressed as . Since many of the in (8) are expected to be zero [cf. Prop. LABEL:pro:sparsemu], can be also written as


where is the subset of congested lines at time . In other words, is a linear combination of few ’s. Given that are sparse, matrix is expected to be sparse too. Typically, only a few transmission lines become congested over a short market period (say one day): In the California ISO (CAISO) for example, only two transmission lines are typically congested [25]. Hence, it can be assumed that the overlap significantly, or that the locations of the non-zero entries of remain relatively time-invariant. Since ’s are linear combinations of those few ’s corresponding to congested lines, is expected to exhibit low rank. The invertibility of implies should have low rank too.

It will be useful also to recognize that the factorization in (7) is scale-invariant: If satisfies (7), so does the pair for all . To waive this inherent ambiguity, the maximum diagonal entry of is assumed to be unity. Due to this normalization, matrix should satisfy and .

Leveraging these properties, one could recover by solving the optimization problem:

(9a) (9b)

where is the -(pseudo)norm of a matrix counting its non-zero entries, and is a parameter balancing the sparsity between the two matrices. Problem (9) finds the sparsest pair that satisfies model (7) and the structural constraints of . Different from [22], the rank of is not penalized here, since and the invertibility of enforce anyway.

Minimizing the -norm is in general NP-hard [26]. Following advances in compressive sensing [27], the -norm will be surrogated by the -norm to yield the convex problem


where , , and . Two observations are in order regarding (10).

First, since that the diagonal entries of are strictly positive, in (9a) has been replaced by the sum of the off-diagonal entries of in their absolute values, that is

where the first equality comes from the non-positive off-diagonal entries of , and the rest from properties of the trace.

Second, ideally should be enforced to be strictly positive definite, i.e., . Nonetheless, current optimization algorithms cannot guarantee the minimizer to lie in the interior of the feasible set. On the other hand, imposing admits the trivial solution . As a remedy, the term has been added in the cost of (10) to confine in the interior of the positive semidefinite cone .

By eliminating , (10) can be equivalently transformed to


The strict convexity of guarantees that (11) and hence (10) have unique minimizers.

When comprise of both MEC and MLC, model (6) is more pertinent than the exact model in (7). The non-convex problem in (9) could be then replaced by


for . The approach in (12) aims at minimizing the least-squares distance between and , while seeking sparse and a low-rank . However, minimizing the -norm and the matrix rank constitutes an NP-hard problem. In the same spirit (9) was surrogated by (10), the hard problem in (12) is approximated by the convex problem


where and serves as a convex surrogate for . A solver and recovery results from (13) can be found in [22]. Given that MCCs are typically announced independently, our focus will be henceforth on model (6).

4 Batch Topology Recovery Scheme

Although (10)-(11) could be solved by commercial software for relatively small problems, interior point-based solvers cannot handle and larger than a few hundreds. There are two main optimization challenges here: the objective term and the feasible set . Regarding the former, not only it is non-differentiable, but also involves a linear transformation of . Note that is the intersection of the positive definite cone and a shifted version of the positive cone. Albeit projection over each of these cones is relatively easy, there is no closed-form solution for projecting on .

Given these challenges, an algorithm based on the alternating direction method of multipliers (ADMM) is derived next. ADMM targets solving problems of the form [28]


where and are convex functions; and are convex sets; and model the linear equality constraints coupling variables and . In its normalized form, ADMM assigns a Lagrange multiplier for the equality constraint and solves (14) by iterating over the recursions


for some . To apply ADMM and end up in efficient updates for (10), we first replace variable with three copies , , and , to yield the equivalent problem

(16a) (16b)

Let , , and , denote the Lagrange multipliers corresponding to (16b), (16c), and (16d), respectively. Partitioning variables into and , ADMM iterates through the next three steps.

At the first step of iterate , the variable is updated given and by solving


The solution of (17) is provided in closed form as .

During the second step, ADMM updates the second block of variables given and . Yet the optimization decouples over the three variables. Specifically, variable is updated as the solution of


whose minimizer is , where the minimum operator is understood entry-wise.

Variable is updated as the minimizer of


which can be readily found as follows [22, Lemma 1]: Define the operator for some as


where is the eigenvalue decomposition of the symmetric matrix . Then, the solution to (19) is


Variable is updated by solving


Problem (22) is separable across the entries of , admitting the closed-form minimizer:


where is the soft thresholding operator defined as

applied entry-wise in (23).

In the third step, the Lagrange multipliers are updated as

0:  Price matrix and .
1:  Initialize and .
2:  Initialize and .
3:  for  do
4:     Update from (17).
5:     Update from (18).
6:     Update from (21).
7:     Update from (23).
8:     Update multipliers from (24).
9:  end for
Algorithm 1 Batch Topology Recovery Scheme

The algorithm is summarized as Alg. 1, and its convergence is guaranteed for all  [28].

5 Grid Topology Tracking

The topology recovery scheme of Section 4 presumes that:
    (c1) the power system topology remains unchanged, and
    (c2) prices are available for the entire period .
In reality, future power grids may be reconfigured frequently for dispatching and maintenance, while real-time LMPs are expected to be announced at a fast rate over thousands of buses; hence, rendering conditions (c1)-(c2) unrealistic.

To cope with these challenges, we first modify the recovery scheme of (10) to address (c1). Specifically, rather than enforcing the constraint , our idea here is to look for sparse yielding a small least-squares error by solving

for . Upon eliminating , the last minimization can be shown to be equivalent to


where , and is the so termed Huber function


Notice that in (11), the entries of are arguments of the absolute value cost. In contrast, the objective in (25) penalizes small entries of with a quadratic cost, and large ones with the absolute value cost.

To cope with (c2), solvers for streaming rather than batch pricing data are developed next. The desiderata here is online schemes where topology estimates are updated every time a price vector is publicized. Advances from online optimization are particularly suitable for this task [29]. Tailored to big data processing, many online convex optimization algorithms aim at solving problems of the form


where depends on the -th datum, and is a regularizer, i.e., a function leveraging prior information on .

Tailoring our grid topology recovery task to the online optimization setup, consider the general problem


By selecting , problem (28) yields (11); whereas, for , (28) is equivalent to (25). Apparently, is price-dependent, and the other two terms in the objective of (28) can be thought of as regularizers for . To solve (28), we resorted to the online ADMM algorithm of [28] that cycles through:


Comparing (29) with its batch counterpart in (15), the iteration index in (29) coincides with the data index , while the first step in (29a) entails only the current together with the proximal term for some . Regarding its convergence, the algorithm attains sublinear regret in terms of both the cost and the constraint violation [28, Th. 4].

Building on (16), introduce copies of to express (28) as

(30a) (30b)

Similarly to (16), let and be the Lagrange multipliers corresponding to constraints (30b) and (30c), respectively.

As soon as the -th price vector is announced, a cycle of the online ADMM of (29) is initiated. In its first step, is updated via (29a), which upon completing the squares yields


where . Problem (31) could be reformulated and solved as a linearly-constrained quadratic program. Interestingly, the minimizer of (31) can be found in closed form for both choices of . Specifically, if , the next claim is shown in the Appendix:

Proposition 1.

Given , the minimizer


is given by , where , and the operator is defined as applied entry-wise.

By Prop. 1, can be found as a rank-one update of


where . The key observation here is that having a single -norm rather than [cf. (11)] enabled the simple update of (33).

For the Huber cost, the next claim is shown in the Appendix:

Proposition 2.

Given , the minimizer


for is given by , where , and is defined as


applied entry-wise.

Based on Proposition 2, when , the minimizer of (31) is


During the second step of iteration , matrices are updated similarly to (18)-(21) as


At the third step, the Lagrange multipliers are updated as


To summarize, grid topology recovery using inexact streaming pricing data can be performed using Algorithm 2.

0:  .
1:  Initialize .
2:  Initialize .
3:  for  do
4:     Acquire price vector .
5:     Update using (33) or (36).
6:     Update from (37) and (38), respectively.
7:     Update from (39) and (40), respectively.
8:  end for
Algorithm 2 Topology Recovery Tracking Scheme

6 Experimental Validation

Figure 1: Topology of the IEEE 30-bus system [30].

The novel topology recovery approaches are tested next using real load data on the IEEE 30-bus benchmark [30]. The latter comprises 18 load buses, 6 generators, and 6 zero-injection buses. The transmission network consists of 41 lines with ratings ranging from 16 to 130 MVA as listed in [30].

Generator Block offers [MWh,$/MWh]
1 (30,26) (20,36) (20,44) (10,50)
2 (20,21) (20,28) (20,35) (20,43)
13 (15,38) (15,42) (10,47)
22 (10,16) (10,27) (10,41) (10,54) (10,66)
23 (15,34) (15,40)
27 (30,35) (15,39)
Table 1: Generation Offers
   0.001 0.01 0.1 1 10
0.001 0.9 1.2 1.5 2.5 5.9
0.01 0.9 1.0 1.4 2.5 5.9
0.1 0.8 1.0 1.6 2.4 5.9
1 0.4 2.8 1.9 2.7 5.9
10 6.5 5.6 5.9 6.0 6.0
Table 2: Average bus degree attained for

Regarding offers, the benchmark provides generation capacities and quadratic generation costs [30]. To comply with market practices, the generation costs were first approximated by convex piece-wise linear functions yielding the block offers of Table 1. The original costs were scaled up by 10 to reflect current wholesale electricity cost levels. To model small fluctuations in costs, the nominal offers of Table 1 were shifted by a deviation uniformly distributed in  $/MWh.

For consumption, apart from the 18 load buses, generator buses 2 and 23 have load demands too, resulting in a total of 20 loads. The IEEE 30-bus benchmark provides a single realization of load demands. To simulate multiple realistic demands, we used the actual load data publicized for the Global Energy Forecasting (GEF) competition 2012 [31]. These data are the hourly energy consumptions over 20 sites. To match the load levels of the IEEE 30-bus grid, all loads were scaled down by a factor of 7. The 20 demand sequences from December 23, 2007, were assigned to buses so that the average consumption per bus matched the demand specified by the benchmark. Hourly loads were perturbed by a zero-mean Gaussian variation having standard deviation 10 times smaller than the nominal value to account for 5-min load fluctuations.

Real-time prices were generated by solving (2) for one day, i.e., 288 5-min intervals, and MCCs were announced separately. Lacking any day-ahead market information, the system was assumed to be dispatched entirely through the real-time market. Out of the 288 dispatches, 3 were infeasible and 45 experienced no congestion (occurred primarily over nighttime). Our experimental validation utilized the remaining MCC price vectors. It is worth stressing that only lines (1,2), (15,23), and (6,28) became congested.

(a) Actual Laplacian matrix.
(b) Laplacian matrix found by Alg. 1 for .
Figure 2: Laplacian matrix for the IEEE 30-bus system.

Upon collecting prices over an entire day, Alg. 1 was run on a 2.4 GHz Intel Core i7 (4GB RAM) laptop computer using MATLAB. Before running the algorithm, parameters were selected. Although such parameters are typically tuned using cross-validation, this methodology becomes cumbersome for our problem. Assuming the average node degree for the grid of interest to be known, were tuned so that the estimate had the same average degree. Given the scale ambiguity, the algorithm outcome was normalized by its maximum diagonal entry, and entries with absolute value smaller than were set to zero.

Algorithm 1 was run for taking the values . Regarding , the convergence rate for the objective (constraint violation) is proportional (inversely proportional) to  [28]. For the problem at hand, setting was empirically observed to provide a good trade-off. Since the average degree of the IEEE 30-bus grid is 2.68, the estimated node degrees obtained in Table 2 hint that could be both set to 1. The actual and the recovered Laplacian matrix for the IEEE 30-bus benchmark are shown in Fig. 2.

Figure 3: Tracking lines using streaming pricing data for January 2008.

To evaluate the online scheme, tests on real-time prices collected over January 2008 were conducted. Consumption data were generated by scaling GEF competition loads so that the maximum daily per-site value was 1.6 times the benchmark demands [31]. Tracking ability was tested by simulating a grid reconfiguration on January 15: lines (2,6) and (23,24) were exchanged for lines (2,7) and (23,26), respectively. Among the 8,928 intervals, infeasible dispatches and dispatches without congestion were ignored yielding 7,220 effective clearings. Alg. 2 with the update of (36) was initialized to the batch solution obtained from Alg. 1. Parameters and were set to yielding sublinear regret [28], while was set to 1. Figure 3 depicts the tracking behavior of Alg. 2. The estimated normalized reactance for line (10,17), i.e., entry , remained relatively invariant. Line (2,7) was initially erroneously detected as active, yet it was adjusted after Jan. 15, while reactance (2,6) approached zero. Interestingly, the replacement of line (23,24) by (23,26) was promptly detected.

7 Conclusions

Grid topology recovery using publicly available energy prices was the subject of this work. Upon exploiting the way real-time LMPs are obtained, recovery approaches with complementary strengths were developed. Advances in compressive sampling and online convex optimization proved to be useful for grid topology tracking. Experimental validation using real consumption data on a benchmark grid corroborated the risk of unveiling the power network structure. Numerical tests using a month-long price dataset showed the possibility of tracking grid reconfigurations too. The recovery performance could be enhanced further in envisioned smart grids: In competitive markets, rapidly changing offers and bids could probe the dispatch linear program in a richer way, while market data announced at higher rates could provide even more information. Leveraging heterogeneous market data and characterizing identifiability are interesting research directions.


[Proof of Proposition 1] Strict convexity of implies that (32) admits a unique minimizer. First-order optimality conditions assert that belongs to the subdifferential of evaluated at . By definition, the subdifferential of at is , where


is the -th entry of , and denotes the -th row of