An Extended Integral Unit Commitment Formulation and an Iterative Algorithm for Convex Hull Pricing
aa r X i v : . [ ee ss . S Y ] J un i An Extended Integral Unit CommitmentFormulation and an Iterative Algorithm for ConvexHull Pricing
Yanan Yu, Yongpei Guan, Yonghong Chen
Abstract —To increase market transparency, independent sys-tem operators (ISOs) have been working on minimizing upliftpayments based on convex hull pricing theorems. However,the large-scale complex systems for ISOs bring computationalchallenges to the existing convex hull pricing algorithms. In thispaper, based on the analysis of specific generator features inthe Midcontinent ISO (MISO) system, besides reviewing integralformulations for several special cases, we develop two integralformulations of a single generator that can capture these features.We then build a compact convex hull pricing formulation basedon these integral formulations. Meanwhile, to improve the com-putational efficiency, we propose innovative iterative algorithmswith convergence properties, plus a complementary algorithm,to obtain a convex hull price. The computational results indicatethat our approach leads to an exact convex hull price on MISOinstances with and without transmission constraints and thesolutions can be obtained within minutes. Index Terms —Convex Hull Pricing, Iterative Algorithm, Inte-gral Formulation. N OMENCLATURE
A. Set and Dimension G Set of all generators. T Set of operation time span. S Set of start-up states, i.e., S = { h(hot), w(warm),c(cold) } . B Set of buses.
B. Parameters D t Load at time t (MW). F is Start-up cost of generator i in state s ∈ S ($). G it No load cost of generator i at time t ($). Q kit Piecewise linear production cost interception co-efficient of generator i in segment k at time t ($). H kit Piecewise linear production cost slope coefficientof generator i in segment k at time t ($/MWh). P it /P it Min/max generation amount of generator i at time t (MW). L i /ℓ i Minimum-up/-down time limit of generator i (h). L i Maximum-up time of generator i (h). V uit /V eit Ramp-up/-down rate of generator i at time t (MW/h). V uit /V eit Start-up/shut-down ramp rate of generator i attime t (MW/h). T wi /T ci Down-time limit for warm/cold start. The start-upcost is hot-start cost if the shut-down time is less
Y. Yu and Y. Guan are with the University of Florida, Gainesville, FL32611. Y. Chen (Consultant Advisor) is with MISO, Carmel, IN 46032. Apreliminary study of this work is shown in an earlier version of this paper [1]. than T wi , warm-start cost if it is longer than T wi and less than T ci , and cold-start cost if it is longerthan T ci (h). T si Minimum number of OFF periods for generator i in start-up state s ∈ S (h). T si Maximum number of OFF periods for generator i in start-up state s ∈ S (h). S ′ i ( t ) /S i ( t ) Time-dependent start-up/shut-down cost for gen-erator i ($). S ′ i /S i Constant start-up/shut-down cost for generator i ($). κ t /̟ t Indicator parameter showing whether min-up/-down time constraint is enforced for a generatorat time t (“ ” if it is yes and “ ” otherwise). a j /b j Slope/intercept of segment j in a piecewise linearfunction ($/MWh / $). N Number of pieces in a piecewise linear cost func-tion.
C. Decision Variables u it On/Off status of generator i at time t . v it Start-up status of generator i at time t . e it Shut-down status of generator i at time t . δ sit Indicator variable, which is if generator i startsat time t in state s ∈ S . x it Generation amount of generator i at time t (MW). f it Production cost of generator i at time t ($). φ s Production cost of a generator at time s ($). w t Indicator variable, which is if an initial-ongenerator shuts down for the first time at time t +1 . y tk Indicator variable, which is if a generator startsup at time t and shuts down at time k + 1 . z tk Indicator variable, which is if a generator shutsdown at time t + 1 and starts up at time k . θ t Indicator variable, which is if a generator shutsdown at time t + 1 and never starts up again. q stk Generation amount of a generator at time s whenit starts up at time t and shuts down at time k + 1 (MW). φ stk Production cost of a generator at time s when itstarts up at time t and shuts down at time k + 1 ($). ¯ γ t ( b ) Convex hull price at time t (in bus b ) ($/MWh). β t ( b ) LMP at time t (in bus b ) ($/MWh). i I. I
NTRODUCTION
In the current U.S. day-ahead electricity market operatedby ISOs, unit commitment and economic dispatch (UCED)problems are solved to determine the generation amount ofeach generator and the market-clearing price. The locationalmarginal prices (LMP) for market clearance is calculatedbased on the optimal dual vector corresponding to the system-wide constraints in the ED problem, in which commitmentdecision variables are fixed at their optimal values. Since thecommitment statuses are fixed, the LMPs cannot cover thestart-up and no-load costs of the generators [2]. Accordingly,uplift payments are introduced to make up the possible lossand motivate the generation side market participants to complywith the commitment schedule provided by ISOs. Since theuplift payment cost is not transparent, ISOs aim to minimizethe uplift payment cost for their daily operations [3], [4].To minimize the uplift payments, a convex hull pricingapproach has been recently introduced and received significantattention [5]–[8]. This pricing approach aims to minimize theuplift payments over all possible uniform prices. This approachhas the advantage to provide an optimal uniform price for thesystem. However, on the other hand, it requires to obtain anoptimal Lagrangian dual multiplier for a mixed-integer UCEDproblem and has been computationally intractable for marketimplementation [2].To overcome this challenge, researchers and practitionershave explored two different approaches to derive prices thatcould provide good approximations for the exact convex hullprice. One approach is to develop efficient algorithms tosolve the Lagrangian dual problem of the UCED problem.For instance, in [9], [10], the authors developed an extreme-point sub-differential method to strengthen the traditionalLagrangian relaxation approach to obtain the price. In [11],a subgradient simplex cutting plane method is developed tosolve the dual problem by iteratively cutting off non-optimalsolutions. Those algorithms show some benefits, while theycan not guarantee convergence in polynomial time. Mean-while, researchers turn to the other approach, called the primalapproach. This approach utilizes the property that, with theconvex hull description of each generator and the convexenvelope of the cost function, we can solve a relaxed UCEDproblem, which is a linear program, and derive the convexhull price based on the optimal dual vector associated withthe system-wide constraints [6], [7]. The work in [7] gets anapproximation of the exact convex hull price by deriving atight linear program relaxation of the convex hull descriptionfor each generator. Readers are also referred to [12]–[14] forother tight formulations.In practice, MISO has implemented an approximation ofconvex hull pricing based on Lagrangian relaxation, namedextended locational marginal prices (ELMPs) [2], and recentlythe primal formulation in [7] has been shown promising inreducing uplift payment [15]. However, a lot of practicalfeatures like ramping rate constraints [16], time-dependentstart-up costs, and so on, are not captured in [7], which makesthe formulation not tight for quite a few generators.In this paper, we follow the spirit of [6], [7] and study the generalization of convex hull pricing problem with thefocus on the real-world MISO system, which has more featuresthan those described in the literature, including time-variantcapacity, time-variant ramp-up/-down limit, flexible min-up/-down time limit, max-up time limit constraints, and time-dependent start-up costs. Those features are important toimprove the model accuracy, while they make it more difficultto describe the convex hull formulation. Therefore, it bringschallenges to the primal convex hull formulation framework.Besides, MISO operates one of the world’s largest energymarkets with more than $29 billion in annual gross marketenergy transactions [17]. It is required for the market-clearingprocess to solve large-scale problems with millions of decisionvariables [18]. Efficient algorithms are desired to achieve highcomputational efficiency. The main contributions of this paperare as follows:1) Based on the specific features of the MISO system,including four types of thermal generators, we reviewconvex hull formulations for two types of generators anddevelop two new convex hull formulations for the othertwo types of relative more complicated generators. Forthe class of generators with most complicated physicaland operational restrictions, we develop the correspond-ing integral formulation, which can capture time-variantcapacity, time-variant ramp-up/-down limit, flexible min-up/-down time limit, and max-up time limit constraints, aswell as time-dependent start-up costs altogether, besidestraditional physical restrictions. These features have notbeen captured altogether by any formulations in theliterature. Our integral formulation could lead to an exactconvex hull price for the MISO system by solving a linearprogram.2) To solve large-size MISO instances, we develop an inno-vative iterative algorithm, as well as its variant, to speedup the process to solve the problem. We prove that theiterative algorithm converges as the number of iterationsincreases. Furthermore, we develop a complementary al-gorithm to utilize several processors to solve the problemtogether, which can provide MISO the flexibility to get abetter solution within a given time limit.3) We test the algorithms on MISO instances. The numericalexperiment results show that our proposed algorithmslead to an exact convex hull price and a minimum upliftpayment for all of the test instances within an acceptabletime limit, which is one hour, the time limit given forexecuting pricing run for MISO’s day-ahead market, forthe cases with and without transmission constraints.The remainder of this paper is organized as follows. First, inSection II, we review the convex hull pricing framework andprimal formulation method. Next, in Section III, we adopt theintegral formulations in the literature that can be used in MISOconvex hull pricing problem. Then, in Section IV, we developtwo new integral formulations that can capture the specialcharacteristics of MISO instances. After that, in Section V, wedescribe our innovative iterative algorithm and its variants, aswell as the convergence proof of this algorithm. In Section VI,we report computational results on MISO instances. Finally, ii we conclude our research in Section VII.II. C ONVEX H ULL P RICING AND U PLIFT P AYMENT M INIMIZATION
In this section, we briefly review the convex hull pricingproblem as described in [5]–[7] and give the correspondingprimal formulation.The system optimization problem for a |T | -period UCEDproblem run by an ISO can be abstracted as follows: Z ∗ QIP = min X i ∈G C i ( x i , u i , v i , e i ) (1a)st. Ex ≥ p, (1b) ( x i , u i , v i , e i ) ∈ X i , ∀ i ∈ G , (1c)where C i ( x i , u i , v i , e i ) represents the total cost for generator i ,including start-up, shut-down, no-load, and generation costs,constraints (1b) represent system-wide constraints includingload balance and angle-eliminated transmission constraints interms of the shift factors [19], constraints (1c) represent thefeasible region of each generator i .After solving the UCED problem, an ISO obtains thegeneration amount ( ¯ x ) and commitment status ( ¯ u, ¯ v, ¯ e ) ofeach generator. The ISO also determines the price based onthe dual vector γ (referred to as price γ for brevity) for allparticipants. For a given price γ offered by the ISO, the self-scheduling profit maximization problem for each generator i can be described as follows: Z i ( γ ) = max γ ⊺ E i x i − C i ( x i , u i , v i , e i ) (2a)st. ( x i , u i , v i , e i ) ∈ X i . (2b)The uplift payment is defined as follows [6]: U = P i ∈G Z i ( γ ) − ( γ ⊺ E ¯ x − X i ∈G C i (¯ x i , ¯ u i , ¯ v i , ¯ e i ))+( γ ⊺ E ¯ x − γ ⊺ p ) , (3)where the first two items represent the difference between themaximum profit obtained through self-scheduling given price γ and the profit obtained following ISO’s schedule, and thelast item is defined as “Financial Transmission Rights (FTR)uplift” in [5].As stated in [7], the convex hull pricing problem corre-sponding to problem (1) is presented as the following formu-lation (4). The exact convex hull price can be derived as ¯ γ ⊺ E based on the optimal dual vector ¯ γ of (1b) in formulation (4)as introduced in [19] and [20, Chapter 8.11] (referred toas convex hull price ¯ γ for brevity). The minimum upliftpayment is guaranteed under the convex hull price ¯ γ . Notehere, conv ( X i ) represents the convex hull formulation for set X i , and ¯ C i ( x i , u i , v i , e i ) is the convex envelope for the generalcost function of each generator i . Z ∗ QLP = min X i ∈G ¯ C i ( x i , u i , v i , e i ) (4a) (¯ γ ) st. (1b) , ( x i , u i , v i , e i ) ∈ conv ( X i ) , ∀ i ∈ G . (4b) III. F ORMULATIONS FOR S EVERAL S PECIAL C ASES
There are different types of generators in practice withinMISO. We first present a traditional 3-bin UC formulation asa base for building our model as follows: Z ∗ QIP = min f,x,u,v,e,δ X i ∈G X t ∈T X s ∈ S F is δ sit + S i e it + G it u it + f it ! (5a)s.t. (1b) ( f i , x i , δ i , u i , v i , e i ) ∈ X i , ∀ i ∈ G , (5b)where the objective function is to minimize the total cost, inwhich the four items represent start-up, shut-down, no-load,and generation costs. The feasible region X i of generator i can be expressed as follows: X i = { f i , x i , δ i , u i , v i , e i ∈ R |T| × R |T| × R |T| × B |T| × B |T|− × B |T |− | f it ≥ H kit x it − Q kit u it , ∀ k ∈ [1 , N ] Z , ∀ t ∈ T , (6a) δ sit ≤ T is X j = T si e i ( t − j ) , ∀ s ∈ S / { c } , t ∈ T , (6b) X s ∈S δ sit = v it , t ∈ T , (6c) u it − u i ( t − = v it − e it , ∀ t ∈ [2 , |T | ] Z , (6d) t X j = t − L i +1 v ij ≤ u it , ∀ t ∈ [ L i + 1 , |T | ] Z , (6e) t X j = t − ℓ i +1 v ij ≤ − u i ( t − ℓ i ) , ∀ t ∈ [ ℓ i + 1 , |T | ] Z , (6f) x it ≥ P i u it , ∀ t ∈ T , (6g) x it ≤ P i u it , ∀ t ∈ T , (6h) x it − x i ( t − ≤ V ui u i ( t − + V ui v it , ∀ t ∈ T , (6i) x i ( t − − x it ≤ V ei u it + V ei e it , ∀ t ∈ T } , (6j)where P i , P i , V ui , V ui , V ei , V ei represent time-invariant parame-ters. Constraints (6a) represent the piecewise linear approxi-mation of the convex cost functions. Constraints (6b) and (6c)use indicator variables to represent each start-up type, whichdepend on the number of time periods the generator has beenoff before it was started-up. Constraints (6d) represent thelogic relationships. Constraints (6e) to (6j) represent the min-up/-down time restrictions (i.e., (6e)-(6f)), generation lowerand upper bounds (i.e., (6g)-(6h)), and ramp-up/-down rates(i.e., (6i)-(6j)).Now we introduce convex hull results for two types ofgenerators in MISO below.
1) The set of generators with constant start-up costs andwithout ramping constraints (labeled as set G ): Amongvarious types of generators, low capacity fast-start gas gen-erators are relatively easy to model in terms of physicalconstraints. They can start-up quickly, and the ramping ca-pability is larger than the gap between upper and lowergeneration limits. Also, the start-up costs are constant and notdependent on the down-time periods before the start-up. Thus, v constraints (6b), (6c), (6i), and (6j) are redundant and (6) canbe simplified as follows: Z ∗ QIP = min x,u,v,e X i ∈G X t ∈T ( S ′ i v it + S i e it + G it u it + f it ) (7)s.t. (1b) , ( f i , x i , u i , v i , e i ) ∈ X i , ∀ i ∈ G . X i = { f i , x i , u i , v i , e i ∈ R |T | × R |T | × B |T | × B |T |− × B |T |− | (6a) , (6d) − (6h) } . (8)The convex hull (as described in [7] and [21]) for X i is D i = { f i , x i , u i , v i , e i ∈ R |T | × R |T | × R |T | × R |T |− × R |T |− | (6a) , (6d) − (6h) ,v i ≥ , e i ≥ } . (9)
2) The set of generators with constant start-up costs andstart-up ramping constraints (labeled as set G ): For anothergroup of gas-fired generators with less restrictive physicalconstraints, the ramping rates during the stable region arelarger than the gap between the upper and lower generationbounds. But the start-up ramping rates ( V ui ) are binding. Forthis case, the convex hull description is (as described in [21]) D i = { f i , x i , u i , v i , e i ∈ R |T | × R |T | × R |T | × R |T |− × R |T |− | (6a) , (6d) − (6h) , (9) ,x it ≤ P i u it + ( V ui − P i ) v it , ∀ t ∈ [2 , |T | ] Z } . (10)IV. F ORMULATIONS FOR THE G ENERAL
MISO I
NSTANCES
For the generators in MISO, besides special generators G and G as described in III, there are several generator featureswhich reflect the market needs and capture more details inpractice. In this section, we first introduce three new featuresin general MISO instances. Then based on generator types, wederive two convex hull descriptions capturing these features. A. Features in the general MISO instances1) Max-up time limit:
For some generators in the MISOmarket, there are restrictions on maximum time periods thatthe generator can stay online because of machine deterioration.To accommodate this, the max-up time constraints can bedescribed as follows: t + L i X j = t +1 v ij ≥ u i ( t + L i ) , ∀ t ∈ T (11)
2) Flexible min-up/-down time limit:
In MISO, the min-up/-down time limit is set to be time-variant to resolve offerdata conflicts. For example, participants may submit a must-run offer for a generator for hours − and − withthe min-down time limit as hours. This will force UCEDto commit this generator between and even if it is costly.So MISO developed a set of rules to ignore min-up/-downtime limit if there are such conflicts. In this example, the min-down time limit is relaxed to be between hours and so that the generator will not be committed if it is costly.It will force market participants to submit proper offers andprevent market manipulation. If they do want to run through,they should submit a must-run offer for all hours. If they do want MISO to determine on/off in between, they should haveat least hours in between. To accommodate this, the refinedmin-up/-down time constraints in the -bin formulation (suchas G and G ) can be described as follows (generator index i is omitted for brevity): t X j = t − L +1 κ t v j ≤ u t , ∀ t ∈ [ L + 1 , |T | ] Z , (12) |T | X j = t − ℓ +1 ̟ t v j ≤ − u t − ℓ , ∀ t ∈ [ ℓ + 1 , |T | ] Z . (13)
3) Time-variant parameters:
In MISO, market participantsare allowed to offer capacity and ramp rates varying by thehour. The time-variant parameters make the convex hull morecomplicated and have rarely been studied in the literature.
B. Convex hull results for general MISO instances1) The set of generators in G with max-up time limitconstraints described in IV-A1 (labeled as set G ): For thistype of generators, we derive and prove the convex hulldescription in Theorem 1.
Theorem 1.
The convex hull description for the 3-bin modelwith max-up time restriction can be described as follows: D i = { f i , x i , u i , v i , e i ∈ R |T | × R |T | × R |T | × R |T |− × R |T |− | (6a) , (6d) − (6h) , (9) − (11) } . (14) Proof.
Based on Proposition in [22], for the min-up/-down time only polytope (without ramping constraints) withmax-up time restrictions, the convex hull description of thefeasible binary variables ( u, v, e ) is { u, v, e ∈ R |T | × R |T | × R |T | | (6d) − (6f) , (9) , (11) } . Then, based on Lemma in [21],the addition of capacity constraints (6g)-(6h), start-up rampingconstraints (10), and linear function (6a) does not affectintegrality. Thus, the conclusion holds.
2) The set of generators with all the additional restrictionsdescribed in IV-A1-IV-A3(labeled as set G ): For this type ofgeneral generators, the model is built in a high-dimensionalspace by introducing binary variables y tk ( z tk ) to indicatea generator is “ON” (“OFF”) from time periods t to k ,in order to keep track of its “ON” and “OFF” intervals.Accordingly, a new decision variable q stk , corresponding toeach “ON” interval [ t, k ] Z , is introduced to represent thegeneration amount at time s, t ≤ s ≤ k . Now, we describethe integral formulation (referred to as EUC formulation) andgive the detailed illustrations in the proof of Theorem 2. Weassume the generator has been initially on for s time periodsbefore time . Thus, the generator cannot shut down until time t + 1 , with t = [ L − s ] + , due to min-up constraints. Thedetailed formulation is shown as follows: min X k ∈K ,k< |T | S ( k + s ) w k + X tk ∈T K ,k< |T | S ( k − t + 1) y tk + X kt ∈KT S ′ ( t − k − z kt + X tk ∈T K k X s = t φ stk (15a) s.t. X k ∈K w k = 1 , (15b) − w t ′ { t ′ ∈K} + X tk ∈KT ,t = t ′ z t ′ k − X kt ∈T K ,t = t ′ y kt ′ + θ t ′ =0 , ∀ t ′ ∈ [ t , |T |− Z , (15c) X tk ∈T K ,t = t ′ y t ′ k − X kt ∈KT ,t = t ′ z kt ′ = 0 , ∀ t ′ ∈ [ t + ℓ t + 1 , |T | ] Z , (15d) P s w k ≤ q stk ≤ P s w k , ∀ s ∈ [ t, k ] Z , ∀ tk ∈ T K , (15e) P s y tk ≤ q stk ≤ P s y tk , ∀ s ∈ [ t, k ] Z , ∀ tk ∈ T K , (15f) q ttk ≤ V ut y tk , ∀ tk ∈ T K , (15g) q ktk ≤ V ek w k , ∀ tk ∈ T K , k ≤ |T | − (15h) q ktk ≤ V ek y tk , ∀ tk ∈ T K , k ≤ |T | − (15i) q s − tk − q stk ≤ V es w k , q stk − q s − tk ≤ V us w k , ∀ s ∈ [ t + 1 , k ] Z , ∀ tk ∈ T K , (15j) q s − tk − q stk ≤ V es y tk , q stk − q s − tk ≤ V us y tk , ∀ s ∈ [ t + 1 , k ] Z , ∀ tk ∈ T K , (15k) φ stk − a j q stk ≥ b j w k , ∀ s ∈ [ t, k ] Z , j ∈ [1 , N ] Z , ∀ tk ∈ T K , (15l) φ stk − a j q stk ≥ b j y tk , ∀ s ∈ [ t, k ] Z , j ∈ [1 , N ] Z , ∀ tk ∈ T K , (15m) w, z, y ≥ , θ t ≥ , ∀ t ∈ [ t , |T | − ℓ t − Z , (15n)where the objective function is to minimize the total cost,including shut-down, start-up and generation costs, con-straints (15b)-(15d) keep track of the generator’s “ON” and“OFF” statuses ( { t ′ ∈K} indicates that w t ′ is included in (15c)only when t ′ is in set K ), and constraints (15e)-(15m) describethe economic dispatch restrictions in the “ON” interval, includ-ing the generation upper and lower bounds (15e)-(15f), start-upand shut-down ramping restrictions (15g)-(15i), general ramp-up and ramp-down restrictions (15j)-(15k), and finally thepiecewise linear convex function f stk ( q stk ) representation (15l)-(15m). Meanwhile, considering the new features in termsof max-up and flexible min-up/-down time limit, we define K = [ t , min { L − s , |T |} ] Z as the subscript set of vari-able w , T K ( KT ) as all possible on-intervals (off-intervals)satisfying min-up and max-up (min-down) time restrictions.More specifically, we let T K = T K ∪ T K , where T K collects all possible first “ON” intervals [ t, k ] Z with t = 1 and k ∈ [ t +1 , min { L − s , |T |} ] Z and T K collects all following“ON” intervals [ t, k ] Z , after a shut-down has happened, with t ∈ [ t + ℓ t + 1 , |T | ] Z and k ∈ [min { t + L t − , |T |} , min { t + L − , |T |} ] Z . In this expression, we have i) L t = L if κ t = 1 and L t = 1 otherwise, and ii) ℓ t = ℓ if ̟ k = 1 and ℓ t = 1 otherwise to reflect the flexible min-up/-down time limit asdescribed earlier in Section IV-A2. Similarly, we have KT collecting all “OFF” intervals [ k, t ] Z with k ∈ [ t , |T |− ℓ t − Z and t ∈ [ k + ℓ t + 1 , |T | ] Z . The new feature in terms oftime-variant parameters is considered by allowing parameters P s , P s , V ut , V ek , V es to be dynamic in time. Theorem 2.
The convex hull description for the general singlegenerator MISO UC is as follows: D i = { w i , y i , z i , θ i , q i , φ i | (15b) − (15n) } , (16) and there exists an optimal solution to (15) binary with respectto decision variables w , y , and z . Proof. The detailed proof is provided in Appendix A.Finally, it can be observed that the convex hull descriptionfor the initial-off generators can be obtained similarly.V. I
TERATIVE A LGORITHMS FOR C ONVEX - HULL P RICINGFOR
MISO I
NSTANCES
Based on the above convex hull descriptions (9), (10), (14),and (16), we derive the general convex hull pricing formu-lation ( P ) as follows. For notation brevity, we let G = G / ( G ∪ G ∪ G ) represent all unclassified generators usingformulation (15). Since we have the convex hull descriptionsfor each type of generators and convex envelope for all costfunctions, we can provide the exact convex hull price andminimize the uplift payment by solving the following linearprogram: ( P ) : Z ∗ LPO = min f,x,u,v,e,w,y,z,θ,q,φ X i ∈G ∪G ∪G X t ∈T g it + X i ∈G g ′ it (17a)s.t. (1b) , ( f i , x i , u i , v i , e i ) ∈ D i , ∀ i ∈ G , ( f i , x i , u i , v i , e i ) ∈ D i , ∀ i ∈ G , ( f i , x i , u i , v i , e i ) ∈ D i , ∀ i ∈ G , ( w i , y i , z i , θ i , q i , φ i ) ∈ D i , ∀ i ∈ G , where the cost function g it of each generator i in G ∪ G ∪ G can be expressed as g it = S ′ i v it + S i e it + G it u it + f it (18)and the cost function g ′ it of each generator i in G can beexpressed as g ′ it = X k ∈K ,k< |T | S i ( k + s ) w ik + X tk ∈T K i ,k< |T | S i ( k − t + 1) y itk + X kt ∈KT S ′ i ( t − k − z ikt + X tk ∈T K k X s = t φ sitk . (19)The formulation ( P ) is sufficient to solve the convex hullpricing problem. However, a large number of variables andconstraints in the EUC formulation (i.e., (15)) increase thecomputational complexity and lead to a long solving timewhen the number of generators in G is large. To solve thelarge-scale problem ( P ), we develop an iterative algorithm,in which a relaxation problem where D i is replaced by D i and g ′ it is replaced by g it for each generator i ∈ G in ( P )is solved in the first step. Then the constraints D i in EUCformulation are added gradually when needed to tighten therelaxation. Meanwhile, the objective function g it is replacedby g ′ it when D i is added back. Our approach can provide avery tight approximation of the convex hull price, with mostof the cases converging at the optimal solution. Meanwhile thealgorithm terminates in a short time. A. Methodology background
Before describing the detailed algorithms, we first introduceLemma 1 and Theorem 3 to provide a theoretical foundationshowing that the uplift payment amount will decrease and i converge under our algorithm. For notation brevity, we use x ∈ R n to represent all variables in (5), which contains binaryvariable set x and continuous variable set x . The traditional3-bin MILP UC formulation (5) can be abstracted as follows: Z ∗ QIP = min { c ⊺ x | Ex ≥ p, x ∈ X } , (20) X = { x ∈ B n , x ∈ R n | Ax ≤ b } . Lemma 1.
For A ∈ R m × n , A ′ ∈ R m × n , E ∈ R p × n , b ∈ R m , b ′ ∈ R m , p ∈ R p and c ∈ R n , we consider the integeroptimization problem (20) , its tightened linear programmingrelaxation problem ( P C ) in which A ′ x ≤ b ′ dominates Ax ≤ b ,and the Lagrangian relaxation D C corresponding to a dualvector γ as follows: ( P C ) : Z C = min { c ⊺ x | A ′ x ≤ b ′ , Ex ≥ p, x ∈ R n } , (21) ( D C ) : Z C ( γ )= min { c ⊺ x + γ ⊺ ( p − Ex ) | A ′ x ≤ b ′ , x ∈ R n } . (22) Given an optimal dual vector ¯ γ for constraints Ex ≥ p in (21) ,if (i) x ∗ = { x ∗ , x ∗ } is an optimal solution to problem Z C (¯ γ ) ,(ii) x ∗ are all binaries,then the uplift payment given price ¯ γ can be calculated as U = Z ∗ QIP − Z C (¯ γ ) ( ⊠ ).Proof. Based on the definition of uplift payment in (3), theprofit following the ISO’s schedule can be calculated as P ISO = ¯ γ ⊺ E ¯ x − Z ∗ QIP ( ♦ ), which is the second item in (3),and the maximum profit P i ∈G Z i (¯ γ ) obtained through self-scheduling problem (2) given price ¯ γ can be derived by solvingthe following problem, where c ⊺ x is an abstract form of C i ( x i , u i , v i , e i ) in (2): X i ∈G Z i (¯ γ ) = − min { c ⊺ x − ¯ γ ⊺ Ex | A ′ x ≤ b ′ , x ∈ B n } . (23)Based on formulations (22) and (23), it is clear that − P i ∈G Z i (¯ γ ) ≥ Z C (¯ γ ) − ¯ γ ⊺ p ( ⋆ ), since x is relaxed tobe continuous between and in (22). Meanwhile, whenconditions (i) and (ii) hold, x ∗ is a feasible solution to (23).Thus, − P i ∈G Z i (¯ γ ) ≤ Z C (¯ γ ) − ¯ γ ⊺ p ( ♣ ). Combining ( ⋆ ) and( ♣ ), we have P i ∈G Z i (¯ γ ) = − ( Z C (¯ γ ) − ¯ γ ⊺ p ) ( > ).Therefore, U = P i ∈G Z i (¯ γ ) − P ISO +¯ γ ⊺ E ¯ x − ¯ γ ⊺ p = − ( Z C (¯ γ ) − ¯ γ ⊺ p ) − (¯ γ ⊺ E ¯ x − Z ∗ QIP )+¯ γ ⊺ E ¯ x − ¯ γ ⊺ p = Z ∗ QIP − Z C (¯ γ ) ,where the first equation holds due to (3), the second equationholds based on ( > ) and ( ♦ ), and the last one is directly derivedby eliminating the same items. Remark 1.
Note here that the prices ¯ γ are calculated bysolving a linear program, and thus, they are not linked to theoptimal integer solutions in the original problem. Remark 2.
When Z ∗ QIP is not obtained in practice, the upliftpayment is calculated using equation U = Z ∗ QIP − Z C (¯ γ ) , i.e.,( ⊠ ), by replacing Z ∗ QIP with its best upper bound Z ′ QIP obtainedby the MIP solver. The derived uplift payment using Z ′ QIP provides an upper bound of the optimal uplift payment, since Z C (¯ γ ) does not change for a given ¯ γ and Z ′ QIP ≥ Z ∗ QIP . Theextra amount of uplift payment is equal to Z ′ QIP − Z ∗ QIP . Theorem 3.
For A ∈ R m × n , A ∈ R m × n , E ∈ R p × n , b ∈ R m , b ∈ R m , d ∈ R p , and c ∈ R n , considering the integeroptimization problem (20) , its linear programming relaxation z = min { c ⊺ x | Ax ≤ b, Ex ≥ p, x ∈ R n } (for notation brevity,we define X L = { x ∈ R n | Ax ≤ b } ), and two tightened linearprogramming relaxation problems ( P C ) , ( P C ) shown below, ( P C ) : Z C = min { c ⊺ x | Ex ≥ p, x ∈ X } , (24) ( P C ) : Z C = min { c ⊺ x | Ex ≥ p, x ∈ X } , (25) in which conv ( X ) ⊆X ⊆ X ⊆ X L with X = { x ∈ R n | A x ≤ b } and X = { x ∈ R n | A x ≤ b } . Underthese two formulations, we use the dual vector for constraints Ex ≥ p to derive the market clearing price, written as ¯ γ and ¯ γ . Then the uplift payment under ¯ γ will be no larger thanthat under ¯ γ if conditions (i) and (ii) in Lemma 1 hold.Proof. According to Lemma 1, the uplift payments underformulations (24) and (25) can be calculated as U = Z ∗ QIP − Z C (¯ γ ) and U = Z ∗ QIP − Z C (¯ γ ) when conditions (i) and (ii)hold. Based on strong duality in Lagrangian relaxation, wehave Z C = Z C (¯ γ ) and Z C = Z C (¯ γ ) . Since X ⊆ X , wehave Z C ≤ Z C . Thus, we can conclude U ≥ U . Remark 3. If X = conv ( X ) , then ¯ γ is a convex hull priceand Z ∗ QIP − Z C (¯ γ ) is the minimum uplift payment.B. The detailed algorithm We consider the following ( P
1) as the starting point: ( P
1) : Z C = min f,x,u,v,e X i ∈G X t ∈T g it s.t. (1b) , ( f i , x i , u i , v i , e i ) ∈ D i , ∀ i ∈ G , ( f i , x i , u i , v i , e i ) ∈ D i , ∀ i ∈ G , ( f i , x i , u i , v i , e i ) ∈ D i , ∀ i ∈ G ∪ G . In this formulation, we keep the convex hull descriptions forgenerators in G ∪ G ∪ G in problem ( P ), and relax theconvex hull set D i for each generator in G to be D i . Asstated above, adding EUC formulations (15) to all of thegenerators in G will lead to the perfect formulation, but itwill increase the computational complexity significantly. Theiterative algorithm is designed to select “necessary” generatorsin G and only add (15) of “necessary” generators into ( P ),which can effectively tighten the formulation and improve Z C with a less computational burden. The detailed algorithm isshown in Algorithm 1. Here we provide the explanation.Following the conditions in Lemma 1, we first solve ( P )and get the dual vector ¯ γ of equation (1b). Then we solve thefollowing ( P ), which is a Lagrangian relaxation problem of( P ) based on the given price ¯ γ : ( P
2) : Z C (¯ γ )= ¯ γ ⊺ p + min f,x,u,v,e P i ∈G P t ∈ T ( g it − ¯ γ ⊺ E i x it ) s.t. ( f i , x i , u i , v i , e i ) ∈ D i , ∀ i ∈ G , ( f i , x i , u i , v i , e i ) ∈ D i , ∀ i ∈ G , ( f i , x i , u i , v i , e i ) ∈ D i , ∀ i ∈ G ∪ G . To satisfy conditions (i) and (ii) in Lemma 1, we need toget an integral solution to ( P ). In ( P ), since there are no ii coupling constraints among generators, this problem is equiva-lent to maximizing profit for each generator independently andsumming them together. Assume the optimal solution to ( P )is ( f ∗ , x ∗ , u ∗ , v ∗ , e ∗ ) . For each generator j ∈ G ∪ G ∪ G ,the values of u ∗ j , v ∗ j , e ∗ j will be binaries for any given price ¯ γ , since D j , D j , and D j provide the convex hull descriptions.For each generator j ∈ G , we may get fractional values in u ∗ j , v ∗ j , e ∗ j since D j is a relaxation of D j . In this algorithm,we consider the constraints in D j as a cutting plane group.When a fractional solution to a particular generator j ∈ G is obtained in ( P ), we add this cutting plane group D j to( P ) and ( P ) and replace the objective function g jt with g ′ jt .In this way, we can ensure that this generator has integralsolutions when we solve ( P ) with any given price ¯ γ . The set Γ (representing the set of generators in G using the constraintsin D j ) is updated accordingly by adding this generator. Sincethe dual vector ¯ γ could change after updating the formulationfor certain generators in G , we need to solve ( P ) again andget an updated ¯ γ . We thus create an inner loop to repeat thisprocess until we could not find an updated price ¯ γ .Next, we further decrease the uplift payment by tightening( P ) utilizing Theorem 3. The optimal solution to the updated( P ) could be infeasible for constraints in (15) for somegenerator i in G / Γ since D i is relaxed to D i . Given theoptimal solution ( ˆ f i , ˆ x i , ˆ u i , ˆ v i , ˆ e i ) of each generator i in ( P ),we can identify infeasible generators by solving the following( P ) for each generator i ∈ G / Γ with fractional ˆ u i , ˆ v i , or ˆ e i : ( P
3) : min w i ,z i ,y i ,θ i ,q i ,φ i (28a)s.t. ( w i , z i , y i , θ i , q i , φ i ) ∈ D i , ˆ x is = X tk ∈T K ,t ≤ s ≤ k q stk , ˆ f is = X tk ∈T K ,t ≤ s ≤ k φ stk , ∀ s ∈ T , (28b) ˆ v is = X kt ∈KT ,t = s z kt , ˆ e is = X kt ∈KT ,k = s z kt + θ s , ∀ s ∈ T , (28c) ˆ u is = X tk ∈T K ,k ≥ s w k + X tk ∈T K ,t ≤ s ≤ k y tk , ∀ s ∈ T , (28d)where constraints (28b)-(28d) build the linear mapping be-tween the optimal solution to EUC formulation (15) and thatof the traditional 3-bin MILP UC formulation (6). If the above( P ) is infeasible for generator i , then this generator i willbe included in set Γ . We update problems ( P ) and ( P )by adding the constraints in D i and replacing g it with g ′ it ,which tightens ( P ) by cutting off the infeasible solutionand meanwhile ensure integral solutions in ( P ) for thosegenerators. Note here that the generators with an integralsolution to ˆ u i , ˆ v i , and ˆ e i are always feasible in ( P ) since theconvex hull D contains all of the feasible integral solutions.Similarly, a new dual vector ¯ γ could be obtained aftersolving the updated ( P ), and it may lead to extra generatorswith fractional solutions in updated ( P ). Thus we create anouter loop to repeat this process iteratively to ensure that wehave an optimal solution feasible to the perfect formulation( P ). The iterative algorithm terminates when we obtain ansolution feasible to ( P ), and an integral solution under thecorresponding ¯ γ is obtained in ( P ). The resulting upliftpayment under price ¯ γ is U = Z ∗ QIP − Z C (¯ γ ) . We can observe that this algorithm satisfies the conditionsin Theorem 3 and thus, the uplift payment is non-increasing ineach iteration. Meanwhile, the uplift payment is bounded be-low by the minimum uplift payment. Therefore, the algorithmconverges based on the monotone convergence theory. Remark 4.
Note here that although the outer loop could guar-antee the constraints of our integral formulation (15) satisfied,the objective of the final ( P ) in the outer loop could be stilldifferent from that of ( P ), because our integral formulation hasa different objective function (19) , as compared to that of theoriginal formulation (18) . Thus, the uplift payment calculationin the algorithm may not reflect the actual uplift paymentamount, i.e., U = Z ∗ QIP − Z C (¯ γ ) in general under this case.Accordingly, this does not guarantee the monotone decrementof the uplift payment. Therefore, solving ( P ) is necessary tohave Theorem 3 held to ensure convergence.C. Complementary algorithm It can be first observed that the outer loop, which tightens( P ) and the inner loop, which ensures conditions in Lemma 1held in ( P ), can be exchanged. The exchange may resultin a different converging path and lead to a different result.We denote these two as IA1 and IA2. To further tighten theresults from IA1/IA2, by taking advantage of the independenceof ( P ) to the other parts of the algorithm, we developa complementary algorithm, which can be implemented inparallel to check and select the potential candidates from G / Γ in the resulting ( P ) after IA1/IA2 to add the corresponding D i , which could further tighten the formulation.Assuming there are M + 1 parallel computing nodes withone master and M slave nodes, the steps of the complementaryalgorithm, denoted as IAC1 and IAC2, are as follows:Step 1: Run Algorithm IA1/IA2 on the master node. After itfinishes, record the resulting ( P ), price ¯ γ , set Γ , and Z C (¯ γ ) ;Step 2: Divide the generators in G / Γ into M groups: G ′ , . . . , G ′ M . The resulting ( P ) and G ′ m are distributed toeach slave node m , ∀ m ∈ { , . . . , M } ;Step 3: For each slave node m ∈ { , . . . , M } , eachgenerator j ∈ G ′ m is updated in the resulting ( P ) in sequencefollowing Steps 12-13 in IA1/IA2. If the optimal objectivevalue is improved by adding a generator j , the resulting ( P )in the master node will be updated accordingly and solved toget improved uplift payment U and price ¯ γ ;Step 4: The stopping rules are flexible. The algorithm willstop if one of the following rules is satisfied: i) the masternode receives n updated generators from the slave nodes, ii)the time limit T limit is reached, or iii) all of the slave nodesterminate.The final updated convex hull price ¯ γ ′ and Z C (¯ γ ′ ) can beused to calculate the uplift payment as U ′ = Z ∗ QIP − Z C (¯ γ ′ ) .VI. C OMPUTATIONAL E XPERIMENTS ON
MISOI
NSTANCES
In this section, we report the performance of our proposedformulations and algorithms for the MISO system. More iii
Algorithm 1:
An iterative algorithm (IA1)
Data:
The parameters of the system
Result:
Convex hull price ¯ γ and uplift payments U Initialization: initialize ( P ) and construct set Γ = ∅ ; Solve MILP (5) with the optimal objective value Z ∗ QIP ; Relax variables u, v, e to be continues; do Solve ( P ) and get the optimal solution ( ˆ f , ˆ x, ˆ u, ˆ v, ˆ e )and dual vector ¯ γ ; do Set count number m = 0 ; for j ∈ G / Γ do Given dual vector ¯ γ , solve ( P ) and get theoptimal solution ( f ∗ j , x ∗ j , u ∗ j , v ∗ j , e ∗ j ) and theoptimal objective value Z C (¯ γ ) ; if there exists fractional values in terms of u ∗ j , v ∗ j , e ∗ j in the optimal solution to ( P ) then Γ ← Γ ∪ { j } , m ← m + 1 ; Replace the objective function g jt in ( P )and ( P ) with g ′ jt ; Replace the corresponding constraint setin ( P ) and ( P ) with D j ; Solve the updated ( P ) and get the optimalsolution ( ˆ f , ˆ x, ˆ u, ˆ v, ˆ e ) and dual vector ¯ γ ; while m != 0 ; Set count number n = 0 ; for i ∈ G / Γ do if there exists fractional values in terms of ˆ u i , ˆ v i , ˆ e i in the optimal solution to ( P ) then Solve ( P ) for generator i and check thefeasibility; if Infeasible then Γ ← Γ ∪ { i } ; n ← n + 1 ; Replace the objective function g it in ( P )and ( P ) with g ′ it ; Replace the corresponding constraint setin ( P ) and ( P ) with D i ; while n != 0 ; The estimated convex hull price is ¯ γ and the upliftpayment is U = Z ∗ QIP − Z C (¯ γ ) .specifically, we test the performance of the proposed tightlinear programming formulation, the iterative algorithms, andtheir complementary parts. Each MISO instance includes over , generators with their corresponding bidding informationprovided. Among all of the generators, the percentages of G , G , G , and G are around , , , and . Werandomly select instances each with a -hour operationplanning horizon and compute the corresponding convex hullprice and uplift payment. All test instances were run on a dual -processor Intel(R) Xeon(R) CPU E - v @ . GHz
GB with Gurobi . . as the optimization solver. Thedefault MIP optimality gap is set to be 1e-3. We report the results for the following models and algo-rithms: • MIP: the traditional -bin MILP UC model (5); • LMP: energy price is calculated based on the dual vec-tor corresponding to the system-wide constraints in theUCED problem when commitment statuses are fixed attheir optimal values; • TLP: use the approximated convex hull pricing formula-tion ( P ) to obtain the price; • IA1, IA2: approximated convex hull pricing formulationusing the algorithms described in Section V to obtain theprice; • IAC1, IAC2: approximated convex hull pricing formula-tion using the iterative algorithm plus its complementaryalgorithm, implemented in a single processor. The set-ting for the complementary algorithm is as follows: 1) M = |G / Γ | , 2) n = 2 , and 3) T limit = 200 s ; • OPT: use the exact convex hull pricing formulation ( P )to obtain the convex hull price.We tested the cases with and without transmission con-straints, respectively. The results for the cases without trans-mission constraints are reported in Table I. In the table, column“Solution ($)” represents the optimal cost for the systemoptimization model, with “MIP” representing the optimal ob-jective for the MILP UC model (5) and others representing theoptimal objective for the corresponding linear programmingrelaxation, column “Uplift payment ($)” represents the totaluplift payment generated from each approach, column “Time(s)” represents the computational time in terms of seconds foreach approach. Since IAC1 and IAC2 are complementary toIA1 and IA2, we only report additional time required to reachthe optimal solution. That is, we use “ ⋄ ” to indicate that IAC1or IAC2 stops after the stopping rule is satisfied, and “(+ δ )”after “ ⋄ ” to indicate that the complementary algorithm takes δ extra seconds after the iterative algorithm terminates to reachthe optimal solution. “(+0)” shows that the result of IA1 or IA2is already optimal. Column “Save ($)” represents the upliftpayment savings as compared to the “LMP” approach. Thesavings of the TLP approach are calculated as the differenceof uplift payments between the TLP approach and the LMPapproach. Since IA1, IA2, IAC1, and IAC2 tighten the TLPmodel, the savings of these four methods are represented as theextra savings beyond the TLP approach. For example, in case“C1”, the TLP approach can save $3 , as compared to theLMP approach. The IA1 approach can save extra $299 , whichmeans the IA1 approach can save $3 ,
093 + $299 = $3 , intotal as compared to the LMP approach. The “OPT” providesa benchmark to indicate the efficiency of each approach. Forinstance, our proposed algorithms achieve an exact convex hullprice and a minimum uplift payment if the uplift paymentderived from our approach is the same as that from the “OPT”algorithm. Column “Optimal” indicates if an exact convexhull price and a minimum uplift payment are achieved byeach approach. Finally, column “Diff ($/MWh)” represents theaverage absolute deviation between the exact convex hull andLMP prices. For the cases without transmission constraints,the “Diff ($/MWh)” is calculated as P t ∈T | ¯ γ t − β t | / |T | with x the convex hull prices ¯ γ t and LMP prices β t , ∀ t ∈ T .From Table I, we can first observe that the TLP approachcan effectively reduce the uplift payments with a shorter com-putational time as compared to the LMP approach. The pricecan be obtained without solving the UC problem. Second, wecan observe that our proposed alternative algorithms IA1 andIA2 are very effective. They provide extra savings based onthe TLP approach. Both IA1 and IA2 can solve nine out of theeleven total cases into optimality. For the remaining two cases,at least one of IA1 and IA2 can achieve the minimum upliftpayment and the exact convex hull price. In practice, ISOs canrun both IA1 and IA2 to improve the chance to obtain the exactconvex hull price. More importantly, the algorithm ran veryfast and can terminate within minutes for all cases. Third,for the cases whose optimal values are not reached by IA1 (orIA2), the complementary algorithm can effectively find out the“key” generators and add them into IA1 (or IA2) to ensure theoptimality. Among the eleven cases at most two generators areadded after terminating IA1 (or IA2), which also verified thecompactness of IA1 and IA2. What is more, the stopping rulesfor IAC1 and IAC2 can be flexible in terms of adding extragenerators into the formulation, which can provide ISOs theflexibility to get a proper price within a time limit. Finally, wecan observe that the average absolute deviations between theconvex hull prices and LMP prices range from . to . .The results for the cases with transmission constraints (C10and C11) are reported in Table II. For the cases with transmis-sion constraints, the price is a vector with dimension |T |×|B| ,where |B| is the number of buses in the MISO system, whichis over . In column “Diff ($/MWh),” we report the av-erage absolute deviation P t ∈T P b ∈B | ¯ γ tb − β tb | / ( |T | × |B| ) ,where ¯ γ tb and β tb are convex hull prices and LMP prices attime t ∈ T in bus b ∈ B . The results are also very promising,and similar observations are obtained. From Table II, we canobserve that both IA1 and IA2 achieve the minimum upliftpayment and exact convex hull price for the first case. Forthe second case, IA2 achieves the minimum uplift payment,and IA1 achieves it with the help of the complementaryalgorithm. Compared with the results for Cases C10 and C11without transmission constraints, the savings for the caseswith transmission constraints are relatively larger, and thedifferences between the exact convex hull prices and LMPprices are also relatively larger. Meanwhile, it does not takemuch longer to solve the cases with transmission constraintswith all finished in minutes.Finally, from Tables I and II, by comparing the solvingtime of our proposed algorithms with the OPT approach,we can observe that our proposed algorithms are much morecomputationally efficient.VII. C ONCLUSION
This paper has provided a customized compact formulation,efficient iterative algorithms, and implementation techniques tosolve the convex hull pricing problem for MISO instances. Tocapture the real-world generator characteristics, we categorizefour groups of generators and develop two new convex hulldescriptions, in which the EUC formulation provides the con-vex hull description for the generators with most complicated
TABLE IT
EST RESULTS FOR
MISO
WITHOUT TRANSMISSION CONSTRAINTS
Case Model Solution($) Uplift Payment($) Time(s) Save($) Optimal Diff($/MWh)C1 MIP 39,986,855 - - - - 0.25LMP - 3,521 36 - NTLP 39,980,331 428 15 3,093 NIA1 39,986,726 129 241 +299 YIA2 39,986,726 129 247 +299 YIAC1 39,986,726 129 ⋄ (+0) +299 YIAC2 39,986,726 129 ⋄ (+0) +299 YOPT 39,986,726 129 11,214 +299 ⋆ C2 MIP 55,311,277 - - - - 0.83LMP - 13,242 35 - NTLP 55,291,978 3,187 15 10,055 NIA1 55,309,361 1,916 177 +1,271 YIA2 55,309,361 1,916 194 +1,271 YIAC1 55,309,361 1,916 ⋄ (+0) +1,271 YIAC2 55,309,361 1,916 ⋄ (+0) +1,271 YOPT 55,309,361 1,916 10,445 +1,271 ⋆ C3 MIP 69,299,295 - - - - 0.99LMP - 18,114 35 - NTLP 69,290,109 2,890 15 15,224 NIA1 69,297,643 1,652 240 +1,238 YIA2 69,297,517 1,778 132 +1,112 NIAC1 69,297,643 1,652 ⋄ (+0) +1,238 YIAC2 69,297,643 1,652 ⋄ (+160) +1,238 YOPT 69,297,643 1,652 18,975 +1,238 ⋆ C4 MIP 50,763,422 - - - - 0.52LMP - 8,948 34 - NTLP 50,754,612 2,033 18 6,915 NIA1 50,762,469 953 431 +1,080 YIA2 50,762,469 953 390 +1,080 YIAC1 50,762,469 953 ⋄ (+0) +1,080 YIAC2 50,762,469 953 ⋄ (+0) +1,080 YOPT 50,762,469 953 15,185 +1,080 ⋆ C5 MIP 49,946,355 - - - - 0.34LMP - 7,637 35 - NTLP 49,861,587 6,046 16 1,591 NIA1 49,945,716 639 493 +5,407 YIA2 49,945,716 639 769 +5,407 YIAC1 49,945,716 639 ⋄ (+0) +5,407 YIAC2 49,945,716 639 ⋄ (+0) +5,407 YOPT 49,945,716 639 16,260 +5,407 ⋆ C6 MIP 58,880,776 - - - - 0.76LMP - 36,630 38 - NTLP 58,861,963 3,889 15 32,741 NIA1 58,880,049 727 475 +3,162 YIA2 58,880,049 727 597 +3,162 YIAC1 58,880,049 727 ⋄ (+0) +3,162 YIAC2 58,880,049 727 ⋄ (+0) +3,162 YOPT 58,880,049 727 16,098 +3,162 ⋆ C7 MIP 57,307,363 - - - - 1.01LMP - 51,424 37 - NTLP 57,288,519 1,832 14 49,592 NIA1 57,306,102 1,261 584 +571 NIA2 57,306,120 1,243 560 +589 YIAC1 57,306,120 1,243 ⋄ (+38) +589 YIAC2 57,306,120 1,243 ⋄ (+0) +589 YOPT 57,306,120 1,243 11,935 +589 ⋆ C8 MIP 69,977,540 - - - - 0.37LMP - 11,472 35 - NTLP 69,950,201 5,230 17 6,242 NIA1 69,977,111 429 758 +4,801 YIA2 69,977,111 429 780 +4,801 YIAC1 69,977,111 429 ⋄ (+0) +4,801 YIAC2 69,977,111 429 ⋄ (+0) +4,801 YOPT 69,977,111 429 10,993 +4,801 ⋆ C9 MIP 47,889,206 - - - - 0.60LMP - 8,875 38 - NTLP 47,860,497 4,042 14 4,833 NIA1 47,887,537 1,669 188 +2,373 YIA2 47,887,537 1,669 234 +2,373 YIAC1 47,887,537 1,669 ⋄ (+0) +2,373 YIAC2 47,887,537 1,669 ⋄ (+0) +2,373 YOPT 47,887,537 1,669 10,363 +2,373 ⋆ C10 MIP 59,195,531 - - - - 0.68LMP - 11,613 36 - NTLP 59,193,235 1,899 13 9,714 NIA1 59,194,229 1,302 108 +597 YIA2 59,194,229 1,302 115 +597 YIAC1 59,194,229 1,302 ⋄ (+0) +597 YIAC2 59,194,229 1,302 ⋄ (+0) +597 YOPT 59,194,229 1,302 9,584 +597 ⋆ C11 MIP 49,628,808 - - - - 0.38LMP - 9,918 38 - NTLP 49,620,385 1,448 17 8,470 NIA1 49,627,991 817 372 +631 YIA2 49,627,991 817 115 +631 YIAC1 49,627,991 817 ⋄ (+0) +631 YIAC2 49,627,991 817 ⋄ (+0) +631 YOPT 49,627,991 817 16,269 +631 ⋆ TABLE IIT
EST RESULTS FOR
MISO
WITH TRANSMISSION CONSTRAINTS
Case Model Solution($) Uplift Payment($) Time(s) Save($) Optimal Diff($/MWh)C10(T) MIP 61,717,153 - 584 - - 3.49LMP - 1,667,967 68 - NTLP 61,596,521 92,541 69 1,575,426 NIA1 61,602,290 87,824 1,182 +4,717 YIA2 61,602,290 87,824 1,240 +4,717 YIAC1 61,602,290 87,824 ⋄ (+0) +4,717 YIAC2 61,602,290 87,824 ⋄ (+0) +4,717 YOPT 61,602,290 87,824 81,630 +4,717 ⋆ C11(T) MIP 50,071,094 - 271 - - 2.19LMP - 476,190 58 - NTLP 50,020,529 24,538 41 451,652 NIA1 50,030,415 23,498 512 +1,041 NIA2 50,030,417 23,495 626 +1,044 YIAC1 50,030,417 23,495 ⋄ (+39) +1,044 YIAC2 50,030,417 23,495 ⋄ (+0) +1,044 YOPT 50,030,417 23,495 31,857 +1,044 ⋆ physical and operational restrictions. More importantly, totackle the computational complexity issues brought by thelarge-scale MISO system, iterative algorithms with conver-gence properties have been proposed. The iterative algorithmsinnovatively take advantage of the uplift payment calculationprocess and the cutting plane technology, which can achieve ahigh-quality solution in a short time. The test results verifiedthe efficiency of our method in both cost savings and solvingtime reduction. In the next step, we will explore other waysto solve this large-scale convex hull pricing problem, such asdecomposition approaches described in [23].A CKNOWLEDGEMENTS
The authors thank the editor and referees for their sinceresuggestions on improving the quality of this paper.A
PPENDIX AP ROOF OF T HEOREM Proof.
In this proof, we show that EUC formulation is the dualformulation of a primal dynamic programming formulation(similar to [24]). Here, we let T = |T | . First, for the dynamicprogramming algorithm, we define two value functions V ↓ ( t ) and V ↑ ( t ) for each time period, in which V ↓ ( t ) represents thetotal net cost (i.e., the total cost minus the revenue) from time t to the end when the generator shuts down at time t + 1 , and V ↑ ( t ) represents the total net cost from time t to the end whenthe generator starts up at time t . Basically, we use V ↓ ( t ) and V ↑ ( t ) to keep track of the starting points of the “OFF” and“ON” intervals. Accordingly, we let C ( t, k ) represent the netcost for the “ON” period [ t, k ] Z , where k ≥ t . The dynamicprogramming equations are as follows: V ↓ ( t ) = min k ∈ [ t + ℓ t +1 ,T ] Z { S ′ ( k − t −
1) + V ↑ ( k ) , } , ∀ t ∈ [ t , T − ℓ t − Z , (29a) V ↓ ( t ) = 0 , ∀ t ∈ [ T − ℓ t , T − Z , (29b) V ↑ ( t ) = min k ∈ [ t + Lt − , min { t + L − ,T − } ] Z { S ( k − t +1)+ C ( t, k )+ V ↓ ( k ) , ˆ C ( t, T ) } , ∀ t ∈ [ t + ℓ t + 1 , T − L t ] Z , (29c) V ↑ ( t ) = C ( t, T ) , ∀ t ∈ [max { T − L t +1 ,T − L +1 } ,T ] Z , (29d)where equations (29a) represent that when the generator shutsdown at time t + 1 , it can either start up again at time k with k − t − ≥ ℓ t (satisfying min-down time limit) and k ≤ T or keep offline until the end. Equations (29c) indicatethat when the generator starts up at time t , it can either keeponline until time k and shut-down at time k + 1 with L t ≤ k − t + 1 ≤ L (satisfying min-up and max-up time limit)and k ≤ T − or keep online until the end while satisfyingthe max-up time limit ( ˆ C ( t, T ) represents the net cost if t ≥ T − L + 1 , and + inf , otherwise). Equations (29b) and (29d)describe the backward initial conditions. The objective is Φ= V ↑ (0):=min t ∈ [ t , min { L − s ,T − } ] Z n S ( t + s )+ C (1 , t )+ V ↓ ( t ) , e C (1 , T ) o , (30)where e C (1 , T ) represents the net cost if L − s ≥ T , and + inf , otherwise.This dynamic program can be written as the followingequivalent formulation: max Φ (31a) ( w t ) s.t. Φ ≤ S ( t + s ) + C (1 , t ) + V ↓ ( t ) , ∀ t ∈ [ t , min { L − s , T − } ] Z , (31b) ( w t ) Φ ≤ C (1 , T ) , if L − s ≥ T, (31c) ( z kt ) V ↓ ( k ) ≤ S ′ ( t − k −
1) + V ↑ ( t ) , ∀ t ∈ [ k + ℓ t +1 , T ] Z , ∀ k ∈ [ t , T − ℓ t − Z , (31d) ( θ t ) V ↓ ( t ) ≤ , ∀ t ∈ [ t , T − ℓ t − Z , (31e) ( θ t ) V ↓ ( t ) = 0 , ∀ t ∈ [ T − ℓ t , T − Z , (31f) ( y tk ) V ↑ ( t ) ≤ S ( k − t + 1) + C ( t, k ) + V ↓ ( k ) , ∀ k ∈ [ t + L t − , min { t + L − , T − } ] Z , ∀ t ∈ [ t + ℓ t +1 , T − L t ] Z , (31g) ( y tT ) V ↑ ( t ) ≤ C ( t, T ) , ∀ t ∈ [max { t + ℓ t +1 ,T − L +1 } , T − L t ] Z , (31h) ( y tT ) V ↑ ( t )= C ( t,T ) , ∀ t ∈ [max { T − L t +1 ,T − L +1 } ,T ] Z , (31i)where in general C ( t, k ) = min k X s = t φ s (32a)s.t. P s ≤ x s ≤ P s , ∀ s ∈ [ t, k ] Z , x t ≤ V ut , x k ≤ V ek , (32b) x s − x s − ≤ V us , x s − − x s ≤ V es , ∀ s ∈ [ t + 1 , k ] Z , (32c) φ s ≥ a j x s + b j , ∀ s ∈ [ t, k ] Z , j ∈ [1 , N ] Z , (32d)where constraints (32b) represent the capacity and start-up/shut-down restrictions, constraints (32c) represent theramping constraints, and constraints (32d) indicate that theobjective function is represented by a piecewise linear functionwith N pieces.By taking the dual of model (32), we can remove the “max”sign and insert the dual formulation into model (31) and thenobtain an overall dual linear program corresponding to theoriginal dynamic program as follows with C ( t, k ) as a linkingdecision variable between (31) and the dual of (32): max Φ (33a)s.t. (31b) − (31i) , (33b) ( p tk ) C (1 , k ) ≤ M , ∀ k ∈ [ t +1 , min { L − s ,T − } ] Z , (33c) ( p tk ) C ( t, k ) ≤ M , ∀ k ∈ [ t + L t − , min { t + L − , T − } ] Z , i ∀ t ∈ [ t + ℓ t +1 , T − L t ] Z , (33d) ( p tk ) C ( t, T ) ≤ M , ∀ t ∈ [max { t + ℓ t +1 ,T − L +1 } ,T ] Z , (33e) ( p tk ) C (1 , T ) ≤ M , if L − s ≥ T, (33f)Constraints in the dual formulation of (32) , (33g)where M , M , M , M represent the dual objective functionof the economic dispatch problem under four different casesbased on the value of t and k , with the detailed representationomitted for brevity.By taking the dual of the derived linear program (33)and cleaning the model, we can obtain model (15) (EUCformulation) as described in Section IV-B2.To prove that D i (i.e., the feasible region for formula-tion (15)) is a convex hull description for single UC, we needto show for any arbitrary linear objective function with it asa feasible region (denoted as formulation LD ), there exists acorresponding optimal solution binary with respect to decisionvariables w , y , and z . To show this, based on the optimal solu-tion obtained from the dynamic programming framework, wecan construct a corresponding solution ( w ∗ , y ∗ , z ∗ , θ ∗ , q ∗ , φ ∗ ) for LD , where w ∗ represents the first shut-down, y ∗ repre-sents the “ON” interval, z ∗ represents the “OFF” interval, θ ∗ represents the final shut-down, and q ∗ and φ ∗ representthe generation amount and cost in the “ON” interval, inthe optimal schedule provided by the dynamic programmingapproach. It is easy to check ( w ∗ , y ∗ , z ∗ , θ ∗ , q ∗ , φ ∗ ) is feasiblefor LD since all physical constraints related to generationamount are derived from the economic dispatch problem (32)and the constraints related to binary variables are derivedfrom the dynamic programming framework. Meanwhile, theobjective value of the constructed solution to LD is equal tothe optimal value of the dynamic programming V ↑ (0) . Then,by the strong duality theorem, ( w ∗ , y ∗ , z ∗ , θ ∗ , q ∗ , φ ∗ ) , binarywith respect to decision variables w , z , and y , is an optimalsolution to LD . Thus, the conclusion holds. Also, there existsan optimal solution to (15) binary with respect to decisionvariables w, y, and z , since the objective function is linear andthere exists an optimal solution at the extreme point.R EFERENCES[1] Y. Yu, Y. Guan, and Y. Chen, “An integral formulation and convex hullpricing for unit commitment,”
Tech. Rep., University of Florida , 2019.[Online]. Available: https://arxiv.org/abs/1906.07862[2] C. Wang, P. B. Luh, P. Gribik, T. Peng, and L. Zhang, “Commitmentcost allocation of fast-start units for approximate extended locationalmarginal prices,”
IEEE Trans. Power Syst. , vol. 31, no. 6, pp. 4176–4184, Nov. 2016.[3] F. S. Bresler, “Energy and ancillary services uplift in PJM,”
FERC UpliftWorkshop , 2014.[4] W. Sauer, “Uplift in RTO and ISO markets,” Fed. Energy Regul.Commission, Washington, DC, USA, Tech. Rep., 2014.[5] P. Gribik, W. Hogan, and S. Pope, “Market-clearing electricity pricesand energy uplift,”
Cambridge, MA , Dec. 2007.[6] D. A. Schiro, T. Zheng, F. Zhao, and E. Litvinov, “Convex hullpricing in electricity markets: Formulation, analysis, and implementationchallenges,”
IEEE Trans. Power Syst. , vol. 31, no. 5, pp. 4068–4075,Sep. 2016.[7] B. Hua and R. Baldick, “A convex primal formulation for convex hullpricing,”
IEEE Trans. Power Syst. , vol. 32, no. 5, pp. 3814–3823, Sep.2017.[8] Z. Yang, T. Zheng, J. Yu, and K. Xie, “A unified approach to pricingunder nonconvexity,”
IEEE Trans. Power Syst. , vol. 34, no. 5, pp. 3417–3427, Sep. 2019. [9] G. Wang, U. Shanbhag, T. Zheng, E. Litvinov, and S. Meyn, “Anextreme-point subdifferential method for convex hull pricing in energyand reserve markets – part I: Algorithm structure,”
IEEE Trans. PowerSyst. , vol. 28, no. 3, pp. 2111–2120, Aug. 2013.[10] G. Wang, U. V. Shanbhag, T. Zheng, E. Litvinov, and S. Meyn, “Anextreme-point subdifferential method for convex hull pricing in energyand reserve markets – part II: Convergence analysis and numericalperformance,”
IEEE Trans. Power Syst. , vol. 28, no. 3, pp. 2121–2127,Aug. 2013.[11] C. Wang, T. Peng, P. B. Luh, P. Gribik, and L. Zhang, “The subgradientsimplex cutting plane method for extended locational marginal prices,”
IEEE Trans. Power Syst. , vol. 28, no. 3, pp. 2758–2767, Aug. 2013.[12] G. Morales-Espaa, J. M. Latorre, and A. Ramos, “Tight and compactMILP formulation for the thermal unit commitment problem,”
IEEETrans. Power Syst. , vol. 28, no. 4, pp. 4897–4908, Nov. 2013.[13] A. Frangioni and C. Gentile, “An extended MIP formulation for thesingle-unit commitment problem with ramping constraints,” , 2015.[14] B. Knueven, J. Ostrowski, and J. Wang, “The ramping polytope and cutgeneration for the unit commitment problem,”
INFORMS J. on Comput. ,vol. 30, no. 4, pp. 739–749, Nov. 2018.[15] MISO, “ELMP III White Paper I,” Jan. 2019.[16] Y. Chen and C. Wang, “Enhancements of extended locational marginalpricing - advancing practical implementation,”
CIGRE Grid of theFuture, Atlanta
IEEE Trans. Power Syst. , vol. 31, no. 6, pp. 4732–4743,Nov. 2016.[19] MISO, “Production cost model fundamentals,” 2008.[20] A. J. Wood and B. F. Wollenberg,
Power Generation, Operation, andControl . Hoboken, NJ, USA: Wiley-Interscience, 2014.[21] C. Gentile, G. Morales-Espa˜na, and A. Ramos, “A tight MIP formu-lation of the unit commitment problem with start-up and shut-downconstraints,”
Eur. J. on Comput. Optim. , vol. 5, no. 1, pp. 177–201,Mar. 2017.[22] M. Queyranne and L. A. Wolsey, “Tight MIP formulations for boundedup/down times and interval-dependent start-ups,”
Math. Program. , vol.164, no. 1, pp. 129–155, Jul. 2017.[23] F. Zhao, E. Litvinov, and T. Zheng, “A marginal equivalent decom-position method and its application to multi-area optimal power flowproblems,”
IEEE Trans. Power Syst. , vol. 29, no. 1, pp. 53–61, Jan.2014.[24] Y. Guan, K. Pan, and K. Zhou, “Polynomial time algorithms andextended formulations for unit commitment problems,”