IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 4, NOVEMBER 2011

Efficient Estimation of Critical Load Levels Using Variable Substitution Method Rui Bo, Member, IEEE, and Fangxing Li, Senior Member, IEEE

Abstract—An accurate estimation of critical load levels (CLLs) is beneficial in evaluating the confidence level for the predicted locational marginal prices (LMPs) in power markets. This paper proposes an efficient algorithm to quickly identify CLLs under the with-loss-DCOPF framework. The proposed method approximates the trajectory of DCOPF solutions with respect to the system load variations. The quadratic coefficients of the trajectory can be efficiently obtained by solving a small-scale nonlinear system of equations formed by the variable substitution method, which does not require additional OPF runs and is thus more computationally efficient than existing methods such as binary search and interpolation. The obtained trajectory can be used for a quick calculation of OPF solutions without re-running OPF. Further, through trajectory extrapolation, new congestions and other binding constraints can be identified simultaneously with the CLL predictions. Particularly, the proposed method is not restricted to a specific DCOPFbased model and can be easily applied to other DCOPF-based variants. Numeric studies on a revised IEEE 30-bus system are presented to illustrate the high computational efficiency and accuracy of the approach. It is further tested on the IEEE 118-bus system using four different DCOPF-based models to illustrate its model adaptability. Index Terms—Congestion prediction, critical load level (CLL), DCOPF, energy markets, locational marginal pricing (LMP), marginal units, optimal power flow (OPF), power markets, power system planning, variable substitution.

I. INTRODUCTION OCATIONAL marginal pricing (LMP) has been the dominant pricing mechanism in deregulated power systems. It plays a key role in market operation and settlement. LMP-based study has also been increasingly integrated into power system planning such as value-based transmission planning for transmission grid operators [1], where LMPs and economic benefits are typically predicted through exhaustive hourly simulation under a projected timeframe and scenarios. In deregulated markets, market participants perform price prediction on a regular basis to optimize their participation in order to maximize profits. These studies and activities extensively rely on efficient

L

Manuscript received November 22, 2010; revised March 11, 2011; accepted April 18, 2011. Date of publication May 27, 2011; date of current version October 21, 2011. This work was supported by the National Science Foundation under grant ECCS 1001999. Disclaimer: The views expressed in this paper are solely of the authors and do not necessarily represent those of Midwest ISO. Paper no. TPWRS-00943-2010. R. Bo is with Midwest Independent Transmission System Operator, St. Paul, MN 55108 USA (e-mail: [email protected]). F. Li is with the Department of Electrical Engineering and Computer Science, The University of Tennessee, Knoxville, TN 37996 USA. Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TPWRS.2011.2149548

Fig. 1. NYISO “CAPITAL” Zone LMP w.r.t. load based on real-time market data for 10–11 AM on July 24, 2009.

and reliable computation of LMPs. However, LMP tends to be very unstable in certain special scenarios. For instance, it has been reported that LMP may demonstrate a step change phenomenon at certain load levels, which are referred to as critical load levels (CLLs) [2], [3]. The step change pattern has been identified from real operation data as illustrated in Fig. 1, which shows an LMP versus load levels graph based on the 5-min data retrieved from the NYISO website. The data points represent the actual LMP of the CAPITAL zone in NYISO from 10–11 AM on July 24, 2009. The data suggests a step-change pattern of LMP with respect to load change around certain load levels, as illustrated by the dotted curve. The LMP sensitivity with respect to load can be effectively calculated using the perturbation method [4], [5], and therefore, LMP at nearby load levels may be estimated. However, in the presence of CLLs, one should be cautious in interpreting the LMP sensitivity as it is subject to step change when the perturbation is very close to a CLL. Therefore, it is crucial to locate the CLLs along with the prediction of LMPs. In [6], the importance of CLL is further demonstrated based on load uncertainty and the probabilistic LMP model, which reflects the stochastic feature of LMPs and has a significant potential impact in decision making. Methods that can be used for estimating CLLs have been reported in recent literature [3], [7], [8]. In [3], CLLs are calculated through parametric analysis for the lossless DCOPF model. The proposed simplex-like algorithm is very computationally effective. It has been reported that the algorithm can achieve performance that is many times faster than the binary

0885-8950/$26.00 © 2011 IEEE

BO AND LI: EFFICIENT ESTIMATION OF CRITICAL LOAD LEVELS USING VARIABLE SUBSTITUTION METHOD

