Probabilistic Forecast of Real-Time LMP and Network Congestion
aa r X i v : . [ s t a t . A P ] J un Probabilistic Forecasting of Real-Time LMP andNetwork Congestion
Yuting Ji,
Student Member, IEEE,
Robert J. Thomas,
Life Fellow, IEEE, and Lang Tong,
Fellow, IEEE
Abstract —The short-term forecasting of real-time locationalmarginal price (LMP) and network congestion is considered froma system operator perspective. A new probabilistic forecastingtechnique is proposed based on a multiparametric programmingformulation that partitions the uncertainty parameter space intocritical regions from which the conditional probability distribu-tion of the real-time LMP/congestion is obtained. The proposedmethod incorporates load/generation forecast, time varying oper-ation constraints, and contingency models. By shifting the com-putation cost associated with multiparametric programs offline,the online computation cost is significantly reduced. An onlinesimulation technique by generating critical regions dynamicallyis also proposed, which results in several orders of magnitudeimprovement in the computational cost over standard MonteCarlo methods.
Index Terms —Congestion forecast, electricity price forecast,locational marginal price (LMP), multiparametric programming,probabilistic forecast.
I. I
NTRODUCTION
As more renewable resources are integrated into the powersystem, and the transmission system operates closer to itscapacity, congestion conditions become less predictable andlocational marginal prices (LMP) more volatile. The increasedcongestion and LMP uncertainties pose significant challengesto the operator and market participants, which motivates usto consider the problem of short-term forecasting real-timeLMP/congestion in the presence of generation, demand, andoperation uncertainties.The benefit of accurate LMP and congestion forecasts istwofold. For market participants, accurate forecast of real-time prices is valuable in risk management, bidding strategydevelopment, and demand side participation. The forecastprices allow market participants to make adjustments in ad-vance to ensure competitive transactions. For system operators,on the other hand, forecast of transmission congestion isimportant for congestion management and system planning.European transmission system operators, for instance, useIntraday Congestion Forecast (IDCF) to improve real-timesecurity assessment [2]. Intuitively, an LMP forecast shouldelicit generation participation at times of potential shortagethus alleviating future congestions. Similarly, an LMP forecastcan be used for demand response that results in shifting part
The authors are with the School of Electrical and Computer Engineering,Cornell University, Ithaca, NY 14853, USA (e-mail: [email protected];[email protected]; [email protected]).This work is supported in part by the U.S. Department of Energy Office ofElectricity Delivery and Energy Reliability Consortium for Electric ReliabilityTechnology Solutions (CERTS) program and the National Science Foundationunder Grant CNS-1135844. Part of this work appeared in [1]. of the load from peak to valley. An example of applying LMPforecast for inter-regional interchange is presented in [3].Currently, some system operators are providing real-timeprice forecasts. The Electric Reliability Council of Texas(ERCOT) [4] offers a 1-hour ahead real-time LMP forecast,updated every 5 minutes. The Alberta Electric System Oper-ator (AESO) [5] provides two short-term price forecasts withprediction horizons of 2 hours and 6 hours, respectively.Most LMP forecasting schemes provide only point forecast ,which gives a single quantity as the prediction. For systemswith high levels of uncertainty, a point forecast is rarely accu-rate, and impacts of prediction error are difficult to quantify.A more attractive alternative is a probabilistic forecast thatprovides a full characterization of the LMP distribution.Significant technical challenges exist for probabilistic fore-casting real-time LMP and congestion. First, reasonably accu-rate models for real-time dispatch and LMP are needed. Sec-ond, network operating conditions and uncertainties need tobe incorporated in real time. Finally, the forecasting algorithmneeds to be simple and scalable for sufficiently large systems.These challenges are daunting for external market partici-pants who do not have access to network operating conditionsand confidential information on bids and offers that influenceLMPs. On the other hand, if it is the system operator providingthe forecast, as in the case of ERCOT or AESO, the barrierto efficient and accurate forecast is lowered.
A. Summary of Contributions
In this paper, we consider the real-time LMP and congestionforecasting problem from an operator perspective. We focuson probabilistic forecasting that, at time t , provides the con-ditional probability distribution at time t + T of the LMPvector and associated congestion, given the system state attime t . Here T is referred to as the prediction horizon whichis considered in the range of T from 1 to 6 hours for short-termforecasts.The key idea behind the proposed approach is the use ofmultiparametric programming that partitions the uncertaintyspace into critical regions with each region attached to aunique LMP and congestion pattern. Thus, the problem ofprobabilistic forecasting reduces to computing the probabilitiesof random parameters falling in the set of critical regions.When loads or generations (treated as negative loads) are ran-dom, their forecasts are incorporated to generate probabilisticLMP and congestion forecasts.The proposed technique also provides several new fea-tures not present in existing methods. For example, it can incorporate system contingency models and allow systemconstraints to vary with time. The latter feature is relevantbecause network topology and thermal limits may be changedin real time by the system operator depending on operatingconditions. In terms of the generation cost, the proposal canbe applied to a linear (or piece-wise affine) function anda quadratic function. Computationally, the proposed methodshifts a majority of the computation offline, which significantlyreduces online computation cost.An alternative algorithm that dynamically generates criticalregions is also proposed. Because load and stochastic gen-eration processes are physical processes, they are boundedand tend to concentrate around the mean trajectory. Theirrealizations thus only fall in a few critical regions instead ofall over the entire parameter space. By generating the criticalregions that contain such realizations, the computation costis reduced by several orders of magnitude comparing withstandard Monte Carlo techniques. B. Related Work
Much of the existing work deals with point forecasts ofLMP by market participants who do not have access to real-time operating conditions and confidential offers and bids. Forthese techniques, historical data on LMP, load, and congestionsdrive the forecasting engine. Literature on these techniquesabounds. See [6] and references therein. For probabilisticforecasting techniques by market participants, see approachesin Global Energy Forecasting Competition 2014 [7].There are several prior studies on LMP/congestion fore-casting from the system operator perspective. The proposedtechnique in [8] employs an online Monte Carlo samplingtechnique that, for each Monte Carlo sample path, solves anoptimal power flow (OPF) problem, which is computationallyexpensive. Monte Carlo technique was also used in [9] wherea reduction of the random variable dimension is made usinga nonhomogeneous Markov chain model based on a partitionof the system state space.A particularly relevant prior work is [10] where the authorsconsider the problem of LMP/congestion forecasting fromthe vantage point of an external observer who has access topublicly available historical data only. Our work, in contrast,considers the forecasting problem from the vantage point ofa system operator who has access to the system operatingcondition at the time of forecasting. In terms of forecastingmethodology, the main difference between our approach andthat in [10] lies in the different uses of conditioning in evaluat-ing the conditional probability distribution of LMP/congestion.The authors of [10] introduce and exploit the decompositionof a multi-dimensional load space into critical regions (calledsystem pattern regions) that are estimated using historicaldata . The work in [10] aims to address the following issue:Given a possible future point L in a multi-dimensional loadspace, what is the probability distribution of the estimatedcritical regions that contain L ? Since each critical regioncorresponds to a specific LMP/congestion, the technique in[10] gives a heuristic estimate of the probability distribution The estimated critical regions are therefore random quantities. of LMP/congestion by conditioning on load L at some futurepoint in time.In contrast to [10], our objective is to forecast directly theprobability distribution of future LMP/congestion, conditionalon the current system operating point . Because a systemoperator has access to all private and public information aboutsystem conditions, the critical regions are computed exactlyvia a multiparametric program. This allows us to incorporateload and generation forecasts and obtain the (conditional)probability distribution of future LMP/congestion directly.Several techniques have been proposed to approximate theLMP distribution at a future time. In [11], a probabilisticLMP forecasting approach was proposed based on attachinga Gaussian distribution to a point estimate. The advantage isthat the technique can incorporate various point forecastingmethods. The disadvantage, on the other hand, is that networkeffects are not easy to incorporate. The authors of [12] and[13] approximate the probabilistic distribution of LMP usinghigher order moments and cumulants. These methods arebased on representing the probability distribution as an infiniteseries involving moments or cumulants. In practice, computingor estimating higher order moments and cumulants are verydifficult; lower order approximations are necessary. C. Organization and Notations
This paper is organized as follows. Section II introducesa model of the real-time economic dispatch and the ex-anteLMP computation. The key (known) results in multiparametricprogramming that form the basis of the proposed probabilisticforecasting approach are discussed in Section III. Details ofthe proposed techniques are given in Section IV and SectionV. Numerical results are presented in Section VI and followedby some concluding remarks in Section VII.The notations used in this paper are standard. For a randomvariable x , its expected value is denoted by ¯ x . We use ˆ θ todenote an estimate of parameter θ . For a random process y t ,the forecast of y t + T using all the information available at time t is denoted by y t + T | t , where T is the prediction horizon. Weuse ¯ y t + T | t to denote the conditional mean of y t + T given y t (and possibly additional information at time t ). The notation z ∼ N ( µ, Σ) means that z is a Gaussian random vector withmean µ and covariance matrix Σ . The indicator function I A ( x ) is one if x ∈ A and zero otherwise.Vectors and matrices are in the real field. I n denotes an n × n identity matrix and an n dimensional column vectorwith all ones, where n is the number of buses in the system. denotes a matrix/vector of zeros with different but compatibledimensions. We use superscript “ ⊺ ” to denote the transpose ofa matrix and the complement of a set I is denoted by I c .II. R EAL -T IME E CONOMIC D ISPATCH AND E X -A NTE
LMPWe consider an ex-ante LMP model that arises from areal-time economic dispatch. Our model is adopted from thestylized model in [14] that captures the process of com-puting LMP by independent system operators. Specifically,the system operator solves a DC-OPF problem for optimaleconomic generation adjustment that meets load and stochastic generation forecasts for the next dispatch interval and satisfiesgeneration and transmission constraints.We describe here a simplified real-time economic dispatchformulation for the ex-ante LMP model. For simplicity, weassume that each bus has a traditional generator, a stochasticgenerating unit and a load. The DC-OPF problem at time t isgiven by min g c ⊺ g subject to ⊺ ( g + w t +1 | t − d t +1 | t ) = 0 ( λ t +1 ) L − ≤ S ( g + w t +1 | t − d t +1 | t ) ≤ L + ( µ + t +1 µ − t +1 ) G − ≤ g ≤ G + (1)where c vector of real-time generation offers; g vector of ex-ante dispatch at time t + 1 ; d t +1 | t vector of 1-step ahead load forecast at time t ; w t +1 | t vector of 1-step ahead forecast of stochastic gen-eration at time t ; S shift factor matrix; G + /G − vector of max/min generation capacities; L + /L − vector of max/min transmission limits; λ t +1 shadow price for the energy balance constraint attime t + 1 ; µ + t +1 /µ − t +1 shadow prices associated with max/min transmis-sion constraints at time t + 1 .The stochastic generation referenced above can be of anyform of renewable integrations including renewable energygeneration, distributed generation, and demand response. Herewe assume that such stochastic generation is non-dispatchablewith possible curtailment. The proposed forecasting methodapplies to other cost functions such as a piece-wise affinefunction or a quadratic function. This will be further addressedin the following sectionsThe LMP value at each bus is the sum of the marginal priceof generation at the reference bus and the marginal congestionprice at the location associated with the active transmissionconstraints. By the Envelope Theorem, the ex-ante LMP π t +1 at time t + 1 is given by π t +1 = λ t +1 − S ⊺ µ + t +1 + S ⊺ µ − t +1 , (2)where λ t +1 , µ + t +1 and µ − t +1 are from (1) at time t .From (2), we note that the LMP is determined by themarginal generator through λ and the congestion patternthrough µ + and µ − . By a congestion pattern we mean the setof transmission lines where power flows have reached theirlimits. The forecasting of congestion pattern is a sub-problemof the forecasting of LMP.III. M ULTIPARAMETRIC P ROGRAMMING
A multiparametric program (MPP) is a mathematical pro-gram indexed by a vector of parameters. The multiparametricprogramming problem is to solve the mathematical programfor all values of the parameter vector (or a parameter space By ex-ante LMP we mean that the computation of real-time dispatch andassociated prices at time t + 1 is based on the estimated system operatingpoint at time t and 1-step ahead load and generation forecasts. of interest). Consequently, an MPP captures simultaneousvariations of distributed parameters in a network, which mo-tivates its application in the problem of LMP forecast. Herewe summarize some of the key theoretical results essentialto the development of the proposed probabilistic forecastmethod. The definitions and results follow [15]. In particular,we are interested in the linear and quadratic programs forwhich the multiparametric programming problems are calledmultiparametric linear programs (MPLP) and multiparametricquadratic programs (MPQP) respectively.Consider a general right-hand side MPP as follows: min x z ( x ) subject to Ax ≤ b + Eθ ( y ) (3)where x is the decision vector, θ the parameter vector, z ( · ) thecost function , y the Lagrangian multiplier vector, and A , E , b are coefficient matrix/vector with compatible dimensions.An MPP solver finds the subset Θ of a given polyhedronset such that the mathematical program (3) is feasible anddetermines the expression x ∗ ( θ ) ∈ X ∗ ( θ ) of the optimizer (orone of the optimizers if it is not unique), where x ∗ ( θ ) is theoptimal solution and X ∗ ( θ ) the set of optimal solutions of (3)given θ .The multiparametric programming analysis builds on theconcepts of optimal partition and critical region. Definition 1 ([15]) . Let J denotes the set of constraint indicesin (3). For any subset I ⊆ J , let A I and E I be the corre-sponding submatrices of A and E , respectively, consisting ofrows indexed by I . An optimal partition of the index set J associated with parameter θ ∈ Θ is the partition ( I ( θ ) , I c ( θ )) where I ( θ ) , { i ∈ J | A i x ∗ ( θ ) = b + E i θ, for all x ∗ ( θ ) ∈ X ∗ ( θ ) } , I c ( θ ) , { i ∈ J | A i x ∗ ( θ ) < b + E i θ, for some x ∗ ( θ ) ∈ X ∗ ( θ ) } . For a given θ ∈ Θ , let ( I , I c ) , ( I ( θ ) , I c ( θ )) . The criticalregion related to the index set I is defined as Θ I , { θ ∈ Θ | I ( θ ) = I } , which is the set of all parameters θ ∈ Θ with the same activeconstraint set I at the optimum(s) of problem (3).By Definition 1, the optimal partition specifies two sets ofconstraints by which the congestion pattern is determined: oneset is a combination of active constraints at the optimum(s)and the other inactive. A critical region is the set of allparameter vectors for a certain optimal partition. The structureof critical region is a key property for the development of theonline simulation technique, so we present it in Section V.Here, we summarize the global property of critical region andits connection with the solution structure of MPP (3). Theconnection between the feasible parameter space and criticalregions is summarized in the lemma below. Lemma 1 ([15]) . The feasible parameter space Θ can bepartitioned into critical regions { Θ i } . If there is no (primalor dual) degeneracy, such a partition is unique. By right-hand side we mean the parameter vector θ is on the right-handside of the constraint inequalities. In this paper, we only consider the cost functions that are linear, piece-wiseaffine or quadratic.
Assume that there is no (primal or dual) degeneracy. In thelinear (or piece-wise affine) cost case, the optimizer is uniqueand the associated Lagrangian multiplier vectors are constant,for all parameters within each critical region. Therefore, eachcritical region is associated with a unique LMP vector. In thequadratic cost case, the LMP vector is a uniquely defined affinefunction of the parameter within each critical region, whichfollows the theorem below.
Theorem 1 ([16]) . If the MPQP problem (3) is not degen-erate , the Lagrangian multiplier vector associated with theoptimal solution x ∗ ( θ ) is a uniquely defined affine function of θ within each critical region. Note that degeneracy may happen for sufficiently largesystems or high dimensional parameters. Although the La-grangian multiplier vector may not be unique in such cases,a consistent tie breaking rule can be introduced to guaranteethe uniqueness.A key property of the optimizer for both MPLP and MPQPis as follows.
Theorem 2 ([15] [16]) . If the optimizer x ∗ ( θ ) of MPLP/MPQP(3) is unique, then it is continuous and piece-wise affine over Θ . In particular, x ∗ ( θ ) is an affine function of θ in each criticalregion Θ i . If there exists multiple optimizers, it is alwayspossible to define a continuous and piece-wise affine function x ∗ ( θ ) ∈ X ∗ ( θ ) for all θ ∈ Θ . This theorem is crucial to the development of the proposedforecasting technique. In essence, we treat the randomness asa vector of parameter and solve the MPPs ahead of time. Theonline computation of the optimal solution is highly simplifiedbecause the piece-wise affine mapping between the optimalsolution and the parameter vector has already been obtainedoffline.The problem of solving MPPs, including the computation ofcritical regions and corresponding piece-wise affine functions,has been well studied. The first method for solving right-hand side MPLPs was proposed by Gal and Nedoma [17]. Ageometric approach for critical region partition was describedin [15]. The piece-wise affine relationship of the parametervector and the optimal solution was proved for MPLP andMPQP in [17] and [16], respectively. In this paper, all cal-culations related to MPPs are performed using MPT3 toolbox[18] except for dynamical generations of critical regions in theproposed online simulation technique presented in Section V.IV. P
ROBABILISTIC F ORECASTING A LGORITHM
Probabilistic forecasting , in contrast to point forecasting,aims to provide the probability distribution of a future LMP. Inparticular, given the estimated system operating point at time t and load and generation forecasts, the probabilistic forecastat time t of the LMP at time t + T is given by the conditionalprobability distribution f t + T | t of the LMP vector. Entries ofthe LMP vector are LMPs at individual buses in the system.Since each congestion pattern can be mapped from an LMP Note that the cost function of the MPQP problem is assumed to be strictlyconvex. Θ Θ Θ Θ Θ Θ Θ Feasible ParameterSpace Θ θ t + T θ t tt + 1 f t + T | t f t Fig. 1. Geometric intuition of the proposed algorithm. vector, we only discuss the probabilistic forecasting of LMPhere; the probability distribution of congestion can be obtainedfrom that of LMP.The key to probabilistic LMP forecasting is to capturespatial and temporal dependencies. Spatial correlations amongLMPs arise naturally from the optimization that governs thereal-time dispatch. Temporal correlations, on the other hand,are the results of time dependencies in load/generation fore-casts. The system randomness may also include occurrencesof random contingencies. In this section, we first give anoverview of the proposed forecasting algorithm. Details onaddressing these dependencies are then discussed.
A. Overview
The basic idea of the proposed probabilistic forecastingtechnique is using multiparametric programming analysis tocharacterize the variation of the real-time LMP with respectto the random load and generation. By formulating the DC-OPF problem (1) as an MPP in the form of (3), the real-timeLMP can be expressed as a function of the random load andgeneration. The distribution of the LMP vector at a futuretime can be obtained from the probabilistic forecasts of thestochastic load and generation.As illustrated in Fig. 1, load and stochastic generation fore-casts in (1) are treated as parameters, denoted by θ = ( d, − w ) ,where d is the vector of parametric load and w the vector ofparametric generation (which is treated as a negative load).The feasible parameter space Θ is partitioned into criticalregions { Θ , · · · , Θ } . Within each region Θ i , the optimaldispatch is associated with the same Lagrange multipliers, andhence a unique LMP vector π i for all θ ∈ Θ i . Given thenetwork parameters, the MPLP solver computes the partition { Θ i } . Correspondences { Θ i , π i } are then obtained by theLagrange multipliers and the LMP model (2). Note that thiscomputation does not depend on the actual realization of loadand generation. Therefore, the computation of the partition andthe correspondences can be obtained offline.Consider now the trajectory of a realization of the randomload and generation process θ t as illustrated in Fig. 1. Giventhe realization θ t and system measurements at time t , we are interested in the conditional probability distribution f t + T | t ( i ) = P [ θ t + T ∈ Θ i | θ t ] . (4)As depicted by the shaded circles in Fig. 1, uncertaintiesassociated with load and generation forecasts increase withtime. At time t , the realization of parameter θ t is knownand thus the distribution f t is an unit vector. But parameter θ t + T may take values from several critical regions where LMPvalues and congestion patterns are different. Therefore, theforecast probability mass function f t + T | t of θ t + T | t may haveseveral non-zero elements.The proposed probabilistic forecasting algorithm involvestwo parts explained in the following subsections: the com-putation of critical regions and the estimation of conditionalprobability distributions, where the former is computed offlineand the latter online. B. Forecast with Varying Operational Conditions
We describe here a baseline formulation from which criticalregions are obtained. We formulate the DC-OPF (1) used tocompute LMP at time t + T as the following right-hand sideMPLP with the uncertainty parameter θ consisting of onlystochastic load and generation. min g c ⊺ g subject to ⊺ − ⊺ S t + T − − S t + T − I n − I n g ≤ ⊺ ⊺ − ⊺ − ⊺ S t + T − S t + T − − S t + T − − S t + T − θ + L + t + T − − L − t + T − G + t + T − − G − t + T − (5)where S t + T − , L + t + T − , L − t + T − , G + t + T − , and G − t + T − areshift factor matrix, vectors of max and min transmission limits,and vectors of max and min generation capacities at time t + T − , respectively.Here we allow time varying but known system parameters such as shift factors, flow limits on transmission lines, andlimits on generations, restricting uncertainties only to loadand generation. The idea is to include deterministically sched-uled events in the forecasting problem. Examples include thescheduled changes in network topology [19], and generationcapacity and transmission limit [20]. C. Forecast in the Presence of Probabilistic Contingencies
The baseline MPLP formulation described in Section IV-Bcan be extended to include the presence of probabilistic con-tingencies for unexpected events. For example, a transmissionline may be tripped in a storm or generation capacity reduceddue to faults. Such uncertainties in the system parametersneed to be handled differently from those associated with thestochastic load and generation. This is the case for linear or piece-wise affine cost functions. If the costfunction in the real-time economic dispatch is quadratic, the forecast distribu-tion f t + T | t of LMP π t + T can be obtained from the (continuous) distributionof θ t + T , and the uniquely defined affine function of the Lagrangian multipliervector within each critical region. Because unexpected changes of system configurations aretypically small probability events, we assume that there are atotal of K possible system configurations at time t + T . Fromhistorical data, we assume that system configuration k happenswith an estimated probability ˆ p k .We solve the baseline MPLP (5) for each system config-uration and obtain critical regions for system configuration k denoted by { Θ ( k ) i } . By the total probability theorem, theprobabilistic forecast of LMP at time t + T is therefore givenby f t + T | t = K X k =1 ˆ p k f ( k ) t + T | t , (6)where f ( k ) t + T | t is the forecast distribution under system config-uration k using critical regions { Θ ( k ) i } for k = 1 , , · · · , K .We illustrate the above idea with an example. Consider thecase when at most one of K contingencies can occur betweentime t and t + T . We model the probabilistic contingencyas independent tosses of a K + 1 faced dice where con-tingency k occurs with probability p k and no contingencywith p = 1 − P Kk =1 p k . We further assume that, once aparticular contingency occurs, it remains until time t + T . Wethen solve each MPLP problem (5) for all K + 1 possiblesystem configurations, including the normal condition. Theprobabilistic forecast of LMP at time t + T is given by f t + T | t = ( f ( k ) t + T | t if contingency k occurs P Kk =0 p k f ( k ) t + T | t otherwise (7)where f ( k ) t + T | t is the forecast distribution under system config-uration k with the feasible space Θ ( k ) t + T | t for k = 0 , , · · · , K .Note that system configuration denotes the normal condition.Having obtained the critical regions, we now consider theproblem of computing the conditional distribution of LMP attime t + T , given the load and generation forecast at time t . D. Probability Distribution Estimation
The estimation of conditional probabilities in (4) dependson statistical models of load and generation. Such models canbe obtained from either models of the load and stochasticgeneration process or specific prediction methods used togenerate load and stochastic generation forecasts.As illustrations, we present a directional Gaussian randomwalk model and an autoregressive (AR) noise model for therandom load and generation processes here. It should be notedthat any statistical model or prediction method can be applied.The purpose of using these models is to gain insights into thebehavior of forecasting performance by taking advantage ofsome of analytically tractable properties.
1) A Directional Random Walk Model:
We first considera random walk model of the load/stochastic generation basedon a given mean trajectory. Such a model represents a case ofminimally informative forecast. Note that the mean trajectorycan be any available forecast. For example, a reasonable meantrajectory is the day-ahead load forecast.Assume that load/stochastic generation θ t follows a randomwalk process with a (known) mean trajectory ¯ θ t : θ t = θ t − + ¯ θ t − ¯ θ t − + ǫ t , (8) where ǫ t ∼ N ( , Σ) . Given the realization θ t at time t , theactual load/generation at time t + T is given by θ t + T = θ t + ¯ θ t + T − ¯ θ t + t + T X i = t +1 ǫ i . (9)Therefore, the distribution of θ t + T conditioning on θ t is: θ t + T ∼ N (¯ θ t + T | t , Σ T ) (10)where ¯ θ t + T | t = θ t + ¯ θ t + T − ¯ θ t is the conditional mean of θ t + T and Σ T = T Σ is the cumulative variance within predictionhorizon T .
2) An AR Noise Model:
The second model we consideris an AR(1) noise model where we assume the deviation ofthe load or generation from the expected value is an AR(1)process. This is a case when the load or stochastic generationis highly structured. In particular, θ t = ¯ θ t + a t , a t = φa t − + ǫ t , (11)where ¯ θ t is the (known) mean trajectory, φ the parameter ofthe AR process, and ǫ t ∼ N ( , Σ) . Given the realization θ t =¯ θ t + a t at time t , the noise at time t + T is given by a t + T = φ T ( θ t − ¯ θ t ) + T − X i =0 φ i ǫ t + T − i , (12)and the actual load/generation at time t + T is given by θ t + T = ¯ θ t + T + φ T ( θ t − ¯ θ t ) + T − X i =0 φ i ǫ t + T − i . (13)Therefore, the distribution of θ t + T conditioning on θ t is: θ t + T ∼ N (¯ θ t + T | t , Σ T ) (14)where ¯ θ t + T | t = ¯ θ t + T + φ T ( θ t − ¯ θ t ) is the conditional mean of θ t + T , and Σ T = P T − i =0 φ i Σ the cumulative variance withinprediction horizon T .To sum up, for both models, the conditional probability of θ t + T falling in critical region Θ i given θ t is: f t + T | t ( i ) = Z Θ i exp {− ( x − ¯ θ t + T | t ) ⊺ Σ − T ( x − ¯ θ t + T | t ) } p (2 π ) n | Σ T | dx, (15)where ¯ θ t + T | t and Σ T are model associated statistics given in(10) and (14).In general, Monte Carlo techniques are necessary to estimatethe conditional probability (15). To accelerate the samplingprocess, importance sampling technique is used. In particular,for each critical region Θ i , instead of drawing values from dis-tribution N (¯ θ t + T | t , Σ T ) , we use N (¯ v (Θ i ) , Σ T ) where ¯ v (Θ i ) is the mean of all vertices of critical region Θ i . The estimateof f t + T | t ( i ) is then given by ˆ f t + T | t ( i ) = 1 N N X j =1 I Θ i ( s j ) g ( s j ) h ( s j ) , (16)where samples { s , · · · , s N } are drawn from N (¯ v (Θ i ) , Σ T ) ,and g ( · ) and h ( · ) are probability density functions (PDFs) ofdistribution N (¯ θ t + T | t , Σ T ) and N (¯ v (Θ i ) , Σ T ) , respectively.Note that the importance distribution N (¯ v (Θ i ) , Σ T ) onlyshifts the mean of the nominal distribution N (¯ θ t + T | t , Σ T ) butkeeps the variance the same.
3) A Quadratic Cost Case:
If the cost function in the DC-OPF (1) is quadratic, the distribution of the LMP π t + T cannotbe estimated by the conditional probabilities of θ t + T fallingin each critical region. Because the LMP in this case is anaffine function of the parameter vector by Theorem 1.Here we derive the conditional distribution f t + T | t at time t of the future LMP π t + T at time t + T given the conditionaldistribution of θ t + T . By Theorem 1 and the LMP formulation(2), for each critical region Θ i , there exists an affine function π i ( · ) : Θ i → Π i such that π i ( θ ) = U i θ + v i , for all θ ∈ Θ i , where Π i is the codomain and U i and v i are associatedcoefficient matrix/vector. Note that these coefficients can beobtained from the affine function of the Lagrangian multipliervector and the LMP formulation (2).Given the conditional distribution N (¯ θ t + T | t , Σ T ) of θ t + T ,the conditional PDF of π t + T is given by f t + T | t ( π ) = N X i =1 I Π i ( π ) f t + T | t, Π i ( π ) , (17)where N is the number of critical regions, and f t + T | t, Π i ( · ) thePDF of N ( U i ¯ θ t + T | t + v i , U ⊺ i Σ T U i ) , which is the conditionaldistribution of π t + T in codomain Π i .In summary, the proposed algorithm treats load/generationas parameters, formulates the DC-OPF (1) as an MPP (3),determines the critical regions, and computes the conditionaldistribution of LMP and congestion using load/generationforecasts and the real-time operation conditions. These systemconditions, such as transmission rate, generation capacity, andnetwork topology, are allowed to vary with time but known tothe operator at the time of forecast. Contingency models arealso incorporated in the proposed technique.V. O NLINE F ORECASTING VIA D YNAMIC C RITICAL R EGION G ENERATION
A limiting factor in the proposed technique above is thecomputational cost associated with the multiparametric pro-gramming. Since the solution structure is characterized bycritical regions, all critical regions that partition the parameterspace have to be calculated. Although such computation canbe made offline, it may not be computationally tractable forlarge systems because the number of critical regions may growexponentially with the number of constraints.In this section, we propose an online Monte Carlo technique,referred to as dynamic critical region generation (DCRG),that significantly reduces the computation cost. The idea isto take advantage of the fact that, in practice, random loadand generation processes are bounded and tend to concentratearound their mean trajectories. As a result, a small fractionof critical regions represents the overwhelming majority ofobserved critical regions. When a parameter falls in a criticalregion that was visited before, the Lagrange multipliers thatare used to generate LMPs can be obtained directly from theaffine mappings associated with that critical region, withoutsolving the DC-OPF (1).
The key idea of DCRG, therefore, is to compute on demandthe critical region and the associated coefficients of the affinemapping of parameter to the LMP. This computation, fortu-nately, is no more than elementary matrix inversions and mul-tiplications. The computation procedure is given by the follow-ing theorem that summarizes known results in MPLP/MPQP[15] [16]. The proof is provided in the Appendix.
Theorem 3.
Given a parameter θ , let ( I , I c ) be the optimalpartition of the index set J of (3) associated with θ . Let A I , E I and b I be, respectively, the submatrices of A , E and subvector of b corresponding to the index set I . Let A I c , E I c and b I c be similarly defined for the index set I c .Assume that MPP (3) is neither primal nor dual degeneratefor all θ . Denote the critical region that contains θ by Θ I .(1) For the MPLP (3) with cost function z ( x ) = c ⊺ x , thecritical region Θ I is given by Θ I = (cid:8) θ (cid:12)(cid:12) A I c A − I ( b I + E I θ ) < b I c + E I c θ (cid:9) . (18) (2) For the MPQP (3) with cost function z ( x ) = x ⊺ Hx where H is positive definite, the critical region Θ I isgiven by Θ I = { θ | θ ∈ P p \ P d } (19) where P p and P d are polyhedra defined by P p = (cid:26) θ (cid:12)(cid:12)(cid:12)(cid:12) A I c H − A ⊺ I ( A I H − A ⊺ I ) − ( b I + E I θ ) < b I c + E I c θ (cid:27) P d = { θ | ( A I H − A ⊺ I ) − ( b I + E I θ ) ≤ } . In applying DCRG to LMP forecasting using online MonteCarlo simulation, we generate samples of random load orgeneration, either based on load/generation models or fromhistorical data. Instead of directly simulating the real-timemarket operation as in [8], we check if the generated samplefalls into a critical region that has been used before. If itdoes, the LMP can be generated using the affine mapping ofLagrangian multipliers directly without solving the DC-OPF(1). If the parameter does not belong to any critical regionin the database, a DC-OPF is solved and a critical regioncontaining the parameter is computed according to (18) forlinear cost case or (19) for quadratic cost case in Theorem 3.VI. E
VALUATION
In this section, we present simulation results to compare per-formances of the proposed probabilistic forecasting algorithmwith some benchmark methods. We first show results of a 3bus system with a linear cost function to gain insights into thebehavior of the proposed algorithm under various scenarios.Simulations using the IEEE 118 bus system with a quadraticcost function are then presented to demonstrate the scalabilityof the proposed algorithm and the effectiveness of the onlineheuristic approach given in Section V.
A. Benchmarks and Performance Measure
We compared the proposed techniques with some existingbenchmarks for forecasting and computation performance.Since, to our best knowledge, there is no probabilistic fore-casting techniques for the operator in the literature, we used
21 3 c $/MWh c $/MWh G +1 = 130 MW G − MW G +2 = 200 MW G − MW Fig. 2. 3 bus system. the direct Monte Carlo simulation method proposed in [8]as the probabilistic forecasting benchmark . We also includedcomparisons of the proposed technique with two well knownpoint forecasting methods to illustrate the performance gain. Inparticular, the deterministic prediction uses the mean trajectoryof load/stochastic generation ¯ θ t + T to calculate LMP andcongestion pattern at time t + T . The certainty equivalenceprediction incorporates measurements at time t and uses theconditional mean trajectory ¯ θ t + T | t .Before presenting numerical results, we introduce a per-formance evaluation metric of probabilistic forecasts. TheLMP /congestion pattern is a discrete random vector. Theprobabilistic forecast of such a random quantity belongs to theso-called categorical forecast, and its performance is measuredby the consistency as well as the statistical concentration ofthe forecast. A standard metric [21] for this type of forecast isthe Brier Score (BS) [22] that measures the average distance(2-norm) between the forecast distribution f t + T | t and thepoint mass distribution at the realized random variable π t + T .Specifically,BS ( f t + T | t ) = E k f t + T | t − δ ( π t + T ) k , (20)where the expectation is taken over all randomness betweentime t and t + T . In (20), f t + T | t is the conditional probabilityvector whose i th entry is given by f t + T | t ( i ) = Pr( π t + T = π i ) ,and δ ( x ) is the unit vector that is one at entry x and zeroelsewhere. This score ranges from for a perfect forecast to for the worst possible forecast.Since BS is a succinct formula to measure the overall per-formance in terms of uncertainty, reliability and resolution, wealso provide a more intuitive assessment — reliability diagram.Reliability diagram is a graph of the observed frequency of anevent plotted against the forecast probability of an event. Itmeasures how closely the forecast probabilities of an eventcorrespond to the actual chance of observing the event. B. Case Study: A 3 Bus System
Consider a 3-bus system as depicted in Fig. 2. Generatorincremental costs and capacity limits are presented in thefigure. All lines are identical with the maximum capacity of
MW. We want to point out that the direct Monte Carlo simulation approachgenerates exactly the same probabilistic forecast as the proposed technique. Note that the discreteness of LMP is only for linear or piece-wise affinecost functions.
TABLE IC
RITICAL REGIONS , LMP
S AND CONGESTION PATTERNS FOR THEBASELINE .Critical Region LMP Congestion1 (0, 130) (10, 10, 10) (0, 0, 0)2 (130, 170) (15, 15, 15) (0, 0, 0)3 (170, 200) (10, 20, 15) (1, 0, 0) Time Index(5-min Interval) B r i e r S c o r e Alg-DAlg-CAlg-P (a) Random walk model
Time Index(5-min Interval) B r i e r S c o r e Alg-DAlg-CAlg-P (b) AR(1) noise model
Fig. 3. Impact of load statistical models.
1) Baseline:
We first evaluated the baseline algorithm withthe two load models described in Section IV-D. Note thatonly the load at bus 2 was stochastic, which gave the onedimensional parametric linear program.Table 1 shows the critical regions of the parametric linearprogram, and the associated LMPs and congestion patterns. Inthis case, the parameter load d at bus 2 was partitioned intothree segments D = { (0 , , (130 , , (170 , } .We used a straight line ranging from MW to
MW with MW increments as the mean trajectory of load ¯ d t . The coefficient φ in AR(1) noise model was set at . .The independent noise sequence ǫ t in both models followedthe standard normal distribution, i.e. , N (0 , . Monte Carlosimulations were used to obtain estimated BSs in Fig. 3.From Fig. 3, we observed that the proposed probabilisticforecasting algorithm “Alg-P” consistently outperforms thedeterministic “Alg-D” and certainty equivalence “Alg-C” pre-dictors in both load models. The superior performance of theproposed technique in these two extreme models (a minimallyinformative model and a highly structured model) shows itscapability of incorporating different load forecasting methodsand its forecasting power. Two interesting phenomena areworth closer examinations. First, for both load models, peaksoccurred at the boundaries of two neighboring critical regionswhen the mean load ¯ d t is at t = 10 and at t = 30 . Atthe boundary point, the probability of d t + T falling in eitherthe left or the right critical region was the same. Roughly half In Table I, the triple for LMP contains price values at bus 1-3. Forcongestion, the triple summarizes status of line 1-2, line 1-3, and line 2-3, respectively, where “ ” indicates positive congestion (the flow reaches theline limit in the positive direction), “ − ” negative congestion, and “ ” nocongestion. Forecast probability O b s e r v ed f r equen cy σ =1.5 σ perfect forecast σ =2 σ Fig. 4. Reliability diagrams for critical region 1 with different load forecastuncertainties. of the time Alg-D predicted the LMP correctly, the other halfwas completely wrong. From the definition in (20), the BS ofAlg-D should be .Second, a more subtle point, the scores for the random walkmodel showed a slight asymmetry with respect to boundariesof neighboring critical regions: the BSs of Alg-P and Alg-Cat the second peak were lower than that of the first peak, andthe ranges of non-zero score neighborhood of the second peak(from time to ) of all three algorithms were wider thanthat of the first peak (from time to ). This phenomenonarose primarily from the process of generating the sampletrajectory of load. Since the entire sample trajectory wasgenerated at once, the deviation | d t − ¯ d t | from the mean growsover time which leads to a bigger variance of the time crossingthe second boundary than that of the time crossing the firstboundary . In other words, comparing to the probability ofcrossing the first boundary at t = 10 , the probability ofcrossing the second boundary at t = 30 is lower, but theprobability in its neighborhood is higher. Therefore, the scoresof Alg-P and Alg-C at the second peak were lower, and theranges of non-zero score neighborhood bigger.To evaluate the robustness and reliability of the proposed al-gorithm, we tested different levels of load forecast uncertainty.Specifically, we varied the variance of the noise in the randomwalk load model: by default, ǫ ∼ N (0 , σ ) , where σ = 1 , andthen we used the distribution N (0 , σ ) with different valuesof σ . The reliability diagrams for the probabilistic forecastof critical region 1 using different load forecast uncertaintiesare given in Fig. 4. The results show good resolution at theexpense of reliability.
2) Forecast with probabilistic generation outage:
We nowconsidered the case of a generating unit outage with a partialloss of capacity, assuming that the maximum capacity of thegenerator at bus 1 can be reduced to
MW with probability p . Other settings were the same as those in the first scenario.The critical regions under this configuration became D out = { (0 , , (100 , } . To predict a future price, we consideredall critical regions { D , D out } that load d t + T may fall in,where D refers to the critical regions in Table I under normalconditions. For the outage frequency p , we chose two levels: p = 0 . and p = 0 . . The random walk model was adoptedto generate load profiles. As benchmarks, both deterministic Time Index(5-min Interval) B r i e r S c o r e Alg-D p=0.01Alg-C p=0.01Alg-P p=0.01Alg-D p=0.001Alg-C p=0.001Alg-P p=0.001
Fig. 5. Impact of outage frequency p . and certainty equivalence forecasting algorithms also tookcontingencies into consideration.From Fig. 5, the behaviors of all three algorithms with thesame outage frequency p were similar to that in the first sce-nario. Boundary effects and peak asymmetries were observedwith contingencies as well. Compare the performances of eachalgorithm with different outage frequencies: the smaller theoutage rate the better the forecast. C. Case Study: IEEE 118 Bus System
The IEEE 118 bus system was used to show the scalabilityof the proposed technique and the effectiveness of the heuristicalgorithm described in Section V. The simulations results werefocused on the complexity comparison between the proposedtechniques and the probabilistic forecast benchmark.We introduced 12 wind generators (roughly of buses)with ∼ renewable penetration . The wind generatorswere located at bus 25, 26, 90, 91, 100, 103, 104, 105, 107,110, 111, and 112. The selection of these locations wereintended to represent two wind farms, the small one has 2wind generators, and the large one has 10 wind generatorsconcentrated on a few neighboring buses. All wind generatorswere assumed to be identical with the maximum capacity of MW. Denote the wind generation space by the hypercube W ∈ R . We imposed the maximum capacity of MW ontransmission line 8, 126, and 155. The load profile, generatorcapacities and cost functions, and line and bus labels werereferred to as in “case118” in MATPOWER [23]. Note thatthe cost function in this system was quadratic, thus the LMPwas an affine function of the wind production within eachcritical region.The mean trajectory ¯ w ( i ) of each wind generator i wasassumed to be linear, for i = 1 , · · · , . In particular, thetrajectory starts from penetration level, i.e. , ¯ w ( i ) =35 . , at time , and ends at penetration level, i.e. , ¯ w ( i ) = 70 . , at time . The increment was assumed tobe constant, i.e. , . MW. LMP values at and penetration levels are given in Fig. 6. We observed that higherrenewable penetration reduced LMP at most buses, but raisedthe LMPs at bus 93, 94, 95, and 96. The reason of suchnonintuitive increase was the congestion of transmission line By x % renewable penetration we mean the mean value of total windgeneration is x % of the total electricity load ( MW in this system).Note that load was assumed to be deterministic in this case.
Bus Index L M P
10% Penetration20% Penetration
Fig. 6. LMPs at and renewable penetration levels.
33 34 35 36 37 38 39 40 41
LMP@Bus49 P D F Predicted DistributionGaussian Approximation
10 15 20 25 30 35 40
LMP@Bus90 P D F Predicted DistributionGaussian Approximation -1000 0 1000 2000 3000 4000 5000 6000 7000
LMP@Bus94 P D F Predicted DistributionGaussian Approximation -12000 -10000 -8000 -6000 -4000 -2000 0
LMP@Bus100 P D F × -3 Predicted DistributionGaussian Approximation
Fig. 7. Predicted LMP distributions and Gaussian approximations at bus 49,90, 94, and 100.
155 caused by the increased wind production from the bigwind farm. The random walk model was used for the stochasticgeneration profile. The distribution of the independent noiseprocess ǫ t was set to be the standard multivariate Gaussian.Results for the 10-step ahead prediction at time t = 0 areprovided in Fig. 7. The predicted marginal distribution of LMPat bus 49 exhibits the a Gaussian distribution as it is wellfitted by the Gaussian approximation with the sample meanand sample variance as distribution parameters. However, suchGaussian characteristics were not observed on the other threebuses. According to the distributions of LMP at bus 94 and100, the extreme values of LMP occasionally appeared asspikes, which were caused by the network congestions.In the following, we will focus on the complexity perfor-mance of the proposed forecasting techniques. Theoretically, Critical Region Index P r obab ili t y Fig. 8. Conditional distribution of observed critical regions at time . Number of Monte Carlo Simulations E x pe c t ed N u m be r o f DC - O P F s Alg-MCAlg-DCRG
Fig. 9. The expected number of DC-OPF computations vs. the total numberof Monte Carlo simulations. there were at most critical regions, because the constraintsin the DC-OPF (1) for this example included 1 energy balanceconstraint, 6 transmission constraints, and 108 generator ca-pacity constraints. For the given hypercube W of the parameterspace, there were in total critical regions. But for thegenerated 10,000 samples at time , only critical regionswere observed and more than samples fell in criticalregions, as shown in the distribution over the observed critical regions in Fig. 8.Instead of exploring the entire wind production space of-fline, we implemented the online forecasting algorithm DCRGgiven in Section V. Fig. 9 shows the comparison of computa-tional cost between the proposed DCRG algorithm, labeledby “Alg-DCRG”, and the direct Monte Carlo simulations,labeled by “ Alg-MC”. Both algorithms present approximatelylinear growth in the logarithm scales. But the proposed DCRGalgorithm provided more than three orders of magnitudereduction in the number of DC-OPF computations requiredin the simulation.Finally, we provide the computation time comparison forthe three probabilistic forecasting techniques in Table II.All computational times were evaluated by implementing thealgorithms in Matlab environment with the default “quadprog”solver and an external MPT3 [18] toolbox on a desktopwith an Intel Core i7-3770 CPU at 3.4 GHz and 8 GBmemory. No attempts were made to optimize the efficiencyof the algorithms and their simulations. From Table II, wecan conclude that the direct Monte Carlo simulation approachAlg-MC failed to meet the time constraint for real-time LMPforecasting, since 10,000 samples took more than 14 hoursto generate the distribution. The proposed techniques, on theother hand, only took less than half a minute in the online TABLE IIC
OMPUTATION TIME ( IN SECONDS ) FOR
SAMPLES .Offline Online TotalAlg-P 160.60 23.32 183.92Alg-DCRG — 23.59 23.59Alg-MC — 52022.57 52022.57 computation, demonstrating the efficiency for the real-timeLMP forecasting. VII. C
ONCLUSION
This paper presents a new methodology for the short-termforecasting of real-time LMP and congestion for system op-erators. Based on a multiparametric programming formulationof DC-OPF, we have developed an approach that exploits theparametric structure of DC-OPF solutions to obtain conditionaldistributions of future LMP and congestion.For system operators, the proposed online forecasting tech-nique provides a new tool for managing operation risks andsolving stochastic optimization problems. See, for example,[3]. For market participants, on the other hand, congestion andLMP forecasts by the operator provide actionable signal formanaging flexible resources and demand side management.Currently, there are very few probabilistic forecasting tech-niques for system operators beside direct Monte Carlo simu-lations [8]. The approach presented here represents a first steptoward online forecasting in large power systems with signif-icant stochastic components. For future work, there is a needto develop computationally tractable techniques, informativeperformance measure, and a set of practical benchmarks.A
PPENDIX P ROOF OF T HEOREM Proof.
For any parameter θ , denote the optimal primal anddual solutions to MPP (3) by x ∗ ( θ ) and y ∗ ( θ ) respectively.(1) For the MPLP case, if θ is in the same critical region as θ , they have the same optimal partitions, which means that A I x ∗ ( θ ) − b I − E I θ = , (21) A I c x ∗ ( θ ) − b I c − E I c θ < . (22)Because MPLP is neither primal nor dual degenerate, A I hasfull rank, and x ∗ ( θ ) = A − I ( b I + E I θ ) . (23)Substituting x ∗ ( θ ) into (22), we have θ ∈ Θ I .(2) For the MPQP case, the first order optimality conditionsare given by Hx ∗ ( θ ) + A ⊺ y ∗ ( θ ) = , (24) Ax ∗ ( θ ) − b − Eθ ≤ , (25) y ∗ i ( θ )( A i x ∗ ( θ ) − b i − E i θ ) = 0 , ∀ i ∈ J , (26) y ∗ ( θ ) ≥ . (27)From (24), x ∗ ( θ ) = − H − A ⊺ y ∗ ( θ ) . (28) Substituting the result into (26), we have A i H − A ⊺ i y ∗ i ( θ ) + b i + E i θ = 0 , ∀ i ∈ I , (29) y ∗ i ( θ ) = 0 , ∀ i ∈ I c . (30)By the non-degeneracy assumption, the rows of A I arelinearly independent. This implies that A I H − A ⊺ I is a squarefull rank matrix. Therefore, from (29), we solve y ∗ I ( θ ) = − ( A I H − A ⊺ I ) − ( b I + E I θ ) (31)and substitute y ∗ I ( θ ) and y I c ( θ ) into (28) to obtain x ∗ ( θ ) = H − A ⊺ I ( A I H − A ⊺ I ) − ( b I + E I θ ) . (32)Substituting x ∗ ( θ ) from (32) in the primal feasibility con-ditions (25) gives P p and substituting x ∗ ( θ ) from (32) in thedual feasibility condition (27) gives P d . We therefore have θ ∈ Θ I . A CKNOWLEDGEMENT
The authors thank Kory Hedman for pointing out the needfor incorporating time varying constraints and anonymousreviewers for their detailed comments.R
EFERENCES[1] Y. Ji, R. J. Thomas, and L. Tong, “Probabilistic forecast of real-timeLMP via multiparametric programming,” in
Proc. of the 48th HawaiiInternational Conference on System Sciences , 2015, pp. 2549 – 2556.[2] C. Auxenfans, R. Baumann, T. Lapetanovic, P. D. Leener, R. Paprocki,A. Wirth, and D. Klaar. (2014) European transmission system operators’cooperation for enhanced system security, integration of renewablesand market support. [Online]. Available: http://digilib.monenco.com/documents/10157/2529863/C2 113 2014.pdf[3] Y. Ji and L. Tong, “Stochastic coordinated transaction scheduling viaprobabilistic forecast,” in
Proc. of IEEE Power and Energy SocietyGeneral Meeting
International Journal of Forecasting
Proc. of the 3rd InternationalConference on Electric Utility Deregulation and Restructuring andPower Technologies , 2008, pp. 764–770.[9] Y. Ji, J. Kim, R. J. Thomas, and L. Tong, “Forecasting real-timelocational marginal price: A state space approach,” in
Proc. of the 47thAsilomar Conference on Signals, Systems, and Computers , 2013, pp.379–383.[10] Q. Zhou, L. Tesfatsion, and C. Liu, “Short-term congestion forecastingin wholesale power markets,”
IEEE Transactions on Power Systems ,vol. 26, no. 4, pp. 2185–2196, 2011.[11] R. Bo and F. Li, “Probabilistic LMP forecasting considering loaduncertainty,”
IEEE Transactions on Power Systems , vol. 24, no. 3, pp.1279–1289, 2009.[12] A. Rahimi Nezhad, G. Mokhtari, M. Davari, S. Hosseinian, andG. Gharehpetian, “A new high accuracy method for calculation of LMPas a random variable,” in
International Conference on Electric Powerand Energy Conversion Systems , 2009, pp. 1–5.[13] M. Davari, F. Toorani, H. Nafisi, M. Abedi, and G. Gharepatian, “De-termination of mean and variance of LMP using probabilistic DCOPFand T-PEM,”
PECon08 , pp. 1280–1283, 2008.[14] A. L. Ott, “Experience with PJM market operation, system design, andimplementation,”
IEEE Transactions on Power Systems , vol. 18, no. 2,pp. 528–534, 2003. [15] F. Borrelli, A. Bemporad, and M. Morari, “Geometric algorithm formultiparametric linear programming,”
Journal of Optimization Theoryand Applications , vol. 118, no. 3, pp. 515–540, 2003.[16] A. Bemporad, M. Morari, V. Dua, and E. N. Pistikopoulos, “Theexplicit linear quadratic regulator for constrained systems,”
Journal ofAutomatica , vol. 38, no. 1, pp. 3–20, 2002.[17] T. Gal and J. Nedoma, “Multiparametric linear programming,”
Manage-ment Science , vol. 18, pp. 406–442, 1972.[18] M. Herceg, M. Kvasnica, C. Jones, and M. Morari, “Multi-parametrictoolbox 3.0,” in
Proc. of the European Control Conference , Z¨urich,Switzerland, 2013, pp. 502–510, http://control.ee.ethz.ch/ ∼ mpt.[19] K. W. Hedman, S. S. Oren, and R. P. O’Neill, “A review of transmissionswitching and network topology optimization,” in Proc. of IEEE Powerand Energy Society General Meeting
Journal of the American Statistical Association , vol.102, no. 477, pp. 359–378, 2007.[22] G. W. Brier, “Verification of forecasts expressed in terms of probability,”
Monthly Weather Review , vol. 78, pp. 1–3, 1950.[23] R. D. Zimmerman, C. E. Murillo-S´anchez, and R. J. Thomas, “MAT-POWER: Steady-state operations, planning, and analysis tools for powersystems research and education,”