Mixed Non-Linear Economic Dispatch Using G-SDDP with Unified Benders Cuts Jesus Velásquez-Bermúdez, Eng. D. Chief Scientist, DecisionWare - DO Analytics LLC
[email protected] IEEE PES T&D Latin America LIMA 2018 Abstract - This paper presents the implementation of the Generalized Dual Dynamic Programing (GDDP, Velásquez 2002), and its use in electric sector applications. The GDDP is based on the chained application Benders Partition Theory (Benders, 1962) to solve multi-period dynamic problems using the Dynamic Programming approach. In the implementation is used the concept of Unified Benders Cuts) (Velásquez, 2018a). Applications to Electric Systems are presented using Generalized Stochastic Dual Dynamic Programming (GSDDP). Key words: Benders Decomposition, Nested Benders, Dynamic Programming, Dual Dynamic Programming, Discrete Control Theory, Unified Benders Cuts, Mixed Non-Linear Economic Dispatch.
Format Problem LP MIP NLP MINLP
Benders Master LP MIP NLP MINLP
Benders Subproblem LP LP LP LP
DDP
GDDP
YES NO NO NO
YES YES YES YES
In the papers “Benders Decomposition using Unified Cuts”, “GDDP: Generalized Dual Dynamic Programming. Theory & Electric Sector Applications” and “G-SDDP: Generalized Stochastic Dual Dynamic Programming. Theory & Electric Sector Applications” the reader can complement this paper. B.
Benders' Partition Theory
BT considers the problem P: composed by two types of variables: y, the coordination variables, and x, the coordinates.
I. INTRODUCTION A.
Background
This research presents theoretical considerations about the solution of dynamic optimization problems integrating Benders Theory (BT) and Dynamic Programming (DP). The Generalized Dual Dynamic Programming Theory (GDDP) was published by Velásquez in 2002; it can be considered as an extension of Dual Dynamic Programming (DDP) (Pereira and Pinto, 1985, 1991, Velásquez 1999) which is a variation of Nested Benders (NB) methodologies (Birge 1985, Murphy 2013). G-SDDP is the stochastic version of GDDP. The Benders Decomposition review presented by Rahmaniani et al. (2016) doesn’t included reference to a work like GDDP.
P: = { min z = cTx + f(y) | F0(y) = b0 Ax + F(y) = b xR+ , yS }
(1)
BT restricts the model on x to be a linear problem, while it doesn’t impose conditions on y that may be continuous, or discrete. Additionally, the functions f(y) and F(y) may be nonlinear convex functions. The P: problem is partitioned in two coordinated problems: CY: over y and SP(y): over x which is defined as SP(y): = { min Q(y) = cTx | Ax = b - F(y), xR+ }
(2)
The coordinator CY: on y can be formulated as DDP solves dynamic multiple-stage problems representing the future cost function by the hyperplanes (generated using BT) that define it. DDP only solves linear models, because it only considers state variables. GDDP makes a distinction between state and control variables. This distinction permits a more detailed algorithm in which the sub-problems are smaller than in DDP; then GDDP solves a large-scale optimization problem as a coordinated sum of very simple problems; it entails a significant speed off the solution time. DDP limits the relations of the state variables to two consecutive periods; Velásquez (1995) applied BT to a more general dynamic problem which allows direct relationship between the state variables of two non-contiguous periods; this case may be included in GDDP. Finally, GDDP is capable to solve the following problems:
CY: = { min z = f(y) + Q(y) | F0(y) = b0 , yS Q(y) (k)T[b - F(y)] k IT 0 (k)T[b - F(y)] kIN }
(3)
where represents the dual variables vector of the restrictions Ax = b - F(y), IT the set of iterations, an extreme ray on the feasibility region and IN the set of iterations on which no feasibility was obtained. Benders proposed the solution of P: by a hierarchical algorithm that works on two levels: i) the coordination level solves the problem CY: and generates a sequence of yk values; ii) on the
second level, yk is used as a parameter of the sub-problem SP(y): to generate a sequence of feasible extreme points, k, or a sequence of extreme rays, k, of the dual feasible zone of SP(y):, this vectors are used to include cutting planes in CY:.
DDP implies the multilevel decomposition (in the time domain) of multiple problems using BT, generating new problems that again are decomposed using BT (Nested Benders). The reduce form of previous problem is:
Then CY: includes two types of cuts: i) optimality cutting planes (OCP) restrict the feasible zone of y to obtain the optimal y; it has the following structure
DDP: = { min z = t= ctTxt | At xt + Et-1 xt-1 = bt t= xtR+ t=}
Q(y) (k)T[b - F(y)]
k IT
(4)
ii) feasibility cutting planes (FCP) restrict the feasible zone of y to maintain feasible x in SP(y), it has the following structure 0 (k)T[b - F(y)]
kIN
As a starting point, DDP divided the variables in two groups: i) coordination variables associated with stages/periods between 1 and T-1 {x1, x2, ... , xT-1}, and ii) coordinated variables associated with stage T {xT}. The coordinator model for the period {1,T-1} is
(5)
CDT-1: = { min z = t=1,T-1 ctTxt + T(xT-1) | At xt = bt - Et-1 xt-1 t=1,T-1 , xtR+ t=1,T-1) T(xT-1) + kTTET-1 xT-1 kTTbT kIT(T-1) } (11)
For simplicity, in the following sections FCPs are ignored, but they can be included following the same way used for OBCs. C.
DP: Dynamic Programming
Dynamic Programming (DP, Bellman 1952) is a numerical methodology that solves optimization problems related to dynamic techno-economic systems. For this study, we are interested in a DP problem in which all the mathematical elements must be separable in the time domain, with mixednon-linear-convex relations for the state variables (xt) and linear relations for control variables (ut). DP: = { min z = t= [ ct(xt) + dt ut ] | Ft(xt-1, xt ) = ft t= At(xt-1, xt) + Btut = bt t= Gt ut = gt t= utR+ , xtR+ , t= }
(7)
CDt: = { min z = = cTx + t+1(xt) | Ax = b - E-1x-1 = , xR+ = t+1(xt) + jt+1TEtxt t+1j jIT(t) }
(12)
(8)
If the vector yt is divided in state and control variables we can see the mathematical problem as:
(13)
and the sub-problem for each intermediate stage t is SDt(xjt-1): = { min t(xjt-1) = ctTxt + t+1(xt) | Atxt = bt - Et-1xjt-1 , xtR+ t+1(xt) + kt+1TEtxt t+1k kIJ(t,j) } (14) where IJ(t,j) represents the number of cuts included in the coordinator CDt: when we solve the sub-problem SDt(xjt-1). tj is a constant value calculated as tj = {
yt represents, jointly, state variables and control variables
DDP: = { min z = t= (ctTxt + dtTut ) | Ft xt = ft t=T At xt + Et-1 xt-1 + Bt ut = bt t= Gt ut = gt t= utR+ , xtR+ , t= }
|
(6)
A linear model (DDPO:) is considered in DDP and NB
yt = { xt , ut }
For the stage T the linear sub-problem SDT(xkT-1) is defined as
Using a recursive approach, the coordinator for each intermediate stage t (less than T) is
DDP: Dual Dynamic Programming
DDPO: = { Min t=1,T ctT yt | At yt = bt - Et-1 yt-1 , t=1,T ; ytR+ }
where t(xt-1) represents the future costs function considering that the system is on state xt-1 at the beginning of period t, IT(t) the number of cuts included in the coordinator CDt, and kt the vector of dual variables of the restrictions At xt = bt - Et-1 xkt-1 obtained as a solution of the sub-problem SDt(xkt-1). We call these cuts as Future Cost Benders Cuts (FCBC).
SDT(xkT-1): = { min T(xkT-1) = cTTxT ATxT = bT - ET-1xkT-1 , xTR+ }
where t represents the periods, xt the state variables vector and ut the control variables vector. Ft() and At() techno-economic functions, Et , Bt , and Gt are techno-economic matrices, ft , bt and gt are resources vectors, ct() cost function for the state variables and dt are cost vector for the control variables. D.
(10)
jTTbT jtTbt
+ k=1,IJ(t,j) k,jt t+1k
t=T t = 1,T-1
(15)
where k,jt represents the dual variable of the k-th Benders cut in the sub-problem SDt(xjt-1). The index j represents the iteration of the coordinator of SDt(xjt-1) and the index k the cuts that have been included in the process of solution of SDt(xjt-1). II. GDDP: Generalized Dual Dynamic Programming A.
(9)
where algebraic symbols are the same used in DP, and matrices Ft and At, corresponds to Ft() and At() functions.
GDDP Basic Theory
Velásquez (2002) considered the solution of the DDP problem dividing the variables in state variables and control variables; then mathematical problem may be relaxed to include: i) discrete variables in the state variables vector; ii) non-linear
convex relations and non-linear separable convex cost functions between the state variables. The GDDP: problem is: GDDP: = { min z = c(x1 , x2 , … , xT) + t= dtTut | F(x1 , x2 , … , xT) = f A(x1 , x2 , … , xT) + Btut = bt t= Gt ut = gt t= utR+ t= xtSt t= } (16) where ct(x) is the cost function of the state variables, Ft(x) and At(x) vectorial functions, and St the domain of xt.. The solution of GDDP: using BT considers a two-stage process: i) first, we define the state variable xt as the coordination variables to decouple the problem at temporary level solving static subproblems SUt(xt-1,xt): for each t, optimizing the control variables ut and ii) in the second stage, the GDDP solves the coordinator problem, CX:, optimizing the state variables xt. CX: is stated as CX: = { min z = c(x1 , x2 , … , xT) + t= t(xt-1,xt) | F(x1 , x2 , … , xT) = f { min t(xt-1,xt) = dtTut | Btut = bt - A(x1 , x2 , … , xT) Gtut = gt , utR+ } t= xtS t= } (17) where t(xt-1,xt) represents the optimum operation costs in the period t as consequence of the border condition (x1 , x2 , … , xT) and it corresponds to the objective function of the linear subproblems SUt(xt-1,xt): defined as SUt(xt-1,xt): = { min t(xt-1,xt) = dtTut | Btut = bt - A(x1 , x2 , … , xT) , Gtut = gt , utR+ }
(18)
The dual problem of SUt(xt-1,xt): is
B.
Unified Benders Cuts
UBCs are related to the special case when the dual feasibility zone of the problems SUt(xt-1,xt) is static, which implies that the matrices Bt and Gt and the dt vector are time independent; then DSUt(xt-1,xt): = { max t(xt-1,xt) = tT[bt - Et-1xt-1 - Atxt] + tT gt | TB + TG dT }
(21)
In this case, a feasible vector of dual variables { ,} for any DSUt(xt-1,xt): is feasible for all DSUt(xt-1,xt): independent of the value of t. then, a solution of SUt(xt-1,xt) for any value of t can generate supply Benders cuts (SBC) for all periods. Then the coordinator CX: can be expressed as CX: = { min z = c(x1 , x2 , … , xT) + t= t(xt-1,xt) | F(x1 , x2 , … , xT) = f t(xt-1,xt) + (k)TEt-1 xt-1 + (k)TA(x1 , x2 , … , xT) (k)Tbt + (k) Tgt t= kIU } (22) where dual variables { , } are time independent. Then, for each iteration coordinator-sub-problems to solve only one problem SUt(xt-1,xt) generates cuts for all periods. IU represents all iterations at sub-problems level. From an economic point of view, t(xt-1,xt) represents the optimum operation costs as a consequence of the boundary conditions (xt-1 and xt); frequently, it is time independent, mainly in short/medium term planning, due to the fact that in these cases the technology, the topology and the price index may don’t change; then, is common for matrices Gt and Bt because they are related with the technology and topology of the industrial system. In this case, we call the Benders cuts as Unified Benders Cuts (UBC, Velásquez 1995, 2002, 2018a).
where t represents the dual variables of the restrictions Bt ut = bt - A(x1 , x2 , … , xT) and t is the dual variables of Gt ut = gt.
UBCs are very important when: i) the number periods is very large, such as: i) discrete control theory (DCT), that require many short periods to represent the continuous movement and derivatives constraints of the state variables; and ii) stochastic optimization that requires many synthetic scenarios to represent continues probabilistic functions.
Considering the decoupled cuts (Velásquez, 2018b) generated by each sub-problem SUt(xt-1,xt): the coordinator CX: is
III. G-SDDP: Generalized Stochastic Dual Dynamic Programming
DSUt(xt-1,xt): = { max t(xt-1,xt) = tT[bt - A(x1 , x2 , … , xT) + tTgt | tTBt + tTGt dtT } (19)
CX: = { min z = c(x1 , x2 , … , xT) + t= t(xt-1,xt) | F(x1 , x2 , … , xT) = f t(xt-1,xt) + (tk)TEt-1 xt-1 + (tk)TA(x1 , x2 , … , xT) (tk)Tbt + (tk) Tgt t= kIU(t) } (20) where IU(t) represents the number of cuts generated for each sub-problem SUt(xt-1,xt). The cuts that integrated the coordinator CX: are called Supply Benders Cuts (SBC), because they represent the cost supply function of the industrial system parametrized on the values of the state variables, xt. If GDDP: is linear, CX: is linear and has a dynamic structure like DDP:; then CX: may be solved using the DDP theory (Velásquez, 2002). For mixed and/or non-linear problems, CX: must be solved in an integrated form. Velásquez (2018b) shows that integrated coordinator is better than DDP coordinator.
A.
G-SDDP Basic Theory
The G-SDDP problem corresponds to a multi-stage stochastic optimization problem that is linear for the control variables (ut,h) and may be mixed non-linear for the state variables (xt,h). The random scenarios index h can be: i) predefined (read) or ii) be synthetically generated within the optimization process. Stochastic modeling is supported on split variables approach (Jornsten, 1985) which is based on creating multiple parallel problems, one for each scenario, all problems are coordinated using non-anticipative constraints that represent the decision tree; it implies to add the h-index to all variables and constraints of the “core” deterministic problem. The basic structure of the problem, without non-anticipative constraints is
G-SDDP: = { Min z =h h [ ch (x1,h , x2,h , … , xT,h) + t= dt,h Tut,h ] | Ft(x1,h , x2,h , … , xT,h) = fh hN At(x1,h , x2,h , … , xT,h) + Bt,h ut,h = bt,h t= hN Gt,h ut,h = gt,h t= hN ut,hR+ t= hN xt,hSh t= hN } (23) where h represents the probability of the scenario h and the objective function is to minimize the expected cost. The n-index is associated with nodes of the decision tree, then the non-anticipative constraints are uAt,h = ut,h t nN(t), h(n) (24) xAt,h = ut,h t nN(t), h(n) (25) where N(t) represents the set of nodes of decision tree related with period t and B(n) the set of scenarios related with node n. G-SDDP: solution of is based on partition and decomposition using BT, structuring the coordinated solution of multiple smaller subproblems, that considers a two-stage process: i) To define xt,h as the coordination variables to decouple the problem at time-scenario level; solving, for each static random subproblems SUt,h(xt-1,h,xt,h):; and ii) To solve the coordinator problem, CX: optimizing/defining the state variables xt,h. B.
G-SDDP Subproblem
G-SDDP solves for every stage of each scenario, a staticrandom sub-problem SUt,h(xt-1,h,xt,h): SUt,h(xt-1,h,xt,h): { min t,h(xt,h-1,xt,h) = dt,hT ut,h | Bt,h ut,h = bt,h - At,h (x1,h , x2,h , … , xT,h) Gt,hut,h = gt,h , ut,hR+ } t= h
()
The dual problem of SUt,h(xt-1,h,xt,h): is DSUt,h(xt,h-1,xt,h): = { max t,h(xt,h-1,xt,h) = t,hT[bt,h - A(x1 , x2 , … , xT,H) + t,hTgt,h | t,hTBt,h + t,hTGt,h dt,hT } (27) where t,h represents the dual variables vector of the restrictions Bt,h ut,h = bt,h - A(x1,h , x2,h , … , xT,h) and t,h the dual variables vector of Gt,h ut,h = gt,h. C.
G-SDDP Coordinator
The coordinator CX:, that must include state variable constraints, Benders cuts and non-anticipative constraints, is CX: = { min z = h h [ ch(x1,h , x2,h , … , xT,h) + t= t,h(xt-1,h,xt,h) ] | Ft,h(x1,h , x2,h , … , xT,h) = fh { min t,h(xt,h-1,xt,h) = dt,hT ut,h | Bt,h ut,h = bt,h - At,h (x1,h , x2,h , … , xT,h) Gt,hut,h = gt,h , ut,hR+ } t= h uAt,h = ut,h t nN(t) h(n) xAt,h = xt,h t nN (t) h(n) xt,hSt,h t= h } (28)
Considering the decoupled cuts generated by each sub-problem SUt,h(xt-1,h,xt,h): the coordinator CX: is CX: = { min z = h h [ ch(x1,h , x2,h , … , xT,h) + t= t,h(xt-1,h,xt,h) ] | Ft,h(x1,h , x2,h , … , xT,h) = fh t(xt-1,xt) + (t,hk)TEt-1,h xt-1,h + (t,hk)TA(x1,h, x2,h , … , xT,h) (t,hk)Tbt,h + (t,hk) Tgt,h t= kIU(t,h) uAt,h = ut,h t nN(t) h(n) xAt,h = xt,h t nN (t) h(n) xt,hSt,h t= h } (29) where IU(t,h) represents the number of cuts (SBC) generated for each sub-problem SUt,h(xt-1,h,xt,h). IV. G-SDDP Implementation In this case G-SDDP solves the coordinator CX: in an integrated problem; for the integrated solution it is considered a replacement of the pair of indexes by a single equivalent index th, ordered by scenarios and then by period, the new problem has equivalent structure with the coordinator problem solved by the GDDP, except for the continuity in the state variables when changing scenario h to h+1 (Velásquez, 2018b). There are three alternatives to implement Benders’ Cuts: ▪ SBC (Standard Bender’s Cut): corresponds to standard BT that solves a single subproblem integrating all periods. ▪ DBC (Decoupled Bender’s Cuts): corresponds to a variation of BT that solves T problems and generates one decoupled cut for each period; ▪ UBC (Unified Bender’s Cuts): when the mathematical conditions are met, it corresponds to a variation of BT that solves a subset of N subproblems of all subproblems and generates N decoupled cuts for each period. Implementation of DBC and UBC cuts implies two levels of loops: i) the outer loop moves iterations counter of BT; and on ii) the inner loop the subproblems associated with periodsscenarios. The difference between the two alternatives is: ▪ UBC solves N subproblems (1 ≤ N ≤ T×H) and DBC always solves T×H subproblems. ▪ UBC add a cut for each period it is generated by each subproblem, independently of the time-scenario associated to the subproblem (N × T×H cuts) and DBC always add one cut for each time-scenario. When it is possible to use UBC has a better performance than DBC (Velásquez 2018a, 2018b). V. Economic Dispatch Modeling GDDP formulation represents many industrial systems in which state variables (xt) may be associated with: i) the stocks (inventories), ii) dynamic capacities and ii) the state of the production units (on-off); and control variables (ut) with the production and distribution activities through the supply-chain (production factories) which are separable between periods. The implementation of GDDP/G-SDDP methodologies were applied to an electric system selected from the IEEE test cases (Diniz, 2010) composed of 12 hydro plants, 8 reservoirs, 6 thermal plants, 1 “deficit” plant and 1 node of electricity demand; electricity transmission network wasn’t considered.
REFERENCE ELECTRIC SYSTEM
A.
The constraints of LP-ED problem can be stated as (between parenthesis the constraint name):
Actual Plants
~
Hydro Plant
Units
Capacity (MW)
i01
Ui01 - Ui12
85
Ui01-Ui02 Ui03-Ui06 Ui07-Ui08 Ui01 - Ui03 Ui01 Ui01-Ui02 Ui01-Ui04 Ui01-Ui02 Ui01-Ui04 Ui01-Ui03 Ui01-Ui03 Ui04-Ui06 Ui01-Ui02 Ui03-Ui05 Ui01-Ui02 Ui03-Ui04 Ui05-Ui06
38 51 54 180 40 130 200 15 300 125 380 380 20 52 150 180 280
i02
HYDRO-ELECTRIC HYDRO PLANTS - CONFIGURATION BUS Actual Capacity Reservoirs (hm3) i01 17000 i02 2600 i05 200 i06 10000 i08 12500 THERMO-ELECTRIC i09 1000 i10 12400 i12 5000
Linear Economic Dispatch - Constraints
CONSUMER
RESERVOIR
i03 i04 i05 i06 i07 i08 i09 i10 i11
12 8 6 1 24
Hydro Plants Reservoirs Thermal Plants Deficit Plant Hour Planning Horizon
i12
Total Capacity (MW) 1020
Total Capacity (MW)
388 540 40 260 800 30 1200 375 2280 196 1220
8349
Three types of applications were studied: i) Linear Economic Dispatch (LP-ED), ii) Mixed Economic Dispatch (MIP-ED), iii) Mixed–Non-Linear Economic Dispatch (MINLP-ED). The full formulation is in Velásquez (2018c). The following tables present the reduced math formulation for the determinist case. To obtain the stochastic formulation is necessary: 1. To include the scenario index h in: i) all variables and constraints of the determinist case; and ii) in the parameters that are associated with random scenarios; in this case water inflow and demand. 2. To define variables and constraints required to guarantee the non-anticipative character of the model. Index i,p d j k t
Hydro plant Deficit plant Thermal unit Hydro unit Period
Sets I J T ITER UI(i) UP(i)
Description Hydro plants Thermal generation units Time stages Active iterations Hydro generation units installed in plant i Hydro plants immediately upstream of plant i
Operational Limits t i t i t i kUI(i) t i kUI(i) jJ t Turbinated outflow (QQHt,i) Qt,i = kUI(i) QHt,k t i Water balance (WATt,i) Vt,i - Vt-1,i = At,i + pUP(i) (Qt,p + SPt,p) t i – (Qt,i + SPt,i) Hydraulic generation (GQTt,i) GHt,k = k × QHt,k t i kUI(i) Total hydraulic generation (GUTt,i) GHt =iI kUI(i) UHt,k t Electricity demand (ELDt,i) GHt +jJ GTt,j = Lt tT Vmini ≤ Vt,i ≤ Vmaxi 0 ≤ Qt,i ≤ Qmaxi 0 ≤ QHt,k ≤ QHmaxk 0 ≤ GUt,k ≤ GUmaxt,k 0 ≤ GTt,j ≤ GTmaxj
(35) (36) (37) (38) (39) (40) (41) (42) (43) (44)
For all models the objective is the minimum dispatch cost, subject to previous constraints. The LP-ED objective function (45) Min t jJ CTj × GTt,j LP-ED formulation is always feasible for any value of electricity demand due the inclusion of an additional thermal plant with infinite generation capacity (plant d). For the case of level of the reservoirs it is feasible ever than the total inflows plus the initial water storage is enough to satisfy the minimum level of the reservoirs. For simplicity, in the experiments Vmini was defined as 0.0. This is important because the application of BT does not require of Feasible Benders Cuts (FBCs).
Description
B.
Parameter Lt k GHmaxk QHmaxk Vmini Vmaxi Qmaxi At,i GTmaxj PUmini PUmaxi CTj CARR CU1i CU2i
Description Electricity demand Generation efficiency Upper hydro generation bound Upper hydro unit outflow Lower reservoir volume Upper reservoir volume Upper turbined outflow limit Natural water inflow Upper thermal generation limit Pumping station lower limit Pumping station upper limit Thermal plant variable generation cost Thermal plant start cost Pumping cost (linear coefficient) Pumping cost (quadratic coefficient)
Units MWh MWh/(m3/s) MWh m3/s hm3 hm3 m3/s m3/s MW/h hm3/period hm3/period $/MWh $ $/hm3 $/(hm3)2
Variable UHt,k QHt,k GHt Vt,i Qt,i SPt,i GTt,j SRt,j STt,j PUt,i,p PWt,i BPt,i
Description Hydro unit outflow Hydro unit power generation Total hydraulic generation Reservoir operating volume Turbined outflow Spillage outflow Thermal generation Start thermal plant (0 stop, 1 start) State thermal plant (0 off, 1 on) Water pumped from i to p Power used in pumping station i Binary pumping (1 pumping)
Units m /s MW MW hm3 m3/s m3/s MW 0-1 0-1 hm3 hm3 0-1
Mixed-Linear Economic Dispatch - Constraints
The MIP-ED model includes all the restrictions of the LP-ED model plus the constraints required to represent the start/stop process in thermal plants, which requires binary variables that represent this process. Thermal availability (STNt,j y STXt,j) STt,j × GMINt,j ≤ GTt,p ≤ STt,j × GMAXt,j jJ t Start/stop process (SSTt,j) STt,j - STt-1,j ≤ SRt,j tT
The partial objective function of MIP-ED is Min t jJ ( CTj × GTt,j + CARRj × SRt,j ) C.
(46) (47)
(48)
Mixed-Non-Linear Economic Dispatch - Constraints ECONOMIC DISPATCH OF ELECTRIC SYSTEMS MIXED NON-LINEAR MODELS
CONSUMER
RESERVOIR
3
~ HYDRO-ELECTRIC HYDRO PLANTS - CONFIGURATION BUS
THERMO-ELECTRIC
1 12 8 6 1 24
Water Pumping Station Hydro Plants Reservoirs Thermal Plants Deficit Plant Hour Planning Horizon
The hydro-electric system includes a pumping water station to a reservoir; to simulate the operation, a non-linear equation determines the amount of energy needed for pumping, which is fixed on pressure and variable on flow. The MINLP-ED model includes the restrictions of LP-ED plus the water pumping process constraints. The additional equations are: Pumping availability (BPNt,i and BPXt,i) STt,j × GMINt,j ≤ GTt,p ≤ STt,j × GMAXt,j t , i Water balance including pumping (WAPt,i) Vt,i - Vt-1,i = At,i + pUP(i) (Qt,p + SPt,p) t , i – (Qt,i + SPt,i) + pRB(i) PUt,p - pRB(i) PUt,p Power use in pumping (POWt,i) PWt,i = 9.8 × 0.001 × pRB(i) PUt,i t iIPU Pumping power availability (PWNt,i and PWXt,i) PUmini × BPt,i ≤ pRB(i) PWt,i t , iIPU ≤ PUmaxi × BPt,i
The objective function of MINLP-ED is Min t jJ ( CTj × GTt,j + i CU1t,i × PWt,i + CU2t,i × PWt,i2)
(49) (50)
GTt,j
(51) (52)
(53)
VI. GDDP ECONOMIC DISPATCH CONCEPTUALIZATION First, to ensure that all solution generated for the state variables by the master problem always generates feasible solutions in the subproblems, it is required to include in the master problem a restriction on the maximum of hydroelectric generation as a function of the demand of the system. This restriction is: ▪ Maximum hydraulic generation (MGHt) GHt ≤ Lt (54) MGHt may be included as an operational bound. Next step is the partition of variables in state variables and control variables. In this case the state variables are reservoir level, hydro generation, water releases and state of thermal plants; the control variables are thermal generation. LP-ED, UC and MINLP-ED models must be decomposed into multiple problems, masters and subproblems, according the principles of GDDP. The problems are: ▪ M-LP-ED: Master of LP-ED ▪ SP-LP-ED(t): Subproblems of LP-ED, indexed in t ▪ M-MIP-ED: Master of UC ▪ SP-MIP-ED(t): Subproblems of UC ▪ M-MINLP-ED: Master of MINLP-ED ▪ SP-MINLP-ED(t): Subproblems of MINLP-ED The following tables show the variables and the constraints assigned to each problem.
UHt,i QHt,i GHt Vt,I Qt,i SPt,i GTt,j
UHt,i QHt,i GHt Vt,i Qt,i SPt,i SRt,j STt,j
LP-ED: Linear Economic Dispatch Variables Constraints G-SDDP Coordinator Hydro unit outflow QQHt,i Turbinated outflow Hydro unit power generation WATt,i Water balance Total hydraulic generation GQTt,i Hydraulic generation Reservoir operating volume GUTt,i Total hydraulic generation Turbined outflow MGHt Maximum hydraulic Spillage outflow generation G-SDDP Subproblem GUTt,i Total hydraulic generation Thermal generation ELDt,i Electricity demand
UHt,k QHt,k GHt Vt,I Qt,i SPt,i PUt,i,p PWt,i BPt,i
GTt,j
MIP-ED: MIP Economic Dispatch Variables Constraints G-SDDP Coordinator QQHt,i Turbinated outflow WATt,i Water balance Hydro unit outflow GQTt,i Hydraulic generation Hydro unit power generation GUTt,i Total hydraulic generation Total hydraulic generation MGHt Maximum hydraulic Reservoir operating volume generation Turbined outflow STNt,j Maximum thermal Spillage outflow availability Start thermal plant STXt,j Maximum thermal State thermal plant (off, on) availability SSTt,j Start/stop process G-SDDP Subproblem GUTt,i Total hydraulic generation Thermal generation ELDt,i Electricity demand MINLP-ED: Mixed-Non-Linear Economic Dispatch Variables Constraints G-SDDP Coordinator QQHt,i Turbinated outflow GQTt,i Hydraulic generation GUTt,i Total hydraulic generation MGHt Maximum hydraulic Hydro unit outflow generation Hydro unit power generation BPNt,i Minimum pumping Total hydraulic generation availability Reservoir operating volume BPXt,i Maximum pumping Turbined outflow availability Spillage outflow WAPt,i Water balance including Water pumped pumping Power used in pumping station POWt,i Power use in pumping Binary control pumping PWNt,i Minimum pumping power availability PWXt,i Maximum pumping power availability G-SDDP Subproblem Thermal generation GUTt,i Total hydraulic generation ELDt,i Electricity demand
VII. Economic Dispatch Results The models were implemented in GAMS-CPLEX. The implementation doesn’t use high level parallel processing to solve GDDP/G-SDDP subproblems. The experiments were realized in a HP-PC with a CPU i7-7500U (4 cores), 2.9 GHz and 12 MBytes of cache memory. A.
Stochastic Linear Economic Dispatch
Experiments with the stochastic linear economic dispatch (LPED) aim to prove the convergence of G-SDDP and compare the performance of: ▪ G-SDDP using Unified Benders Cuts (G-SDDP) ▪ SDDP ▪ G-SDDP using a DDP as coordinator (G-SDDP-DDP). FULL case corresponds to the deterministic equivalent of the stochastic optimization problem, solved without using large scale methodologies. Far al cases the number of periods is 24. Model FULL G-SDDP SDDP G-SDDP-DDP FULL G-SDDP SDDP G-SDDP-DDP
Scenarios
10
20
GAP
0.00000 0.00011 0.00000 0.00000 0.00017 0.00000
Solution Time (secs) 0.402 4.297 458.00 25.172 0.606 7.005 687.00 53.092
Times FULL 1.00 10.69 1139.30 62.61 1.00 11.56 1133.66 87.607
Times GSDDP
GAMS Options
1.00 106.59 5.85
SLINK SLINK+GUSS SLINK+GUSS
1.00 98.07 7.58
SLINK SLINK+GUSS SLINK+GUSS
The conclusions are: 1. The complexity doesn’t justify large-scale methodologies; but it is useable to compare G-SDDP and SDDP. 2. G-SDDP with unified cuts is approximately 100 times
faster than SDDP G-SDDP with an integrate coordinator is approximately 6 times faster than G-SDDP with a DDP coordinator The G-SDDP approach is more efficient than SDDP
3. 4. B.
Deterministic MIP Economic Dispatch
Additionally to LP-ED, the MIP Economic Dispatch (MIPED) model includes start/stop cost for thermo-electric plants. Experiments with the deterministic MIP-ED aim to: ▪ Prove the convergence of GDDP for mixed coordinator, using DBC and UBC. ▪ Evaluate the performance of GDDP. Methodology FULL GDDP-DBC GDDP-UBC FULL GDDP-DBC GDDP-UBC FULL GDDP-DBC GDDP-UBC FULL GDDP-DBC GDDP-UBC FULL GDDP-DBC GDDP-UBC FULL GDDP-DBC GDDP-UBC
Periods
Real Variables
Binary Variables
Cons traints
24
576
288
769
48
2214
384
1537
96
4225
1152
3073
192
9217
1536
6145
384
18433
3072
12289
768
36865
6144
24577
Solution Time (secs) 0.373 1.277 1.947 0.896 4.673 5.424 1.309 17.546 13.126 27.466 44.582 16.540 1001.147 75.094 33.032 NS 185.228 130.555
Time/ Scenario (sec/und) 0.016 0.053 0.081 0.019 0.097 0.113 0.014 0.183 0.137 0.143 0.232 0.086 2.607 0.196 0.086 NA 0.241 0.170
The conclusions are: 1. GDDP is convergent for any type of Benders cuts. 2. GDDP with Unified Benders Cuts is faster than GDDP with Decoupled Benders Cuts 3. In “large” cases the FULL problem doesn’t get a solution. C.
Stochastic MIP Economic Dispatch
Experiments with the stochastic MIP-ED aim to evaluate the performance of G-SDDP. Scenarios 20 50 100 200 500 1000 2000 4000 8000
Total Variables
Binaries Variables
Constraints
29,924 67,224 134,444 268,844 672,044 1,344,044 2,688,044 5,376,044 10,752,044
5,760 14,400 28,800 57,600 144,000 288,000 576,000 1,152,000 2,304,000
16,221 40,551 81,101 162,201 440,501 811,001 1,622,001 3,244,001 6,488,001
NO-Ceros
Solution Time
Time/ Scenario
69,521 173,801 347,601 695,201 1,738,001 3,476,001 6,952,001 13,904,001 27,808,001
(secs) 12.01 127 387 1092 1321 1570 2524 4634 10852
(secs) 0.601 2.537 3.867 5.461 2.642 1.570 1.262 1.159 1.356
The relation solution time versus number of periods-scenarios is not exponential, seems lineal, that implies that is possible to manage very large problems; additionally, the solution time per period indicates that a learning process occurs, then it decreases in the first part of the scenario axis. STOCHASTIC UNIT COMMITMENT MIXED MODEL – UNIFIED CUTS + INEXACT SOLUTIONS 12000
STOCHASTIC UNIT COMMITMENT MIXED MODEL – UNIFIED CUTS + INEXACT SOLUTIONS
6
TIME (sec)
10000
TIME (sec)
5
8000
4
6000
3
4000
2 1
2000
SCENARIOS
0
0
1000
2000
3000
4000
5000
6000
7000
8000
SCENARIOS
0 0
1000
2000
3000
4000
5000
6000
7000
8000
D.
Stochastic Mixed Non-Linear Economic Dispatch
Two stochastic non-linear dispatch models were considered: ▪ NLP-ED, includes all the variables and constraints related with water pumping but relaxed the minimum lower bound for flow in the pumps; than the problem is continuous. ▪ MINLP-ED, includes all the variables and constraints related with water pumping, but includes a lower bound (greater than cero) for flow, that required a binary variable. In both cases the number of periods is 24 and G-SDDP-UBC. Model Type
Scenarios
1
NLP-ED 20
50
MINLPED
1 5
Dual Bound 2718504 2718499 2718507 2718504 2290453 2290450 229045 I
Primal Bound 2718504 2718499 2718507 2718504 2290453 2290450 2290456
2256641 2256620 12998260 12998260
2256641 2256620 12998260 12998260
0 0 0 0
12302500 12302500
12302500 12302500
0 0
GAP 0 0 0 0 0 0 0
Solution Time (secs) 3.737 4.102 4.905 63.855 38.888 40.075 217.467 NS 95.831 102.706 9.475 24.160 NS 15.901 301.267
Solver QPC/MQPC CPLEX XPRESS MINOS IPOPT CPLEX XPRESS MINOS IPOPT CPLEX XPRESS CPLEX XPRESS BONMIN CPLEX XPRESS
The conclusion is that the G-SDDP is convergent for mixed non-linear convex problems. VIII. Conclusions The conceptual formulation of the GDDP and G-SDDP problem enables development of efficient algorithms based on the partition and the decomposition of the original problem using Benders' Theory and following the conceptualization of Dynamic Programing. The solution of the original problem is found by the coordinated solution of multiple problems of smaller dimension than other large-scale methodologies. In some cases, it is possible to visualize special matrix structures to generate Benders´ cuts for all periods and eliminates the need to solve a problem for each period. This is of special importance when the number of periods is very large, as in the case of discrete optimal control problems and stochastic optimization. Acknowledgments The author is grateful to Juan José Torres the work done in the implementation of the GDDP. This research has been financed by DecisionWare Ltda. and DO Analytics LLC. The GDDP/GSDDP are part of large-scale methodologies implemented in OPTEX Mathematical Modeling System (Velásquez, 2017). References Bellman, R.E. (1952), “On the Theory of Dynamic Programming”, Proceedings of the National Academy of Sciences, 1952 Benders, J. F. (1962) "Partitioning Procedures for Solving Mixed Variables Programming Problems". Numer. Math 4, 238-252. Birge, J. R., (1985). “Decomposition and Partitioning Methods for Multistage Stochastic Linear Programs”. Operations Research 33 (5), 989–1007. Diniz, A. L. (2010) “Test Cases for Unit Commitment and Hydrothermal Scheduling Problems”. Power and Energy Society General Meeting, 2010 IEEE Geoffrion, A. M. (1972) “Generalized Benders Decomposition”. Journal of Optimization Theory and Applications”, Vol. 10, No.4. Murphy, J. (2013). “Benders, Nested Benders and Stochastic Programming: An Intuitive Introduction”. Cambridge University
Engineering Department Technical Report CUED/FINFENG/TR.675, December 2013. Pereira, M. V. F. and Pinto, L. M. G. (1985) "Stochastic Optimization of a Multi-Reservoir Hydroelectric System: A Decomposition Approach". Water Resources Research, 21, 779-792. Pereira, M. V. F. and Pinto, L. M. G. (1991) "Multi-stage Stochastic Optimization Applied to Energy Planning". Mathematical Programming, 52, 359-375. Rahmaniania, R., Crainic, T. G., Gendreau, M., Rei, W. “The Benders Decomposition Algorithm: A Literature Review”, European Journal of Operational Research Volume 259, Issue 3, 16 June 2017. Velásquez, J. M. (1995) OEDM: Optimización Estocástica Dinámica Multinivel. Teoría General”. Universidad Nacional de Colombia, Revista Energética No. 13 (1995) http://www.doanalytics.net/Documents/DW-DT-033-01Optimizacion-Estocastica-Dinamica-Multinivel.pdf Velásquez, J. M. et al. (1999) “Dual Dynamic Programming: A note on Implementation”. Water Resources Research Vol 35. No. 7, 2269-2271, (July 1999).
Velásquez, J. M. (2002) “GDDP: Generalized Dual Dynamic Programming Theory”, Annals of Operation Research, November 2002, Volume 117, Issue 1–4, pp 21–31 Velásquez, J. M. (2017) “OPTEX: The First Automatic Generator of Large Scale Optimization Models”. DOA web-Conference http://www.doanalytics.net/Documents/OPTEX-Presentation.pdf Velásquez, J. M. (2018a) “Benders Decomposition Using Unified Cuts”, White Paper DW-DT-047-2018, http://www.doanalytics.net/Documents/Benders-DecompositionUsing-Unified-Cuts.pdf Velásquez, J. M. (2018b) “GDDP: Generalized Dual Dynamic Programming Theory & Electric Sector Applications”, DW-DWDT-048-2018, http://www.doanalytics.net/Documents/GDDPGeneralized-Dual-Dynamic-Programming-Theory-&-ElectricSector-Applications.pdf Velásquez, J. M. (2018c) “G-SDDP: Generalized Stochastic Dual Dynamic Programming Theory & Electric Sector Applications”, White Paper DW-DT-049-2018, http://www.doanalytics.net/Documents/G-SDDPGeneralizedStochastic-Dual-Dynamic-Programming-Theory-&-ElectricSector-Applications.pdf