search method in locating CLLs given the same precision criterion. In [7], local sensitivities are extended to a broader range for predicting system statuses for the lossless DCOPF model. These approaches, however, cannot address the nonlinearity brought by power losses. For the ACOPF model, the CLLs can be well predicted by an interpolation method proposed in [8]. System statuses such as generation dispatch, line flows, and Lagrange multipliers associated with binding constraints have been numerically verified to follow a quadratic pattern with respect to load variations within two consecutive CLLs. The system statuses are extrapolated and then utilized to estimate the previous and next CLLs around the given operating point. This method generally requires two additional ACOPF runs at two load levels other than the given operating point. In power system planning and market operation, the withloss-DCOPF model has gained prevalent use in practice due to its retention of the simplicity of DCOPF and incorporation of power losses. In comparison with the ACOPF model, the simplicity in modeling of the with-loss-DCOPF framework allows more efficient algorithms with a reasonable accuracy achieved. In this regard, this paper explores the with-loss-DCOPF model and proposes an efficient algorithm for prediction of CLLs. The proposed algorithm does not need any additional OPF runs and can be solved very efficiently while achieving a high accuracy. This paper is organized as follows. First, the conventional with-loss-DCOPF model is briefly reviewed in Section II. Then, the binding constraints which can be employed to identify the trajectory of the optimal solution w.r.t. load variation are analyzed in Section III. In Section IV, the trajectory is assumed to follow quadratic polynomials, and a variable substitution method is proposed to solve for the quadratic coefficients which subsequently facilitate the estimation of CLLs, congestions, and new binding constraints. In Section V, the proposed method is illustrated in detail on a revised IEEE 30-bus system for the conventional with-loss-DCOPF model. In Section VI, the approach is further tested on the IEEE 118-bus system for the same with-loss-DCOPF model as well as three other variants of the model including a more recent fictitious nodal demand DCOPF model. The model adaptability, easy implementation and flexibility, and time complexity analysis of the proposed method are discussed in Section VII. Section VIII concludes the paper. II. REVIEW OF WITH-LOSS-DCOPF MODEL The with-loss-DCOPF model is represented by (1)–(4). To facilitate the study of CLLs, the model is written w.r.t. marginal unit and non-marginal unit output. Marginal units are referred to as the generators which will adjust their generation output to meet any sufficiently small load variation. In contrast, nonmarginal units are those having their output bounded at limits such as capacity rating or ramping capability (1)

2473

(2)

(3) (4) where is the generation of th marginal unit; is the generation of th non-marginal unit; is the per MWh cost of generator ; is the load at bus ; is the delivery factor of bus and is the delivery factor of generator ; is the net injection at bus and is the total injection at bus

; and ;

are mapping functions between the unit and to the bus index in , respec-

index in tively; is the total system loss; is the generation shift factor of line with respect to bus ; is the transmission limit of line ; are the upper and lower generation limits of the th marginal unit, respectively; are the sets of marginal and non-marginal unit(s), respectively; is the set of load buses; and is the set of transmission line(s). A few variants of with-loss-DCOPF models exist [9], [10]. The one described in (1)–(4) is studied in this paper as it aligns with the present practices at a number of ISOs/RTOs [12]–[15]. Nevertheless, the selection of the model representation is only for illustration purposes. The method proposed in this paper is not restricted to any specific modeling approach, and can be potentially applied to other nonlinear DCOPF-based models, as will be demonstrated in Section VI. The energy balance equation, as denoted by (2), is formulated using delivery factors. The delivery factors are defined on a marginal basis and therefore can cause doubled power losses. It has been proven that an offset item needs to be added into the energy balance equation [11]. The delivery factors also contribute to the nonlinearity of the with-loss-DCOPF model due to its interdependency on the generation dispatch. Generation shift factors (GSFs), also known as power transfer distribution factors (PTDFs), are used to formulate the power flow equations in (3). Equations (1)–(4) may differ from other documented with-loss-DCOPF models in that they are written in terms of marginal and non-marginal unit generation outputs, instead of the general generation output. Furthermore, it is a more generic model because it does not require a reliable estimation of power losses from real time data or historical data and therefore can

2474

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 4, NOVEMBER 2011

be used in planning models. It should also be stressed that staircase generation bids, namely, piece-wise linear bid curves, are utilized in the model since it is a common practice in contemporary power markets. This model is a nonlinear programming model and thus the simplex-like method [3] does not fit. The quadratic interpolation method for the ACOPF framework [8] may apply to this withloss-DCOPF model; however, it should be noted that this model is far simpler than its ACOPF counterpart due to the absence of reactive power and bus voltage magnitude. Furthermore, the interpolation method generally requires two additional OPF runs, which will be shown unnecessary in this paper in the with-lossDCOPF modeling context.

In addition, (5)–(6) hold the same form not only at a specific operating point, but also for a range of the system load, as long as there is no change in the binding/unbinding constraint set, e.g., within any two neighboring CLLs. Furthermore, participation factor is introduced to model the load variation pattern. It is defined as follows [3]: (7) where

is the total system load; ; and superscript “(0)” denotes the present or given operating point. Substituting (7) into (5)–(6) gives the characteristic constraints w.r.t. the total system load

III. CHARACTERISTIC CONSTRAINTS OF THE MODEL Among the constraints (2)–(4), the energy balance equation (2) is an equality constraint while the line limits in (3) and generation output limits in (4) are inequality constraints. When the load varies and reaches a CLL, several of the inequalities may become equalities. In other words, certain unbinding constraints may become binding, such as an uncongested line flow reaching its limit or a marginal generation output reaching its upper/lower limit. Furthermore, when the system load varies within the same CLL interval, namely, between the same two adjacent CLLs around the given operating point, all binding constraints will remain binding. For instance, the non-marginal unit set will remain unchanged and the congested lines will remain congested. The binding constraints are called “characteristic constraints” in this paper. The corresponding equalities are written as follows:

(5)

(6)

