Robust Optimal Operation of Virtual Power Plants Under Decision-Dependent Uncertainty of Price Elasticity

JOURNAL OF L A T E X CLASS FILES, VOL. XX, NO. X, FEB. 2025 1 Rob ust Optimal Operation of V irtual Po wer Plants Under Decision-Dependent Uncertainty of Price Elasticity T ao T an, Rui Xie, Member , IEEE , Meng Y ang, and Y ue Chen, Senior Member , IEEE Abstract —The rapid deployment of distrib uted energy re- sources (DERs) is one of the essential efforts to mitigate global climate change. Howev er , a vast number of small-scale DERs are difficult to manage individually , motivating the introduction of virtual power plants (VPPs). A VPP operator coordinates a group of DERs by setting suitable prices, and aggregates them for interaction with the power grid. In this context, optimal pricing plays a critical role in VPP operation. This paper pr oposes a rob ust optimal operation model for VPPs that considers uncertainty in the price elasticity of demand. Specifically , the demand elasticity is found to be influenced by the pricing decision, giving rise to decision-dependent uncertainty (DDU). An impro ved column-and-constraint (C&CG) algorithm, together with tailored transformation and ref ormulation techniques, is developed to solve the rob ust model with DDU efficiently . Case studies based on actual electricity consumption data of London households demonstrate the effectiveness of the proposed model and algorithm. Index T erms —r obust optimization, virtual power plant, decision-dependent uncertainty , optimal pricing, price elasticity N O M E N C L AT U R E Acr onym C&CG Column-and-constraint generation. DDU Decision-dependent uncertainty . DER Distributed energy resource. DR O Distributionally robust optimization. R O Robust optimization. SP Stochastic programming. TOU T ime-of-use. VPP V irtual power plant. Indices and Sets I Set of power nodes in VPP . L Set of power lines in VPP . T Set of periods. K Set of intervals of TOU price. Corresponding author: Rui Xie T . T an, M. Y ang, and Y . Chen are with the Department of Mechani- cal and Automation Engineering, The Chinese University of Hong K ong, HKSAR, China (e-mail: ttan@mae.cuhk.edu.hk, myang@mae.cuhk.edu.hk, yuechen@mae.cuhk.edu.hk). R. Xie is with the Department of Mechanical and Automation Engineering, The Chinese Univ ersity of Hong K ong, HKSAR, China, and also with the College of Electrical and Information Engineering, Hunan University , Changsha 410082, China (e-mail: xierui@hnu.edu.cn). P arameters L j t Predicted load at node j in period t . κ j Ratio of reactive demand to active demand at node j . R ij , X ij Resistance and reactance of line i → j . P i , P i Lower and upper bounds of acti ve power gener- ation at node i . Q i , Q i Lower and upper bounds of reactive power gen- eration at node i . S ij Maximum apparent power flow on line i → j . V i , V i Lower and upper bounds of the voltage magni- tude at node i . C T OU Lower bound of the TOU price. C T OU Upper bound of the TOU price. C RE F t Reference price in period t . r − tk , r + tk Lower and upper bounds of the ratio of TOU price to the reference price in period t and interval k . ρ t Day-ahead market price in period t . ρ + t , ρ − t Real-time purchasing/selling price from/to the power grid in period t . ρ G i Unit generation cost of activ e power at node i . ξ − itk , ξ + itk Lower and upper bounds of price elasticity at node i in period t and interval k . V ariables F irst-Stag e (Day-Ahead) V ariables: p DA it , q DA it Scheduled active and reactive power generation at node i in period t . p DA ij t , q DA ij t Scheduled inflo w activ e and reactiv e power on line i → j in period t . v DA it Square of scheduled voltage magnitude at node i in period t . c T OU t TOU price in period t . p DA 0 t Scheduled active power injection at node 0 in period t . Uncertainty V ariables: ξ it Price elasticity at node i in period t . Second-Stage (Real-T ime) V ariables: l it Real-time demand at node i in period t . p RT it , q RT it Real-time active and reactive po wer at node i in period t . p RT ij t , q RT ij t Real-time inflo w active and reactiv e power on line i → j in period t . JOURNAL OF L A T E X CLASS FILES, VOL. XX, NO. X, FEB. 2025 2 v RT it Square of real-time voltage magnitude at node i in period t . p ∆+ 0 t , p ∆ − 0 t Upward and do wnward power adjustment at node 0 in period t . I . I N T RO D U C T I O N Distributed energy resources (DERs), such as electric v e- hicles, battery energy storage, and rooftop solar panels, have been widely deplo yed to reduce dependence of power grids on fossil fuels and mitigate global warming [1]. In this context, virtual power plants (VPP) have been introduced, which ag- gregate multiple small-scale DERs through intelligent control to function collectively like traditional lar ge po wer plants [2]. During the operation of a VPP , the internal pricing strategies to incenti vize the participation of DERs are important, which need to consider the uncertainty of DER response [3]. T o deal with uncertainty in VPPs, three typical method- ologies are widely adopted, including stochastic programming (SP), robust optimization (R O), and distrib utionally robust optimization (DR O). SP models uncertain factors, e.g., re- new able power generation, as random variables with gi ven probability distributions, and deriv es the optimal solution based on sampling or chance constraints [4]. Although SP is straightforward, it requires exact and detailed probability distributions, which makes it less practical. Moreov er , it is computationally demanding to solve large-scale SP problems with numerous scenarios. RO aims to determine a solution that remains feasible and optimal under a range of possible scenarios [5]. It is easier to implement than SP without re- quiring detailed probabilistic information. DR O e xtends R O by optimizing the decision over a set of probability distributions rather than a set of scenarios [6]. Although DR O can obtain a less conservati ve result than R O, its formulation and solution are significantly more complex. More importantly , SP and DR O allow a certain lev el of risk [7], making it difficult to guarantee 100% confidence in the security of the VPP . For these reasons, RO is a more suitable method to optimize the operation of VPPs where security is of top priority . R O has been widely employed in VPPs. For example, the day-ahead ener gy storage scheduling within a VPP was studied in [8], considering the uncertainty of renewable generation and load. In [9], the v alues of rene wable-only VPPs and grid-scale energy storage systems in energy and reserve markets were compared using a two-stage RO framework. These studies mainly focus on the decision-independent uncertainty (DIU), in which the decisions in the first stage have no impact on the uncertainty set. R O with DIU can be solved efficiently using the Benders decomposition [10], the column-and-constraint generation (C&CG) algorithm [11], etc. These algorithms iterativ ely identify the worst-case scenarios and return the scenarios or the corresponding cutting planes to the first-stage problem until con vergence. Recently , decision-dependent uncertainty (DDU), where the first-stage decisions can affect the uncertainty set, has been widely recognized. R O with DDU is an emerging topic in power systems [12] and there are few studies considering DDU in VPPs. A stochastic adaptive R O model for a VPP’ s scheduling in day-ahead energy-reserve markets was estab- lished in [13], considering the DDU of real-time reserv e deployment requests. The DDU from energy storage operation was modeled in [14] and then integrated into a stochastic optimization method for VPP . DDU substantially adds difficulties to modeling RO prob- lems, because decisions can change the uncertainty set in various ways [15]. A typical example is demand response, a popular research topic in VPPs, in which studies encompass optimal scheduling [16], bidding strategy [17], and real-time pricing [18]. The responsiveness of electricity demand is influ- enced by time-of-use (TOU) pricing [19]. Such responsiveness is usually characterized by a concept called price elasticity [20]. In this paper , the uncertainty of price elasticity and how it is influenced by the pricing decision are considered. In terms of solving R O problems with DDU, traditional methods are found to be intractable, suboptimal, or lacking con vergence guarantees [21]. Therefore, innov ativ e modeling techniques and solution algorithms are required to address RO with DDU. Currently , sev eral variants of C&CG have been dev eloped for this problem. Adapti ve C&CG is an example, which returns the set of active constraints to represent the worst-case scenario identified [21]. Another algorithm is map- ping C&CG, which applies mapping rules to maintain the worst-case scenarios at the vertices of the new uncertainty set [22]. Dual C&CG and transformation-based C&CG apply duality and variable substitution to conv ert an RO with DDU problem to an RO with DIU problem [23], [24]. There are also some other algorithms, including variants of Benders decomposition [25], multi-parametric programming [26], K- adaptability [27], and generic solution algorithms [28]. Ho w- ev er , these general algorithms might be time-consuming and difficult to implement due to the complicated transformation and projection. This paper aims to fill the research gaps above by proposing a robust VPP optimal operation model considering the impact of pricing strategies on the uncertain demand elasticity , and dev eloping a novel and tractable solution algorithm based on the special structure of the problem studied. Our main contribution is two-fold: 1) Robust VPP Optimal Operation Model Under DDU . W e propose a two-stage robust optimal operation model for VPPs to determine their internal T OU tarif fs for mini- mizing the ov erall cost. Particularly , the impact of first- stage pricing decisions on the uncertain price elasticity is considered, which renders DDU. Compared to tradi- tional robust VPP operation models with DIU only , the proposed model accounts for the mutual impact between decisions and uncertainty , and thus gi ves more precise optimal operation and pricing decisions. 2) Impr oved C&CG Algorithm . T o solve the robust optimal operation problem of a VPP with DDU, we dev elop an improv ed C&CG algorithm that leverages an uncertainty- set transformation. Compared with standard R O solution methods that may suf fer from non-con ver gence or subop- timality , the proposed algorithm is guaranteed to con verge to a globally optimal solution. The remainder of the paper is organized as follows. Sec- JOURNAL OF L A T E X CLASS FILES, VOL. XX, NO. X, FEB. 2025 3 tion II introduces the robust optimal operation model for VPPs considering DDU. The solution methodology is presented in Section III. Numerical case studies are conducted in Sec- tion IV. Finally , conclusions are summarized in Section V. I I . R O B U S T O P T I M A L O P E R A T I O N M O D E L O F V I RT UA L P OW E R P L A N T S In this section, we propose a robust optimal operation model for a VPP , considering the impact of pricing decisions on the demand uncertainty . The ov erview of the VPP system is shown in Fig. 1. W e first introduce the first-stage constraints, the decision-dependent uncertainty set, and the second-stage constraints, respectiv ely . Then, we gi ve the ov erall robust optimal operation model. Renewable energy Controllable generators Industrial load Commercial load Residential load Upstream grid Energy flow Cash flow Day-ahead price Real-time price Price Load Impact Decision- dependent VPP system               Fig. 1. The overvie w of the VPP system. A. F irst-Stage Constraints In the first stage, the VPP operator determines the day- ahead setpoints for generators and the TOU prices for demand, subject to power network constraints. Since the VPP resides at the distrib ution lev el, the LinDistFlo w model is adopted. Particularly , the first-stage operational constraints include the following: p DA ij t + p DA j t − L j t = X l : j → l p DA j lt , ∀ i → j ∈ L , ∀ t ∈ T , (1a) q DA ij t + q DA j t − κ j L j t = X l : j → l q DA j lt , ∀ i → j ∈ L , ∀ t ∈ T , (1b) v DA j t = v DA it − 2 R ij p DA ij t − 2 X ij q DA ij t , ∀ i → j ∈ L , ∀ t ∈ T , (1c) P i ≤ p DA it ≤ P i , Q i ≤ q DA it ≤ Q i , ∀ i ∈ I , ∀ t ∈ T , (1d) ( p DA ij t ) 2 + ( q DA ij t ) 2 ≤ ( S ij ) 2 , ∀ i → j ∈ L , ∀ t ∈ T , (1e) ( V i ) 2 ≤ v DA it ≤ ( V i ) 2 , ∀ i ∈ I , ∀ t ∈ T , (1f) C T OU ≤ c T OU t ≤ C T OU , ∀ t ∈ T . (1g) Constraints (1a) and (1b) describe the activ e and reactive power balance at node j ∈ I in period t ∈ T . Here, p DA j t and q DA j t are the scheduled activ e and reactiv e power generation, p DA ij t and q DA ij t are the scheduled inflo w acti ve and reactive power on line i → j , L j t is the predicted activ e load, and κ j is the ratio of reacti ve to acti ve demand. Constraint (1c) models the voltage drop along power line i → j , where v DA it is the square of the voltage magnitude at node i , and R ij and X ij are the line resistance and reactance, respecti vely . Constraint (1d) enforces the generation limits for acti ve and reactiv e power at node i . Constraint (1e) limits the apparent power flow on line i → j to its maximum rating S ij . Constraint (1f) ensures that the square of the voltage magnitude at each node remains within its secure range. Finally , constraint (1g) bounds the TOU price c T OU t . B. Decision-Dependent Uncertainty Set In this paper , we focus on the demand uncertainty; or more specifically , the uncertain price elasticity ξ it indicating ho w demand changes with prices. Based on analysis of real-world data (gi ven in Section IV later), we find that the variation range of ξ it depends on the ratio between the determined TOU price c T OU t and the reference price C RE F t . Therefore, the uncertainty set is formulated as a function of the first- stage decision vector c T OU = ( c T OU t ) t ∈T . W e denote this set as U ( c T OU ) : U ( c T OU ) =      ξ = ( ξ it ) i ∈I ,t ∈T        ∀ i ∈ I , ∀ t ∈ T , ∀ k ∈ K : if r − tk ≤ c T OU t C RE F t ≤ r + tk , then ξ − itk ≤ ξ it ≤ ξ + itk      . (2) In (2), the set of interv als K partitions the possible price ratios. For interval k ∈ K defined by bounds [ r − tk , r + tk ] , there is a corresponding set of elasticity bounds [ ξ − itk , ξ + itk ] for the uncertainty variable ξ it . This uncertainty set is essentially a decision-dependent uncertainty set because the first-stage decision c T OU t influences the range of ξ it . This decision dependency challenges the solution of the RO problem, as rev ealed by [21]. In Section III, we propose a reformulation of the uncertainty set to address this challenge. C. Second-Stage Constraints In the second stage, after observing the actual demand, the VPP operator adjusts the po wer output of generators to maintain real-time balance between electricity supply and demand. The second-stage constraints include the following: l it = L it (1 + ξ it ( c T OU t /C RE F t − 1)) , ∀ i ∈ I , ∀ t ∈ T , (3a) p DA 0 t + p ∆+ 0 t − p ∆ − 0 t = p RT 0 t , ∀ t ∈ T , (3b) p ∆+ 0 t ≥ 0 , p ∆ − 0 t ≥ 0 , ∀ t ∈ T , (3c) p RT ij t + p RT j t − l j t = X l : j → l p RT j lt , ∀ i → j ∈ L , ∀ t ∈ T , (3d) q RT ij t + q RT j t − κ j l j t = X l : j → l q RT j lt , ∀ i → j ∈ L , ∀ t ∈ T , (3e) v RT j t = v RT it − 2 R ij p RT ij t − 2 X ij q RT ij t , ∀ i → j ∈ L , ∀ t ∈ T , (3f) P i ≤ p RT it ≤ P i , Q i ≤ q RT it ≤ Q i , ∀ i ∈ I , ∀ t ∈ T , (3g) ( p RT ij t ) 2 + ( q RT ij t ) 2 ≤ ( S ij ) 2 , ∀ i → j ∈ L , ∀ t ∈ T , (3h) ( V i ) 2 ≤ v RT it ≤ ( V i ) 2 , ∀ i ∈ I , ∀ t ∈ T . (3i) Constraint (3a) defines the realized demand l it based on the predicted load L it , the TOU price c T OU t , and the realized JOURNAL OF L A T E X CLASS FILES, VOL. XX, NO. X, FEB. 2025 4 uncertainty ξ it . Constraints (3b) and (3c) model the real- time po wer p RT 0 t purchased from the upstream grid. This is composed of the day-ahead schedule p DA 0 t and upward p ∆+ 0 t or downw ard p ∆ − 0 t adjustments. Constraints (3d) and (3e) enforce the real-time activ e and reactiv e power balance, respecti vely . The remaining constraints (3f)–(3i) are the real-time network operational constraints. They model the voltage drop, enforce generation limits, bound the line apparent power , and ensure that the square of the voltage magnitude remains within its secure range. D. Overall Rob ust Optimal Oper ation Model The overall robust optimal operation model for the VPP is formulated as a two-stage R O problem. W e use notation x to denote the collection of the first-stage decision v ariables p DA it , q DA it , p DA ij t , q DA ij t , v DA it , and c T OU t , and y for the collection of second-stage variables l it , p RT it , q RT it , p RT ij t , q RT ij t , v RT it , p ∆+ 0 t , and p ∆ − 0 t . The model is as follows: min x ∈X ( X t ∈T ρ t p DA 0 t + max ξ ∈U ( c T OU ) min y ∈Y ( x,ξ ) ( − X t ∈T X i ∈I c T OU t l it + X t ∈T  ρ + t p ∆+ 0 t − ρ − t p ∆ − 0 t  + X t ∈T X i ∈I ρ G i p RT it )) , (4) where the uncertainty set U ( c T OU ) is defined in (2), and the first- and second-stage feasible regions ( X and Y ( x, ξ ) ) are: X =  x = ( p DA it , q DA it , p DA ij t , q DA ij t , v DA it , c T OU t )   (1)  , (5a) Y ( x, ξ ) =  y = ( l it , p RT it , q RT it , p RT ij t , q RT ij t , v RT it , p ∆+ 0 t , p ∆ − 0 t )   (3)  . (5b) The objectiv e function in (4) is composed of two parts. The first term P t ∈T ρ t p DA 0 t is the first-stage (day-ahead) cost of purchasing energy from the upstream grid at day-ahead price ρ t . The second part is the worst-case second-stage (real- time) operational cost. The max operator finds the worst- case scenario of the price elasticity ξ from the decision- dependent uncertainty set U ( c T OU ) . The inner min operator represents the VPP’ s real-time adjustments to minimize total operational costs, including the negati ve income of selling electricity to clients P t ∈T P i ∈I c T OU t l it , the adjustment cost P t ∈T ( ρ + t p ∆+ 0 t − ρ − t p ∆ − 0 t ) for purchasing/selling power from/to the upstream grid based on the real-time prices, and the actual power generation cost P t ∈T P i ∈I ρ G i p RT it . I I I . S O L U T I O N M E T H O D Solving the robust VPP optimal operation model (4) faces two challenges: (i) The nonlinear constraint (3h) hinders the transformation of the second-stage problem based on optimal- ity conditions. (ii) The decision-dependent uncertainty set (2) makes traditional RO algorithms inapplicable. A. Linearization of Quadratic Constraints W e first address the challenge (i). The constraint ( p RT ij t ) 2 + ( q RT ij t ) 2 ≤ ( S ij ) 2 represents a disk region on the ( p RT ij t , q RT ij t ) plane, which can be approximated by its inscribed regular polygons with S sides. Therefore, we use the following constraints to inner approximate ( p RT ij t ) 2 + ( q RT ij t ) 2 ≤ ( S ij ) 2 : p RT ij t cos  2 s − 1 S π  + q RT ij t sin  2 s − 1 S π  ≤ cos  π S  S ij , ∀ s = 1 , 2 , . . . , S. (6) Fig. 2 illustrates the basic idea of this inner approximation. As we can see, the larger the S , the more accurate the approximation. After this transformation, the second-stage problem turns out to be a linear program (LP), which facilitates the subsequent transformation. Fig. 2. The basic idea of inner approximation. B. Uncertainty Set Reformulation Next, to address challenge (ii), we propose an equi valent transformation (7) of the uncertainty set (2). The binary indi- cators z tk identify the interv al [ r − tk , r + tk ] into which the TOU- to-reference price ratio c T OU t /c RE F falls and accordingly determine the corresponding price elasticity ξ it . X k ∈K ξ − itk z tk ≤ ξ it ≤ X k ∈K ξ + itk z tk , ∀ i ∈ I , ∀ t ∈ T , (7a) X k ∈K r − tk z tk ≤ c T OU t c RE F ≤ X k ∈K r + tk z tk , ∀ t ∈ T , (7b) X k ∈K z tk = 1 , ∀ t ∈ T , (7c) z tk ∈ { 0 , 1 } , ∀ k ∈ K , ∀ t ∈ T . (7d) Specifically , constraints (7c) and (7d) ensure that for each time period t there is only one z tk = 1 while all others equal 0. Suppose z t ¯ k = 1 , then constraint (7b) gi ves r − t ¯ k ≤ c T OU t c RE F ≤ r + t ¯ k , and correspondingly (7a) gives ξ − it ¯ k ≤ ξ it ≤ ξ + it ¯ k , which is consistent with (2). JOURNAL OF L A T E X CLASS FILES, VOL. XX, NO. X, FEB. 2025 5 C. Impr oved C&CG Algorithm After the abov e transformation and reformulation, the robust optimal operation model (4) can be restructured in a compact form as shown in (8): min x,z C ⊤ x + max ξ min y E ( x ) ⊤ y , (8a) s.t. Ax ≥ B , (7) , P y ≥ Q ( x ) ξ + R ( x ) , (8b) where x and y collect all the first- and second-stage decision variables, respectively; A , B , C , and P are the corresponding coefficient matrices; the coefficient matrices E ( x ) , Q ( x ) , and R ( x ) are influenced by the first-stage decision x . In the following, we introduce our improv ed C&CG algorithm to solve the problem, which iterati vely solv es a master problem to update the first-stage decision x and solv es a feasibility check problem or a subproblem to identify the worst-case scenario under giv en x . 1) F easibility Chec k Pr oblem and Subpr oblem: First, given the first-stage decision x , we formulate the feasibility check problem (or subproblem) to identify the worst-case scenario that causes infeasibility (or the highest cost). The feasibility check problem is formulated as follows: FC : max ξ min y ,s 1 ⊤ s, (9a) s.t. (7) , P y + s ≥ Q ( x ) ξ + R ( x ) , s ≥ 0 . (9b) The FC problem is always feasible since y = 0 and s = max { 0 , Q ( x ) ξ + R ( x ) } is always a feasible solution. Moreov er, if the optimal solution to (9) satisfies s ∗ = 0 , then for e very ξ in the decision-dependent uncertainty set (7), there exists a feasible solution satisfying the original constraint P y ≥ Q ( x ) ξ + R ( x ) . T o solve the feasibility check problem, we turn the inner minimization problem into its Karush-Kuhn-T ucker (KKT) optimality condition: 1 − π − θ = 0 , (10a) P ⊤ π = 0 , (10b) 0 ≤ π ⊥ ( P y + s − Q ( x ) ξ − R ( x )) ≥ 0 , (10c) 0 ≤ θ ⊥ s ≥ 0 , (10d) where π and θ are the dual variables. Substituting (10) into the FC problem and linearizing the complementary slackness conditions (10c)-(10d) using the Big-M method, we get a mixed-inte ger linear program (MILP) that can be solved efficiently by commercial solvers. If the first-stage decision x passes the feasibility check, i.e., the optimal solution of (9) is s ∗ = 0 , then we move on to solve the subproblem: SP : max ξ min y E ( x ) ⊤ y , (11a) s.t. (7) , P y ≥ Q ( x ) ξ + R ( x ) . (11b) Similarly , we turn the inner minimization problem into its KKT condition: P ⊤ π = E ( x ) , (12a) 0 ≤ π ⊥ ( P y − Q ( x ) ξ − R ( x )) ≥ 0 . (12b) Again, substituting (12) into the SP problem (11), we get an MILP . Solving the FC or SP problem, we can obtain a worst-case scenario ξ ∗ . 2) Repr esentation of W orst-Case Scenario: Conv entionally , we directly return the worst-case scenario ξ ∗ to the master problem. Ho wever , since (7) is a decision-dependent uncer- tainty set, it has been re vealed that this con ventional method may lead to ov er-conserv ativeness or failure of con ver gence [12]. T o address this issue, we propose an innovati ve represen- tation of the worst-case scenario: According to [21], ξ always resides at a vertex of the uncertainty set (7). Therefore, it can be represented by ξ it = X k ∈K z tk  (1 − v it ) ξ − itk + v it ξ + itk  , (13) where v it is a binary variable. Moreover , to reduce conserva- tiv eness, we add the follo wing constraints to the uncertainty set, where Γ T and Γ S are uncertainty budgets. X i | 2 v it − 1 | ≤ Γ S , ∀ t, (14a) X t | 2 v it − 1 | ≤ Γ T , ∀ i. (14b) For each worst-case scenario ξ ∗ , we can determine the corre- sponding v ∗ based on (13). Then, instead of passing ξ ∗ directly to the master problem, we return the following constraints: ξ itm = X k ∈K z tk  (1 − v ∗ itm ) ξ − itk + v ∗ itm ξ + itk  , ∀ i ∈ I , ∀ t ∈ T , (15a) P y m ≥ Q ( x ) ξ m + R ( x ) . (15b) where the inde x m indicates the worst-case scenario in the m -th iteration. W ith (15), when the first-stage decision c T OU t changes, the z tk will change accordingly based on (7b)–(7d). Thus, the scenario ξ m = ( ξ itm ) i ∈I ,t ∈T will be mapped to the verte x of the active respective range in (2) based on (15a). 3) Master Pr oblem: Gi ven a set of v ectors V = { v ∗ m = ( v ∗ itm ; i ∈ I , t ∈ T ) ∈ [0 , 1] I T : m ∈ M} returned from previous iterations from FC or SP , the master problem is formulated by: MP : min x,z ,η,ξ ,y C ⊤ x + η , (16a) s.t. Ax ≥ B , (7b)–(7d) , (16b) η ≥ E ( x ) ⊤ y m , ∀ m ∈ M , (16c) (15) , ∀ m ∈ M . (16d) Howe ver , the master problem cannot be solved directly due to the nonlinear constraint (16c) with a bilinear term E ( x ) ⊤ y m ; more specifically , the term c T OU t l it . According to (3a) and (15a), we hav e l itm = L it + X k ∈K z tk ϕ km ( c T OU t /c RE F − 1) , (17) where ϕ km = L it  (1 − v ∗ itm ) ξ − itk + v ∗ itm ξ + itk  . (18) Therefore, the term c T OU t l it can be equiv alently represented by an expression of the term ( c T OU t ) 2 z tk , where z tk is a binary variable. This can be turned into solvable constraints JOURNAL OF L A T E X CLASS FILES, VOL. XX, NO. X, FEB. 2025 6 by replacing the term ( c T OU t ) 2 z tk in the objectiv e function by ω tk and adding the following constraints with a large enough constant M : ω tk ≥ 0 , (19a) ( c T OU t ) 2 ≤ ω tk + M (1 − z tk ) . (19b) When z tk = 0 , we have ω tk ≥ ( c T OU t ) 2 , so minimization ends up at ω tk = ( c T OU t ) 2 . When z tk = 1 , we hav e ω tk ≥ 0 and ( c T OU t ) 2 ≤ ω tk + M . Since M is a lar ge enough constant, the latter constraint always holds. Therefore, minimization ends up at ω tk = 0 . 4) Overall algorithm: The overall procedure of the pro- posed improv ed C&CG algorithm is presented in Algorithm 1. The con vergence and optimality of Algorithm 1 are stated in Proposition 1, whose proof is giv en in Appendix A. Algorithm 1 : Improved C&CG Algorithm 1: Initiation : Error tolerance τ > 0 ; N = 1 ; Let V be an empty set and assign a large value to U B 0 . 2: Solve the MP problem (16), using the technique in (19) to address the term ( c T OU t ) 2 z tk . Let x N ∗ be the optimal solution and LB N be the optimal value. 3: Solve the FC problem (9) with x N ∗ . If the optimal solution s ∗ = 0 , go to step 4; otherwise, let ξ N ∗ be the worst-case scenario and derive the corresponding v N ∗ by (13), let U B N = U B N − 1 , and go to Step 5. 4: Solve the SP problem (11) with x N ∗ . Let ( ξ N ∗ , y N ∗ ) be the optimal solution and U B N = C ⊤ x N ∗ + E ( x N ∗ ) ⊤ y N ∗ , and deriv e the corresponding v N ∗ by (13). Go to Step 5. 5: if | U B N − LB N | ≤ τ , terminate and output x N ∗ . Otherwise, add v N ∗ into V and let N = N + 1 ; go to Step 2. Pr oposition 1: Let N = 2 I T + 1 , where I is the number of nodes and T is the number of periods. If error tolerance τ = 0 , then Algorithm 1 con ver ges within N iterations and outputs the optimal solution of problem (8). I V . N U M E R I C A L E X P E R I M E N T S In this section, we test the performance of the proposed model and solution method using the IEEE-33 bus system, where 3 additional generators are connected to nodes 2, 3, and 6, respecti vely . All the simulations are conducted on a laptop with Intel(R) Core(TM) Ultra 9 185H 2.50 GHz processor and 32 GB RAM, using MA TLAB R2024A with GUR OBI 11.03. A. Simulation Settings First, we use the real-world energy consumption data from the Low Carbon London project [29] to form the uncertainty set. The dataset contains millions pieces of data. Ho wev er , the amount of useful data is much less than the amount of the entire full dataset, because most of the users were using fixed electricity and during the period of data, T oU was only applied in 2013. For these reasons, the selected data only accounts for a little of the entire data. Also, the referenced electricity price is calculated by a verage real-time electricity price in 2024 Hi gh - se n sitive us e r s - 1 - 0 . 9 - 0 . 8 - 0 . 7 - 0 . 6 - 0 . 5 - 0 . 4 - 0 . 3 - 0 . 2 - 0 . 1 El e c t r icity Pr ice E lasticit y 0 5 10 15 20 25 30 35 M id - se n sitive us e r s - 0. 18 - 0. 16 - 0. 14 - 0. 12 - 0. 1 - 0. 08 - 0. 06 - 0. 04 - 0. 02 0 0 5 10 15 20 25 30 35 40 L ow - se n sitive us e r s - 20 - 15 - 10 - 5 0 10 - 3 0 10 20 30 40 50 60 70 El e c t r icity Pr ice E lasticit y El e c t r icity Pr ice E lasticit y Number of u se r s Number of u se r s Number of u se r s Fig. 3. Up: High-sensiti ve user distribution at electricity price ratio 0.3; Middle: Medium-sensiti ve user distribution at electricity price ratio 3; Down: Low-sensiti ve user distribution at electricity price ratio 16. in AECO cite from PJM data miner [30]. W e consider three types of users: high-sensitive users, mid-sensitive users, and low-sensiti ve users, categorized based on the range of their electricity price elasticity . According to the definition of the uncertainty set (2), the elasticity is mainly influenced by the ratio between TOU price c T OU t and the reference price C RE F t . Therefore, here, we divide the range of TOU-to-reference ratio into five intervals, i.e., [0 . 00 , 0 . 25] , [0 . 25 , 0 . 50] , [0 . 50 , 1 . 00] , [1 . 00 , 4 . 00] , [4 . 00 , 16 . 00] , respecti vely . For each type of users, JOURNAL OF L A T E X CLASS FILES, VOL. XX, NO. X, FEB. 2025 7 we characterize their v ariation range of electricity price elastic- ity for each of these interv als. Fig. 3 sho ws three representativ e results. For e xample, Fig. 3(a) shows the electricity price elasticity of high-sensiti ve users under c T OU t /C RE F t = 0 . 3 . As we can see, the elasticity ranges from − 0 . 97 to − 0 . 11 , which is chosen as the parameters in the uncertainty set. Similarly , we can determine the parameters for all types of users under all TOU-to-reference ratios, which gives the uncertainty set. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 27 28 29 30 31 32 33 26 19 20 21 22 24 25 23 U p s t ream gri d G en erato r L o ad Fig. 4. The topology of the IEEE-33 bus system. B. Benchmark Results Using the uncertainty set generated in the last subsection and the proposed model and method, we can get the results for the IEEE-33 bus systems. The topology of the system is sho wn in Fig. 4. The proposed algorithm con verges in 3 iterations, which takes 115.44 seconds. The time needed is acceptable for the operation of virtual po wer plants. The change of LB and U B , defined in Algorithm 1, during the iteration processes is shown in Fig. 5; and they con verge to -$1645. Particularly , the negati ve value implies that the system is gaining profits by selling electricity . The optimal TOU prices obtained and the reference prices are sho wn in Fig. 6. The a verage predicted load L it of all buses and the av erage worst-case load profiles are given in Fig. 7. As we can see, the optimal T OU price curve follows a similar patten to the load demand curv e in Fig. 7, i.e., the price is relativ ely high during periods with large demand while the price is relati vely low when the load demand is small. This indicates that the TOU prices set by the proposed method can more effecti vely reflect the change of needs in the market than the reference price curv e (blue curve in Fig. 6), which is almost flat. Moreov er, the three load profiles under worst-cases in Fig. 7 are all smoother than the predicted load profile. This indicates that, by realizing and lev eraging the demand flexibility through TOU prices, we can achiev e peak shaving and v alley filling, and thus, reduces the electricity bills. C. Necessity of the pr oposed methods First, to show the necessity of the proposed RO algorithm, a traditional C&CG is applied to optimize the model for comparison. More specifically , instead of returning the v ∗ it to the master problem as in the proposed algorithm, the traditional C&CG algorithm returns the scenario ξ ∗ directly to master problem. The result of the traditional C&CG algorithm is shown in Fig. 8. As we can see, the traditional C&CG algorithm fails to conv erge, with the upper and lower bounds remain nearly unchanged. In addition, the value of LB is ev en larger than the v alue of UB, which is probably caused Fig. 5. Change of UB and LB during the iterations of Algorithm 1. Fig. 6. Optimal TOU and reference prices of a day . by the inaccurate estimation of ξ by ignoring its decision- dependency in the scenarios. This shows the necessity of a new R O algorithm to address the DDU. Next, to show the accuracy of the proposed transformation method in Section III-C3, we compare the result using the proposed transformation method and the one using nonlinear models. The latter is sho wn in Fig. 9. The nonlinear model takes 249.49 seconds to approximately conv erge, with the lower bound equaling -$1662.6 and the upper bound equaling -$1631.7, which are slightly lar ger than the -$1645 obtained by the proposed algorithm. This re veals that the original nonlinear model could result in a suboptimal result. D. Sensitivity Analysis Further , we analyze the impact of different factors to better show the effecti veness of the proposed method. 1) F ixed TOU price: W e test the influence of TOU prices. Instead of allo wing the VPP operator to optimize the TOU prices, we use fixed TOU prices (constant c T OU t , ∀ t ). W e let c T OU t be 1, 2, and 5 times of C RE F t , respecti vely . The results are summarized in T ABLE I. W e can find that, with the optimal TOU prices deri ved by the proposed method, the VPP operator can generate the highest profit (the lo west objectiv e value). This demonstrates the advantages of treating electricity prices as decision variables. JOURNAL OF L A T E X CLASS FILES, VOL. XX, NO. X, FEB. 2025 8 Fig. 7. Predicted and worst-case load profiles. Fig. 8. Result of traditional C&CG Algorithm. T ABLE I R E SU LT S U N D E R D I FFE R E N T T O U P R I C E S c T O U t Optimal C RE F t 2 C RE F t 5 C RE F t Obj ($) -1645 73.65 -260.8 -1178.3 2) Uncertainty budget: The level of uncertainty is also another important influencing factor , which is controlled by the uncertainty budgets Γ T and Γ S . In the benchmark case, we let Γ T = 24 and Γ S = 33 . Here, we test the impact of uncertainty lev el by changing these two parameters, with the results summarized in T ABLE II. As sho wn, the VPP operator’ s profit enhances with both reduced Γ T and Γ S . Howe ver , smaller uncertainty budgets also result in more number of iterations and longer computation time. Still, the time needed is acceptable for the operation of VPP . T ABLE II R E SU LT S U N D E R D I FFE R E N T U N C ER T A I N TY L EV E L S (Γ T , Γ S ) (24 , 33) (20 , 33) (24 , 25) (20 , 25) Obj ($) -1645 -1648 -1667 -1668 T ime (s) 115.44 287.95 276.36 283.02 No. of Iteration 3 4 4 5 V . C O N C L U S I O N This paper considers the optimal operation of virtual power plants under DDU in the price elasticity of demand. W e de- velop a tw o-stage robust VPP pricing model to determine T OU Fig. 9. Result of Nonlinear Model Con vergence Figure tariffs. The uncertain demand is influenced by the T OU pricing decisions. T o solve this RO problem with DDU ef ficiently , we customize an improv ed C&CG algorithm. This algorithm maps the previously selected worst-case scenarios to the vertices of the ne w uncertainty set, and can address a decision-dependent mixed-inte ger uncertainty set. The case studies re veal that compared to the traditional C&CG algorithm, the proposed algorithm offers conv ergence guarantee under DDU. In the future, we may further explore methods to more precisely characterize the decision-dependency of uncertainty , as well as more effecti ve RO algorithms for addressing DDU. R E F E R E N C E S [1] IEA, “Global energy revie w 2025, ” IEA, Paris, France, T ech. Rep., 2025. [Online]. A vailable: https://www .iea.org/reports/ global- energy- re view- 2025 [2] G. Ruan, D. Qiu, S. Si varanjani, A. S. A wad, and G. Strbac, “Data- driv en energy management of virtual power plants: A review , ” Advances in Applied Energy , vol. 14, p. 100170, 2024. [3] Z. Ullah, G. Mokryani, F . Campean, and Y . F . Hu, “Comprehensive revie w of VPPs planning, operation and scheduling considering the uncertainties related to renew able energy sources, ” IET Ener gy Systems Inte gration , vol. 1, no. 3, pp. 147–157, 2019. [4] P . Shinde, I. Kouveliotis-L ysikatos, and M. Amelin, “Multistage stochas- tic programming for VPP trading in continuous intraday electricity markets, ” IEEE T ransactions on Sustainable Ener gy , vol. 13, no. 2, pp. 1037–1048, 2022. [5] X. Wu, H. Xiong, S. Li, S. Gan, C. Hou, and Z. Ding, “Improved light robust optimization strategy for virtual power plant operations with fluctuating demand, ” IEEE Access , v ol. 11, pp. 53 195–53 206, 2023. [6] S. Y u, F . F ang, and J. Liu, “Flexible operation of a CHP-VPP considering the coordination of supply and demand based on a strengthened dis- tributionally robust optimization, ” IET Contr ol Theory & Applications , vol. 17, no. 16, pp. 2146–2161, 2023. [7] J. Cao, B. Y ang, C. Y . Chung, Y . Gong, and X. Guan, “Distribu- tionally robust management of hybrid ener gy station under exogenous- endogenous uncertainties and bounded rationality , ” IEEE T ransactions on Sustainable Energy , vol. 15, no. 2, pp. 884–902, 2023. [8] D. Xiang, K. He, G. Zhang, S. Zhang, Y . Chen, L. Zhao, H. Liu, X. Su, and H. Zheng, “Optimizing VPP operations: A novel robust demand- side management approach with ener gy storage system, ” in 2023 IEEE International Conference on Ener gy Internet (ICEI) . IEEE, 2023, pp. 45–50. [9] H. Nemati, I. Egido, P . S ´ anchez-Mart ´ ın, and ´ A. Ortega, “ Assessing value of renewable-based VPP versus electrical storage: Multi-market participation under dif ferent scheduling regimes and uncertainties, ” 2025. [Online]. A vailable: https://arxiv .org/abs/2507.22496 [10] D. Bertsimas, E. Litvinov , X. A. Sun, J. Zhao, and T . Zheng, “ Adaptive robust optimization for the security constrained unit commitment prob- lem, ” IEEE T ransactions P ower Systems , v ol. 28, no. 1, pp. 52–63, Feb. 2013. JOURNAL OF L A T E X CLASS FILES, VOL. XX, NO. X, FEB. 2025 9 [11] B. Zeng and L. Zhao, “Solving two-stage rob ust optimization problems using a column-and-constraint generation method, ” Operations Resear ch Letters , vol. 41, no. 5, pp. 457–461, 2013. [12] T . T an, M. Y ang, R. Xie, Y . Cao, and Y . Chen, “ Adjustable robust optimization with decision-dependent uncertainty for po wer system problems: A review , ” Energy Con version and Economics , 2025. [13] Y . Zhang, F . Liu, Z. W ang, Y . Su, W . W ang, and S. Feng, “Robust scheduling of virtual power plant under exogenous and endogenous uncertainties, ” IEEE T ransactions on P ower Systems , vol. 37, no. 2, pp. 1311–1325, 2021. [14] N. Qi, L. Cheng, H. Li, Y . Zhao, and H. T ian, “Portfolio optimization of generic energy storage-based virtual power plant under decision- dependent uncertainties, ” Journal of ener gy storag e , vol. 63, p. 107000, 2023. [15] L. Hellemo, P . I. Barton, and A. T omasgard, “Decision-dependent probabilities in stochastic programs with recourse, ” Computational Man- agement Science , vol. 15, no. 3, pp. 369–395, 2018. [16] K. Zhou, N. Peng, H. Y in, and R. Hu, “Urban virtual power plant operation optimization with incenti ve-based demand response, ” Energy , vol. 282, p. 128700, 2023. [17] S. Mei, Q. T an, Y . Liu, A. T rivedi, and D. Srinivasan, “Optimal bidding strategy for virtual power plant participating in combined electricity and ancillary services market considering dynamic demand response price and integrated consumption satisfaction, ” Ener gy , vol. 284, p. 128592, 2023. [18] X. K ong, W . Lu, J. W u, C. W ang, X. Zhao, W . Hu, and Y . Shen, “Real- time pricing method for VPP demand response based on PER-DDPG algorithm, ” Ener gy , v ol. 271, p. 127036, 2023. [19] G. Lu, B. Y uan, S. Zhou, L. W ei, and Z. W u, “ Assessing the effecti veness of time-of-use pricing design: Provincial evidence from China, ” Ener gy Strate gy Reviews , vol. 60, p. 101780, 2025. [20] D. S. Kirschen and G. Strbac, Fundamentals of power system economics . W iley , 2004. [21] Y . Chen and W . W ei, “Robust generation dispatch with strategic re- new able power curtailment and decision-dependent uncertainty , ” IEEE T ransactions on P ower Systems , vol. 38, no. 5, pp. 4640–4654, 2022. [22] R. Xie, P . Pinson, Y . Xu, and Y . Chen, “Robust generation dispatch with purchase of renewable po wer and load predictions, ” IEEE T ransactions on Sustainable Energy , vol. 15, no. 3, pp. 1486–1501, 2024. [23] T . T an, R. Xie, X. Xu, and Y . Chen, “ A robust optimization method for power systems with decision-dependent uncertainty , ” Ener gy Con version and Economics , vol. 5, no. 3, pp. 133–145, 2024. [24] M. Y ang, R. Xie, Y . Zhang, and Y . Chen, “Robust microgrid dispatch with real-time energy sharing and endogenous uncertainty , ” IEEE Tr ans- actions on Smart Grid , 2025. [25] Y . Zhang, F . Liu, Y . Su, Y . Chen, Z. W ang, and J. P . Catal ˜ ao, “T wo-stage robust optimization under decision dependent uncertainty , ” IEEE/CAA Journal of Automatica Sinica , vol. 9, no. 7, pp. 1295–1306, 2022. [26] S. A vraamidou and E. N. Pistikopoulos, “ Adjustable robust optimization through multi-parametric programming, ” Optimization Letters , v ol. 14, no. 4, pp. 873–887, 2020. [27] P . V ayanos, A. Georghiou, and H. Y u, “Robust optimization with decision-dependent information discovery , ” Management Science , 2025. [28] Y . Zhang, Y . Su, and F . Liu, “On decision-dependent uncertainties in power systems with high-share renewables, ” Engineering , 2025. [29] U. P . Networks, “Smartmeter energy consumption data in london households, ” 2014, accessed: 2025-9-20. [Online]. A vailable: https: //innov ation.ukpowernetworks.co.uk/projects/lo w- carbon- london [30] PJM, “PJM data miner 2, ” 2024, settlement verified hourly LMPs of PJM market. [Online]. A vailable: https://dataminer2.pjm.com/list A P P E N D I X A P RO O F O F P RO P O S I T I O N 1 Suppose O ∗ is the optimal value of problem (8) and the solution ( x ∗ , z ∗ , ξ ∗ ) is optimal. W e prov e Proposition 1 by showing three claims: 1) Claim 1: LB N ≤ O ∗ ≤ U B N for any N . Proof of Claim 1: According to the master problem (16), LB N = min x,z C ⊤ x + max ξ,v min y E ( x ) ⊤ y , (A.1a) s.t. Ax ≥ B , (7b)–(7d) , P y ≥ Q ( x ) ξ + R ( x ) , (A.1b) ξ it = X k ∈K z tk  (1 − v it ) ξ − itk + v it ξ + itk  , (A.1c) v ∈ { v 1 ∗ , v 2 ∗ , . . . , v N − 1 ∗ } . (A.1d) Because O ∗ = min x,z C ⊤ x + max ξ,v min y E ( x ) ⊤ y , (A.2a) s.t. (A.1b) , (A.1c) , v ∈ V , (A.2b) and { v 1 ∗ , v 2 ∗ , . . . , v N − 1 ∗ } ⊂ V , problem (16) is a relaxation of problem (8), which implies LB N ≤ O ∗ . According to the subproblem (11), U B N = min x,z C ⊤ x + max ξ min y E ( x ) ⊤ y , (A.3a) s.t. Ax ≥ B , (7) , P y ≥ Q ( x ) ξ + R ( x ) , (A.3b) x = x N ∗ , (A.3c) so O ∗ ≤ U B N . 2) Claim 2: If N 1 < N 2 and Algorithm 1 does not con verge after N 2 iterations, then v N 1 ∗  = v N 2 ∗ . Proof of Claim 2: W e prove it by contradiction. Suppose v N 1 ∗ = v N 2 ∗ , then v N 2 ∗ ∈ { v 1 ∗ , v 2 ∗ , . . . , v N 2 − 1 ∗ } . Using the optimal solution x N 2 ∗ of the master problem (16) in the N 2 -th iteration and the equation in (A.1), we have LB N 2 = min x,z C ⊤ x + max ξ,v min y E ( x ) ⊤ y , (A.4a) s.t. (A.1b) , (A.1c) , (A.4b) v ∈ { v 1 ∗ , v 2 ∗ , . . . , v N 2 − 1 ∗ } , x = x N 2 ∗ . (A.4c) Substituting { v 1 ∗ , v 2 ∗ , . . . , v N 2 − 1 ∗ } by { v N 2 ∗ } , we hav e LB N 2 ≥ min x,z C ⊤ x + max ξ,v min y E ( x ) ⊤ y , (A.5a) s.t. (A.1b) , (A.1c) , (A.5b) v = v N 2 ∗ , x = x N 2 ∗ , (A.5c) where the right-hand side equals U B N 2 according to (A.3) and the optimality of v N 2 ∗ in the subproblem (11). Thus, LB N 2 ≥ U B N 2 . Combining it with Claim 1, we have U B N 2 = LB N 2 and Algorithm 1 conv erges after N 2 iterations, which is a contradiction. 3) Claim 3: Algorithm 1 con verges within N iterations. Proof of Claim 3: For any iteration N , the worst-case scenario ξ N ∗ always resides at a verte x of the uncertainty set (7) [21]. Therefore, v N ∗ can be achieved at a vertex of the set V = [0 , 1] I T . By Claim 2, a verte x of V cannot be found twice unless Algorithm 1 con ver ges. Therefore, the number of iterations cannot exceed the number of vertices of V plus 1, which is N = 2 I T + 1 . According to Claim 3, we can assume Algorithm 1 con- ver ges after the N 0 -th iteration, where N 0 ≤ N . By the con vergence criterion, we ha ve U B N 0 = LB N 0 . Thus, Claim 1 implies LB N 0 = O ∗ = U B N 0 , which proves Algorithm 1 finds the optimal solution. This completes the proof.