where represents the set of congested line(s). ), the Note that when losses are ignored (i.e., with-loss-DCOPF model will be simplified to a lossless DCOPF model, which is a linear programming problem. For a lossless DCOPF model, it can be easily proved that the number of equations in (5)–(6) equals the number of variables . This means that the number of marginal is equal to the number of congested lines units plus 1, i.e., [2]. Therefore, only one solution exists for (5)–(6) for a lossless DCOPF. It implies that the optimal solution can be effectively determined by only the binding constraints, because the feasible region constrained by all binding constraints boils down to a single point. Therefore, for the lossless DCOPF, (5)–(6) are essentially equivalent to . the lossless DCOPF model in (1)–(4) with The above statements also hold true for the with-loss-DCOPF model because of the roughly quadratic nature of the losses.

(8)

(9) By introducing the participation factor , a set of variables, , is converted to a single variable, , such that we can focus on the load variation at the system level. The participation factor is normally defined to be a constant vector to model the load variation pattern around a given operating point and can be assigned different values over time or on a case-dependent basis for various study purposes. Note that the CLLs may change when a different load variation pattern is assumed. One may argue that each load can vary independently in theory and therefore a deterministic load variation pattern should never exist. It is however not the case in practice. In power system operation, loads at the aggregated level do not change independently. Rather, they are correlated through a number of factors such as customer usage pattern, temperature, time of day, and day of week. In planning, bus-level loads under a control area are normally assumed to vary in the same pattern. A pre-defined load variation pattern has been widely employed in various industry practices such as continuation power flow [16], power transfer capability analysis, and long-term economic study. IV. VARIABLE SUBSTITUTION METHOD In (8)–(9), the independent variable is the total system load, , and the dependent variables are . It is beneficial to w.r.t. the variation of . It is deknow the trajectory of sirable to solve these algebraic equations for the marginal unit . However, different from the linear output as a function of lossless model in [3], it is hard to derive an analytical form due

BO AND LI: EFFICIENT ESTIMATION OF CRITICAL LOAD LEVELS USING VARIABLE SUBSTITUTION METHOD

to

and , which present the nonlinearity of w.r.t. . An alternative approach is to obtain an approximate solu. Since the amount of power losses is roughly a tion for quadratic function of the load and is compensated by the marginal units, as seen in the power balance equation, it is intuitive to assume the marginal unit generation follows a quadratic polynomial pattern. In fact, a quadratic pattern has been reported to be a good approximation of the marginal unit generation even under the ACOPF dispatch framework [8]. It has been proved to be a strictly quadratic pattern at conditions when there is a marginal unit at the reference bus [19]. Therefore, we define as to approximate follows a quadratic polynomial function trajectory w.r.t. the load variation: the

(10) where denote the 2nd degree, 1st degree coefficients, and constant part of the th marginal unit generation, respectively. These quadratic coefficients can be obtained in three steps summarized as follows. The first step is variable substitution. Substituting (10) into the left-hand side of (8) yields a funcand , which is denoted tion w.r.t. by . Likewise, the set of functions obtained by substituting (10) into the left-hand side of (9) is denoted by . After the substitutions, the variables are replaced with . Let us a new set of variables such that the assume that there exists a set of satisfies (8)–(9). It corresponding approximation of implies (11) (12) The second step is to construct a sufficient number of independent equations in order to solve for the coefficients. Equations (11)–(12) are expected to hold true at the given operating point and any other load level close enough to as long as no change of binding constraints will occur at that load level. As such, two arbitrary load levels, and , are selected and (11)–(12) are evaluated at the three load levels, , and , respectively. Therefore, a set of nonlinear equations are established as follows:

2475

(16) (17)

(18) The reason that (11)–(12) are evaluated at three load levels is to construct enough equations so that the number of equalities equals the number of variables, i.e., the quadratic polynomial coefficients. and can be any arbitrary Note also that, in theory, load levels even if they are actually outside the range bounded by . the previous CLL and the next CLL of a given load level This is due to the fact that (11)–(12) should always hold true for any load level if we hypothetically remove all unbinding constraints and keep only the binding constraints in the OPF model, i.e., no change of binding constraint will ever occur. Ideally, we and to be near the to-be-predicted want to choose CLL so as to minimize possible numeric errors in computation. However, CLL is yet unknown and to be estimated. In this paper, and are selected to be close to . For example, and . The marginal unit generation output at the given operating point should not be added into the set of equations due to its dependence on (13)–(14). It should also be stressed that the realization of (13)–(18) does not need any information at load levels and such as generation dispatch except the values of the two load levels. Therefore, no additional OPF runs in any form are needed. The third step is to solve the nonlinear system of equations for the unknowns, namely, the quadratic polynomial coefficients. Equations (13)–(18) are nonlinear functions of the new vari. The number of equations is ables ) which is equal to . Therefore, the set the number of variables of equations could be solved by a common nonlinear equation solving algorithms such as Newton’s method, Gauss-Newton is method, and Levenberg-Marquardt method. Typically, a very small number. Thus, the above nonlinear system is a very unknowns. small-scale problem with only Finally, with solved coefficients of marginal unit dispatch, the CLLs can be easily estimated. We can extrapolate the quadratic curve characterized by these coefficients and predict the load levels where the marginal units will reach their upper or lower limit [8]. For instance, solving univariate quadratic equations which denotes the distance from the present (19) gives load level to the load level where the th marginal unit will reach its upper limit

(13)

(14) (15)

(19) The line flow on a non-congested line is a linear function of the marginal unit generation and therefore is also a quadratic function of the total

2476

Fig. 2. System topology of a revised IEEE 30-bus system.

system load, . Extrapolating the quadratic line flow pattern and solving it at the maximum line limit will provide the load level distance between the present load level and the load level where a line flow will reach its limit. The complete process of using the variable substitution method to predict congestions and CLLs is as follows. 1) Establish (13)–(18) and solve for coefficients of the approximated marginal unit generation, . 2) Extrapolate the quadratic representation of the marginal unit outputs to find the load level distance between the present load level and the level where a new binding constraint occurs. Perform a similar calculation for line flow through each non-congested transmission branch [8]. 3) Find the largest negative load level and smallest positive load level. These two load levels are the previous CLL and the next CLL , respectively [8]. Steps 2)–3) are the same procedures as those for the ACOPF framework reported in [8]. Note that new congestions can be identified simultaneously when a CLL is identified. It should be stressed that the proposed method does not utilize any modelspecific features and it is easy to implement. Therefore, it is a potentially universal approach and may be applied to various types of DCOPF-based nonlinear models. V. NUMERIC STUDY ON A REVISED IEEE 30-BUS SYSTEM The proposed method will be illustrated on a revised IEEE 30-bus system [17], [18], which contains realistic transmission line limits, bus area assignment, and relocated generators. The detailed system data can be found in [18]. The system topology is shown in Fig. 2. The biddings of the six generators are set at 10, 15, 30, 35, 40, and 45, respectively, all in $/MWh. A. Comparison of Quadratic Approximation and Linear Approximation The method is evaluated on the conventional with-lossDCOPF model presented in Section II. Fig. 3 compares the benchmark data of the marginal unit generation at Bus 22 with the corresponding quadratic approximation results from the proposed method, which is performed at 209 MW, the base-case load level. The benchmark data can be obtained from

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 4, NOVEMBER 2011

Fig. 3. Benchmark data, quadratic approximation, and linear approximation of the generation output of marginal unit at Bus 22.

Fig. 4. Differences in percentage of the quadratic approximation and the linear approximation of the marginal unit output at Bus 22 w.r.t. the benchmark data.

repetitive DCOPF runs in the load range from 209 MW to 234 MW. In this load range, Line 6–8 is always congested and the marginal units are the generators at Buses 22 and 27. As a comparison, Fig. 3 shows the results of a linear approximation approach, which ignores the power losses and yields a linear representation for the marginal unit generation. The differences between the benchmark data and approximation results are shown in Fig. 4. Figs. 3 and 4 suggest that the quadratic approximation achieved by the proposed method is a reliable approximation to the benchmark data, and outperforms linear approximation. Although it appears in Fig. 3 that linear approximation is still an acceptable candidate, it may not work in other cases where outstanding nonlinearity of system statuses is observed. For instance, Fig. 5 shows the benchmark data and the quadratic approximation results of the power flow through Line 1–2, as well as their differences in percentage. Evidently, linear approximation of the nonlinear line flow may not yield results of acceptable accuracy. In Fig. 5, the differences appear to be increasing when the load grows. Nevertheless, the quadratic approximation is still highly reliable, given that the difference between the benchmark data and quadratic approximation results is less than 0.05%, while the load variation spans 12% of the base case load (209 MW). Figs. 3–5 illustrate the high accuracy of the results from

BO AND LI: EFFICIENT ESTIMATION OF CRITICAL LOAD LEVELS USING VARIABLE SUBSTITUTION METHOD

2477

TABLE I PREVIOUS AND NEXT CRITICAL LOAD LEVELS FOR THE GIVEN OPERATING POINT 209 MW

Fig. 5. Quadratic approximation, benchmark data, and their percentage differences of the power flow through Line 1–2.

the proposed quadratic approximation approach. Similar patterns are observed for simulations at other load intervals which do not contain a CLL. It should be stressed that the benchmark data at a number of load levels, as shown in Figs. 3 and 5, are only used for comparison purposes, and not inputs to the proposed method. This is completely different from other methods such as the curve-fitting or interpolation method which requires OPF solutions at a few load levels such that the solution trajectory can be estimated. On the contrary, the proposed method does not need OPF runs at any other load levels. For instance, in the studied case, the proposed method needs only the information of the marginal unit and congested line sets at the given load level 209 MW, around which the previous CLL and next CLL are to be estimated. Equations (13)–(18) are then established, which consist nonlinear equations as there is only one of congested line and correspondingly two marginal units. The two marginal units have totally 6 unknown coefficients. Likewise, if there are two congested lines, the total number of equations will while there are 9 unknown coefficients for be 3 marginal units. Thus, an increase in power system size does not necessarily lead to a larger system of nonlinear equations because it depends on the number of congestions rather than the system size. With coefficients obtained from the solution of nonlinear equations (13)–(18), the coefficients of all the non-congested lines can be subsequently obtained. The resulting approximation of both marginal unit generations and line flows are amazingly good fits to all the actual benchmark data at various load levels as well as the given operating point. B. Verification of the Proposed Method The CLL estimation is performed at the 209 MW load level. Table I shows a comparison between the previous and next CLLs predicted from the proposed method and the actual benchmark values. The actual CLLs can be obtained from repetitive withloss-DCOPF simulations or the binary search method with 0.01 MW precision. New binding constraints at the corresponding CLLs, as correctly identified by the proposed method, are also shown in Table I.

Table I shows that, when compared with the benchmark results obtained from the repetitive DCOPF runs, the estimated CLLs from the proposed method are highly accurate while involving far less computational efforts. (Note: Additional discussions about computational complexity analysis and test results can be found in Section VII-C.) In addition, new binding constraints are identified along with the estimation of CLLs. For instance, by the proposed method, Line 25–27 will reach its flow limit in the reverse direction if the system load increases from the given load level to 234.32 MW, which is very close to the benchmark value of 234.29 MW. At this CLL, there is a new marginal unit at Bus 13. In the case of a load decrease, the first change of binding constraints is expected to occur when the load decreases to 208.46 MW, at which point the marginal unit at Bus 27 will reach its lower limit and is no longer a marginal unit. The disappearance of a marginal unit is due to the relief of congestion through Line 6–8. This example illustrates the ability of the proposed method in congestion prediction. For any given operating point, we can obtain an estimated previous CLL and next CLL by applying the proposed method. Table II shows the estimated previous and next CLLs for a selection of given load levels ranging from 0.6 p.u. to 1.3 p.u. of the base case load, i.e., from 113.52 MW to 245.96 MW. The given load levels are not selected to be evenly distributed in the studied load range because in that way, certain given load levels may be within the same segment bounded by two consecutive CLLs. Thus, the predicted CLLs will be the same. For instance, assume that the given load levels are 0.9 p.u. of the base case load (i.e., 170.28 MW) and 1.0 p.u. (i.e., 189.20 MW), both of which reside in the segment bounded by two adjacent CLLs, 155.67 MW and 194.99 MW. The same previous CLL and next CLL will be found for the cases of given load level 0.9 and 1.0 p.u. Also, a CLL may be missed if none of the given load levels falls into the two segments associated with the CLL. Therefore, for illustrative purpose, one and only one load level is picked from each segment demarcated by CLLs. For each CLL, the estimated value from its left side may differ slightly from the value estimated from its right side. In this regard, the estimated CLLs in Table II are rearranged, as shown in Table III, such that each actual CLL can be conveniently compared with the values estimated from its left and right, respectively. Tables II and III demonstrate the high accuracy of the proposed method in estimating CLLs. The maximum error of all estimated CLLs is 0.1 MW, or 0.05% of the base case load. In addition, the CLLs estimated from left and right, respectively, are very close in values and thus suggest the consistency of the

2478

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 4, NOVEMBER 2011

TABLE II PREVIOUS AND NEXT CRITICAL LOAD LEVELS ESTIMATED FROM DIFFERENT GIVEN OPERATING POINTS

TABLE IV COMPARISONS OF ACTUAL CLLS AND ESTIMATED CLLS FOR MODEL M-1 ON THE IEEE 118-BUS SYSTEM

TABLE III COMPARISONS OF ACTUAL CLLS AND ESTIMATED CLLS

a $0.5 increment; the next more expensive 20 generators have bids from $30 to $49 with a $1 increment; and the remaining most expensive 14 generators have bids from $70 to $83 with a $1 increment. To illustrate that the proposed approach is not restricted to a specific form of DCOPF-based modeling, the approach will be tested using four different with-loss-DCOPF models. One of the models is a recent fictitious nodal demand (FND)-based DCOPF model [11] which converts branch losses into fictitious demand at the two ends of a branch and correspondingly incorporates power losses in power flow formulations. In other words, (3) is revised to include a nonlinear item to model loss effect on line flows. Due to power losses being more reasonably distributed in the system, this model may yield more accurate power flows, generation dispatches, and LMP than the conventional DCOPF model. The four models are summarized as follows. • The first model, M-1, is the conventional with-lossDCOPF model, as presented in Section II. • The second model, M-2, is the FND-based DCOPF model. • The third model, M-3, is a variant of M-1, with (2) being substituted by (20), an alternative form of energy balance equation. • The fourth model, M-4, is a variant of M-2, with (2) being substituted by (20):

results from the proposed method. Note that the results are numerically stable to the given operating point. Numeric study has verified that slight variations of any given operating point in the first column of Table II essentially lead to no difference in the estimated CLLs. Also note that Table II implies a useful application of the quick estimation of OPF solutions including generation dispatch and line flows. For instance, for any interested load level in the range of [78.37 MW, 155.68 MW], no OPF re-run is needed. Instead, the solution of the with-loss-DCOPF model can be quickly estimated by (10) with the assumed load variation pattern. VI. NUMERIC STUDY ON THE IEEE 118-BUS SYSTEM In this section, the proposed method will be tested on the IEEE 118-bus system. The system consists of 118 buses, 54 generators, and 186 branches/transformers. System total load is 4242 MW with 9966.2 MW of total generation capacity. The detailed system data is available in [17]. The OPF related data is derived from [6]. Five line limits are added into the transmission system: 345 MW for Line 69–77, 630 MW for Line 68–81, 106 MW for Lines 83–85, and 94–100, 230 MW for Line 80–98. For illustrative purposes, the generators are sorted in the order of the bus numbers where they are connected. The first 20 generators are assumed to be the cheapest with bids from $10 to $19.5 with

(20)

First, the binary search method is applied for each of the four models to find all CLLs in the load level range between 3393.6 MW to 5090.6 MW (i.e., 0.8 p.u. to 1.2 p.u. of the base-case load) to the precision of 4.242 MW (i.e., 0.001 p.u.). An iterative linear programming (LP) algorithm as presented in [11] is adopted for solving the four DCOPF models. Second, the variable substitution method is employed to estimate CLLs. The actual CLLs and its estimations using models M-1, M-2, M-3, and M-4 are shown in Tables IV–VII. Tables IV and V show that the average estimation error is 1.57 MW for M-1 and 1.65 MW for M-2, and the maximum estimation error is 3.31 MW for M-1 and 2.84 MW for M-2. The errors are quite trivial given that the precision of the actual

BO AND LI: EFFICIENT ESTIMATION OF CRITICAL LOAD LEVELS USING VARIABLE SUBSTITUTION METHOD

TABLE V COMPARISONS OF ACTUAL CLLS AND ESTIMATED CLLS FOR MODEL M-2 ON THE IEEE 118-BUS SYSTEM

CLL is 4.242 MW. Therefore, any estimated CLL with an error less than 4.242 MW should be viewed as technically accurate. In Table IV, it can be seen that the distances between any two consecutive CLLs may vary widely. For instance, the smallest distance is only 29.16 MW, or 0.0069 p.u. of the base case load, between CLLs at 4630.56 MW and 4659.72 MW. The segment bounded by these two CLLs is so small that any possible pattern change in dispatch or LMP in the segment will likely be missed by repetitive DCOPF simulations with a sparser sampling resolution, such as 0.01 p.u., of the base case load. In the opposite scenario, the distance between CLLs at 3642.71 MW and 4081.62 MW is 438.91 MW, a little more than 0.1 p.u., or 10% of the base case load. This segment between the two CLLs is so wide that numerous repetitive DCOPF runs can be avoided because existing studies have shown that the variation pattern of dispatch and LMP do not change in any segment bounded by two consecutive CLLs [8], [11]. For example, the LMPs follow a linear pattern in the segment. Therefore, given the LMP at any given operating point in the segment and the calculated LMP sensitivity [4], we can predict the LMP for any load level in the segment without additional DCOPF runs. The knowledge of CLL locations indicates how far we can apply the local LMP sensitivities to a wider interval without considerable loss of accuracy. Tables VI and VII show that the average estimation error is 1.23 MW for M-3 and 1.12 MW for M-4, and the maximum estimation error is 3.43 MW for M-1 and 9.12 MW for M-2. Again, the errors are very small and even the largest error 9.12 MW is only 0.00215 p.u. of the base case load. Therefore, the estimated CLLs can be considered as highly accurate. In summary, simulation results in Tables IV–VII have consistently demonstrated the high accuracy achieved by the proposed method on the IEEE 118-bus system. It can be seen from Tables IV–VII that the number and the locations of actual CLLs differ considerably among the four models. This occurs when high system losses are present. For example, in the base case of the IEEE 118-bus system, total system losses account for approximately 15% of the total load due to high resistance and high line flow. The losses are so high that they may have a significant impact on the generation dispatch. Because the four models employ different approaches in modeling power losses in the energy balance and line flows

2479

TABLE VI COMPARISONS OF ACTUAL CLLS AND ESTIMATED CLLS FOR MODEL M-3 ON THE IEEE 118-BUS SYSTEM

TABLE VII COMPARISONS OF ACTUAL CLLS AND ESTIMATED CLLS FOR MODEL M-4 ON THE IEEE 118-BUS SYSTEM

equations, shifting from one model to another for the same given load level will likely lead to considerable changes in generation dispatch. The marginal unit generation trajectories are shifted and subsequently the load level where the trajectories hit their boundaries may change. This may result in considerable change in CLL locations. On the contrary, the system loss level is approximately 2.6% of the total load at the base case for the IEEE 30-bus system, and CLLs for different models are quite close. In addition, numeric studies have verified that the CLLs for the four models are much closer for the IEEE 118-bus system with a significantly lower level of losses. Nevertheless, this study leads to an interesting topic about which with-loss-DCOPF model presented here as well as in the literature is the most appropriate representation of the actual AC system. Apparently, this is a possible future work. VII. DISCUSSIONS Some unique features of the proposed variable substitution method will be discussed in this section, including the model adaptability, easy implementation and expandability, and high computational efficiency. A. Model Adaptability As presented in Section VI, the proposed method has been tested on four different DCOPF models. It demonstrates the

2480

wide applicability of the method. In fact, the method is a generic approach for estimating the three coefficients of the assumed quadratic form of every marginal unit generation. It does not utilize any model specific feature and therefore is not restricted to a specific form of DCOPF modeling. For example, the transmission losses can be modeled at the reference bus like in M-1 and M-3, or distributed in the transmission lines like in M-2 and M-4. Therefore, it has the potential of working with any with-loss-DCOPF model. B. Easy Implementation and Flexibility The proposed method is very easy to implement. It needs only a small piece of information as inputs, namely, the set of marginal unit and the set of congested branches at the given operating point only. Then, nonlinear equations (13)–(18) can be established and solved for the quadratic polynomial coefficients. Note that the system of nonlinear equations may be solved without deriving explicit formulas of derivatives. Numerical approaches such as finite-differencing derivatives can be used to ease the use of the solvers. In this case, there is even no need to explicitly derive the nonlinear equations in terms of the quadratic coefficients. In fact, the nonlinear equations could be kept in their original form, namely, as functions of intermediate . This makes the implementation of the proposed variables method even simpler. Although a quadratic polynomial is the main focus of the proposed method, it allows use of higher order polynomials to approximate the trajectory of the DCOPF solutions w.r.t. load variations. This can be easily implemented. For example, if a third-order polynomial is considered, we can simply establish an additional set of nonlinear equations by evaluating (11)–(12) at a new hypothetical load level such as . Then, the new equations, combined with the existing equations, can be used to solve the coefficients. C. Computational Efficiency Speed is always a concern in short-term or online applications. It is even more demanding when multiple load variation patterns need to be evaluated. In this section, the computational efficiency of the proposed method will be compared with binary search and a more efficient method, the three-point interpolation method [8]. The binary search algorithm is a commonly used algorithm in locating the position of items in a sorted array. On average, it times of iterations before locating the item, takes is the total number of items. When it is applied to where find CLLs, the number of items depends on the length of the searched interval and precision requirement. For instance, if there is one and only one CLL in a 100 MW interval with a required precision of 0.1 MW, the total number of possible load . Therefore, the average levels will be . For each number of iterations needed is iteration, a with-loss-DCOPF model needs to be solved. Given the solution at the present operating point, there are 8 additional with-loss-DCOPF runs to be executed. In comparison with binary search, the three-point interpolation method is much more efficient simply because it generally requires as few as two additional with-loss-DCOPF runs. The

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 4, NOVEMBER 2011

TABLE VIII AVERAGE SPEEDUP OF THE PROPOSED METHOD COMPARED WITH BINARY SEARCH METHOD AND THREE-POINT INTERPOLATION METHOD

solutions of the two with-loss-DCOPF runs as well as the run at the given operating point will be used to perform quadratic interpolation to acquire the coefficients of quadratic polynomials [8]. The polynomials will be extrapolated to estimate the CLLs. Therefore, the three-point interpolation method saves a number of with-loss-DCOPF runs in the expense of solving a linear system of equations for quadratic coefficients and solving a set of univariate quadratic equations. The proposed variable substitution method is even more superior to the three-point interpolation method because it requires no additional DCOPF runs other than the solution at the given operating point. The quadratic coefficients can be obtained by solving a small-scale nonlinear system of equations. Therefore, the variable substitution method saves the time of at least two additional with-loss-DCOPF runs in the expense of solving a small-scale nonlinear system of equations. Solving the nonlinear system of equations generally requires much less effort than solving the with-loss-DCOPF model which is a nonlinear programming (NLP) problem. The with-loss-DCOPF model can be solved iteratively, and in each iteration, a LP model needs to be solved. A well-known algorithm for solving LP is Karmarkar’s algorithm, an interior point , where method. Its performance complexity is up to is the number of variables in its special canonical LP form. Generally, is far greater than the total number of generators. In contrast, the nonlinear system of equations can be solved by Newton’s method which linearizes the nonlinear equations and solves the linear system of equations in iterations. In each iteration, the linear system of equations can be solved by direct methods such as LU factorization. The time complexity for LU , where is the size of the problem. factorization is . In this study, the size of the problem is Given that particularly true for large systems, the computational complexity of the variable substitution method is much smaller than that of the binary search method as well as the three-point interpolation method. With the estimated number of with-loss-DCOPF runs as mentioned above, the average speedup of the proposed variable substitution method against binary search and three-point interpolation is shown, respectively, in Table VIII. The speedup concept in [3] is adopted here and defined as the running time of the binary search or three-point interpolation method divided by the running time of the proposed method. While there is no, or limited, efficiency gains for the small PJM five-bus system, significant speedup is observed for the IEEE 118-bus system, namely, 5.0 and 18.2 in comparison with three-point interpolation and binary search, respectively. More important, the speedup increases with the size of the study system. In

BO AND LI: EFFICIENT ESTIMATION OF CRITICAL LOAD LEVELS USING VARIABLE SUBSTITUTION METHOD

addition, the speedup of the proposed method will be even greater in cases when higher precision is required for binary search, and when more than two additional with-loss-DCOPF runs are necessary for the three-point interpolation method (e.g., one or two of the selected trial load levels lies beyond the CLL to be predicted). Note that the speedup may vary with the models, algorithms, and implementation techniques which affect the number of iterations before convergence. VIII. CONCLUSION This paper firstly reviews the with-loss-DCOPF model, which is a nonlinear programming model due to the incorporation of power losses. When the total system load is in the range bounded by two adjacent CLLs, the optimal solution of the with-lossDCOPF model can be essentially characterized by the binding constraints, which are called characteristic constraints in this paper. These equations are employed to identify the trajectory of the optimal solution with respect to load variation. By introducing the load variation participation factor, the characteristic equations can be rewritten as a set of implicit functions of the marginal unit generation w.r.t. only one variable, the system load. Based on intuition and the observations from existing literatures which show that the marginal unit generations are approximately quadratic to the system load, quadratic polynomials are employed to model the marginal unit outputs in the characteristic constraints. Then, nonlinear equations based on the characteristic constraints will be solved to obtain the marginal units’ quadratic polynomial coefficients w.r.t. the total system load. The nonlinear equations can be easily and efficiently solved due to their small dimension. The quadratic coefficients are used to perform trajectory extrapolation, and subsequently, to identify congestions and the previous/next CLLs around the given operating point. The obtained quadratic trajectory can be directly applied for a quick estimation of the OPF solutions, including generation dispatch and line flows, at any load level around a given operating point, as long as the projected load level is within the range bounded by the previous and next CLLs. Thus, no OPF re-run is needed and significant saving in computation is achieved. A case study on a revised IEEE 30-bus system demonstrates the effectiveness of the proposed method. Among all studied scenarios, the maximum error in estimated CLLs is 0.1 MW, or 0.05% of the base case load. Potential congestion and other binding constraints are identified along with the estimation of CLLs. The proposed variable substitution method is also tested on the IEEE 118-bus system for four different DCOPF models including the conventional with-loss-DCOPF model and the FND-based DCOPF model. For each model, every single CLL has been identified with high accuracy. In addition to its high accuracy, the proposed method is also very efficient. Different from existing methods such as the interpolation method in [8] which utilizes OPF solutions at two other load levels beside the given one, the proposed method does not need any information other than those at the given load level, and therefore, requires no additional OPF runs. Coefficients of the quadratic approximations are obtained by solving a typically

2481

very small-scale system of nonlinear equations. Furthermore, the proposed method is not restricted to any specific formulation of the DCOPF model, and thus can be easily applied to variants of the with-loss-DCOPF model. This is demonstrated with the IEEE 118-bus system. The computational efficiency has been discussed through a time complexity analysis. It is shown to be much more efficient than existing methods such as binary search and interpolation. Also, the proposed method is very easy to implement and can be easily adapted to support the high-order polynomials, if indeed necessary, for a very large system. In conclusion, the proposed method is a highly efficient and effective method for estimation of CLLs and prediction of congestion. It is model independent and easy to implement, and therefore, applicable to various nonlinear with-loss-DCOPF models. REFERENCES [1] L. John et al., “The view from the top,” IEEE Power Energy Mag., vol. 7, no. 6, pp. 76–88, Nov./Dec. 2009. [2] F. Li, “Continuous locational marginal pricing (CLMP),” IEEE Trans. Power Syst., vol. 22, no. 4, pp. 1638–1646, Nov. 2007. [3] F. Li and R. Bo, “Congestion and price prediction under load variation,” IEEE Trans. Power Syst., vol. 24, no. 2, pp. 911–922, May 2009. [4] A. J. Conejo, E. Castillo, R. Minguez, and F. Milano, “Locational marginal price sensitivities,” IEEE Trans. Power Syst., vol. 20, no. 4, pp. 2026–2033, Nov. 2005. [5] X. Cheng and T. J. Overbye, “An energy reference bus independent LMP decomposition algorithm,” IEEE Trans. Power Syst., vol. 21, no. 3, pp. 1041–1049, Aug. 2006. [6] R. Bo and F. Li, “Probabilistic LMP forecasting considering load uncertainty,” IEEE Trans. Power Syst., vol. 24, no. 3, pp. 1279–1289, Aug. 2009. [7] Q. Zhou, L. Tesfatsion, and C.-C. Liu, “Global sensitivity analysis for the short-term prediction of system variables,” in Proc. 2010 IEEE Power and Energy Soc. General Meeting, Minneapolis, MN, Jul. 2010. [8] R. Bo, F. Li, and C. Wang, “Congestion prediction for ACOPF framework using quadratic interpolation,” in Proc. IEEE PES General Meeting 2008, Pittsburgh, PA, 2008. [9] B. Stott, J. Jardim, and O. Alsac, “DC power flow revisited,” IEEE Trans. Power Syst., vol. 24, no. 3, pp. 1290–1300, Aug. 2009. [10] Y. Fu and Z. Li, “Different models and properties on LMP calculations,” in Proc. 2006 IEEE Power Eng. Soc. General Meeting, Montreal, QC, Canada, Oct. 2006. [11] F. Li and R. Bo, “DCOPF-based LMP simulation: Algorithm, comparison with ACOPF, and sensitivity,” IEEE Trans. Power Syst., vol. 22, no. 4, pp. 1475–1485, Nov. 2007. [12] A. L. Ott, “Experience with PJM market operation, system design, and implementation,” IEEE Trans. Power Syst., vol. 18, no. 2, pp. 528–534, May 2003. [13] PJM Training Materials (LMP101), PJM. [Online]. Available: http:// www.pjm.com/training/training-material.aspx (accessed in Oct. 2010). [14] NYISO Transmission & Dispatch Operation Manual, NYISO. [Online]. Available: http://www.nyiso.com/public/webdocs/documents/manuals/operations/trans_disp.pdf (accessed in Oct. 2010). [15] Business Practices Manual: Energy and Operating Reserve Markets. [Online]. Available: http://www.midwestiso.org/publish (accessed in Oct. 2010). [16] M. Crow, Computational Methods for Electric Power Systems. Boca Raton, FL: CRC, 2002. [17] Power System Test Case Archive,University of Washington. [Online]. Available: http://www.ee.washington.edu/research/pstca/ (accessed in Oct. 2009). [18] R. D. Zimmerman, C. E. Murillo-Sánchez, and D. Gan, MatPower-A Matlab Power System Simulation Package, School of Electrical Engineering, Cornell Univ.. Ithaca, NY. [Online]. Available: http://www. pserc.cornell.edu/matpower/matpower.html (accessed in Oct. 2010). [19] R. Bo and F. Li, “Marginal unit generation sensitivity and its applications in transmission congestion prediction and LMP calculation,” in Proc. IEEE PES Power Systems Conf. Exhib. (PSCE) 2011, Phoenix, AZ, Mar. 20–23, 2011.

2482

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 26, NO. 4, NOVEMBER 2011

Rui Bo (S’02–M’09) received the B.S. and M.S. degrees in electric power engineering from Southeast University, Nanjing, China, in 2000 and 2003, respectively, and the Ph.D. degree at The University of Tennessee, Knoxville (UTK), in 2009. He worked at ZTE Corporation and Shenzhen Cermate Inc. from 2003 to 2005, respectively. He has been employed at Midwest Independent Transmission System Operator (Midwest ISO), St. Paul, MN, since 2009. His interests include power system operation and planning, power system economics, and computational methods.

Fangxing (Fran) Li (M’01–SM’05) received the B.S. and M.S. degrees in electric power engineering from Southeast University, Nanjing, China, in 1994 and 1997, respectively, and the Ph.D. degree from Virginia Tech, Blacksburg, in 2001. He has been an Assistant Professor at The University of Tennessee (UT), Knoxville, TN, since August 2005. Prior to joining UT, he worked at ABB, Raleigh, NC, as a senior and then a principal engineer for 4 and a half years. During his employment at ABB, he was the lead developer of GridView ™, ABB’s market simulation tool. His current interests include energy market, reactive power, distributed energy control, and renewable energy. Dr. Li is a registered Professional Engineer (P.E.) in North Carolina.