# An Improved Lagrangian Relaxation Algorithm for the Robust Generation Self-Scheduling Problem

Authors: Ping Che; Zhenhao Tang; Hua Gong; Xiaoli Zhao
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

The robust generation self-scheduling problem under electricity price uncertainty is usually solved by the commercial solver, which is limited in computation time and memory requirement. This paper proposes an improved Lagrangian relaxation algorithm for the robust generation self-scheduling problem where the quadratic fuel cost and the time-dependent exponential startup cost are considered. By using the optimal duality theory, the robust generation self-scheduling problem, which has a max-min structure, is reformulated as a minimization mixed integer nonlinear programming (MINLP) problem. Upon the reformulation, the Lagrangian relaxation algorithm is developed. To obtain a solvable relaxed problem, the variable splitting technique is introduced before the relaxation. The obtained relaxed problem is decomposed into a linear programming-type subproblem and multiple single-unit subproblems. Each single-unit subproblem is solved optimally by a two-stage backward dynamic programming procedure. The special cases of the problem are discussed and a two-stage algorithm is proposed. The proposed algorithms are tested on test cases of different sizes and the numerical results show that the algorithms can find near-optimal solutions in a reasonable time.


---

## Body

## 1. Introduction

In the competitive electricity market, the generation self-scheduling problem plays a key role in the planning and operation of electric power systems. An effective generation self-scheduling can help the power generation company to decide the bidding strategies and make generation plans for the subsequent time horizon. The generation self-scheduling problem is to determine the operation of generation units, including the on/off statuses of generation units and the output level of committed generation units, according to electricity prices over the scheduling horizon. The objective is to maximize the total generation profit. Because the formulation of the generation self-scheduling problem includes quadratic fuel cost function, exponential startup cost function, and integer variables indicating the on/off statuses of generation units, the problem is a nonconvex MINLP problem. As there could be a large number of generation units and a long scheduling horizon, the generation self-scheduling problem can become large-scale and computationally challenging. This motivates us to study the solution method of the generation self-scheduling problem.The decision-making of the generation self-scheduling is affected by electricity prices, which are volatile and full of uncertainty in the practical spot market. The uncertainty may be caused by the volatility of fuel prices and the changes in government subsidies, industry regulations, and local weather conditions. Therefore, electricity prices cannot be forecasted accurately and need to be treated as uncertain parameters. Two approaches are usually applied in dealing with uncertain parameters: stochastic programming and robust optimization. The stochastic programming approach relies on the probability distribution of uncertain parameters and typically expresses the possible realizations of uncertain parameters by a finite number of scenarios according to the probability distribution. Examples of using the approach can be seen in [1–7]. The drawback of the approach is that it is difficult to obtain the exact distribution of uncertain parameters in the practice. The robust optimization approach is distribution-free. The approach represents uncertain parameters with a deterministic uncertainty set that contains all possible realizations or an adequate realization range of uncertain parameters. The obtained decision is robust against the variation of uncertain parameters within the uncertainty set. Examples of using the approach can be seen in [8–13]. To make the decision more practical, in this paper, we use the robust optimization approach to model uncertain electricity prices in the generation self-scheduling problem.The solution method of the generation self-scheduling problem has been widely investigated, with the majority concentrated on the deterministic generation self-scheduling problem and the stochastic generation self-scheduling problem. For the deterministic generation self-scheduling problem, the mixed integer linear programming approach was used in [14–16], a Lagrangian relaxation algorithm was proposed in [17], the Lagrangian relaxation method and the MIP method were compared in [18], a particle swarm optimization algorithm was proposed in [19], a genetic algorithm was proposed in [20], a new memetic algorithm was proposed in [21], an ant colony optimization algorithm was proposed in [22], an evolutionary algorithm was proposed in [23], and a survey of solution methods was presented in [24]. For the stochastic generation self-scheduling problem where uncertain electricity prices are modeled by the stochastic programming approach, the mixed integer linear programming approach was used in [1–4], the typical Lagrangian relaxation algorithm was applied in [5], and a sample average approximation algorithm was proposed in [6].For the robust generation self-scheduling problem where uncertain electricity prices are modeled by the robust optimization approach, the solution algorithm is not widely investigated. A critical review of the robust generation self-scheduling was presented in [25]. To provide a better sketch for the research of the robust generation self-scheduling problem, a comparison of some relevant literature is presented in Table 1. From Table 1, we can observe that the robust generation self-scheduling problem was usually solved by commercial solvers. In order to make the problem solvable by the solvers, the model was simplified by piecewise linearly approximating the quadratic fuel cost function [8, 10, 11], reducing the time-dependent exponential startup cost function as a constant [8, 10], omitting the startup cost [9, 12], or omitting the unit commitment decision [12]. The solution approach has two disadvantages. First, the use of commercial solvers is limited in computation time and memory requirement [18]. When the robust generation self-scheduling problem is large-scale, it is not practical to solve the problem by commercial solvers. Second, the simplification of the model will reduce the effect of the solution. Therefore, an effective solution algorithm should be developed to solve the robust generation self-scheduling problem.Table 1
A comparison of some relevant literature.
No. Price Robust Thermal Unit Commitment Quadratic Fuel Cost Exponential Startup Costa Solution Reference 1 Y Y N N MILP Solver [8] 2 Y Y Y - MIQP Solver [9] 3 Y Y N N MILP Solver [10] 4 Y Y N N MILP Solver [11] 5 Y N Y - QP Solver [12] 6 Y Y Y Y LRb Proposed a “-” indicates that the startup cost was not considered. b Lagrangian relaxation algorithm.In this paper, we develop a Lagrangian relaxation algorithm for the robust generation self-scheduling problem. Compared with the existing robust generation self-scheduling literature, this paper provides the following contributions:(1) We consider the quadratic fuel cost and the exponential startup cost in the robust generation self-scheduling problem and develop an effective Lagrangian relaxation algorithm to solve the problem. The algorithm decomposes the problem into a linear programming-type subproblem and multiple single-unit subproblems. Each single-unit subproblem is solved by a two-stage backward dynamic programming procedure. The feasible solution is constructed by a heuristic algorithm.(2) Unlike the typical Lagrangian relaxation algorithm that may lead to an unbounded relaxed problem for the considered problem, we introduce the variable splitting technique to improve the algorithm.(3) A numerical comparison between the proposed algorithm and the MILP solver is reported. Numerical results demonstrate the effectiveness of the proposed algorithm.The rest of the paper is organized as follows. The mathematical formulation of the robust generation self-scheduling problem is presented in Section2. The solution methodology is proposed in Section 3. Numerical experiments are carried out in Section 4. The research is concluded in Section 5.


## 2. Problem Description and Formulation

We consider the robust generation self-scheduling problem under electricity price uncertainty from the viewpoint of a price-taking power producer. The objective is to maximize the total generation profit under the worst-case scenario within the considered range of electricity prices. The fuel cost is a quadratic function of the power generation level and the startup cost is an exponential function of the down time of the generation unit. The schedule needs to satisfy generation unit operation constraints and system electricity demand constraints. We take an hour as the scheduling period in the generation self-scheduling problem under consideration.The remainder of the section is organized as follows. The notations used in the problem formulation are presented in Section2.1. The modeling of uncertain electricity prices is described in Section 2.2. An MINLP model for the robust generation self-scheduling problem is formulated in Section 2.3.


### 2.1. Notations

Parameters n: Number of generation units T: Scheduling horizon R D i / R U i: Ramp-down/ramp-up rate of generation unit i S D i / S U i: Shutdown/startup ramp rate of generation unit i P i L / P i U: Minimum/maximum power output of generation unit i L i / U i: Minimum down-time/up-time of generation unit i D t: Electricity demand in hour t F C i ( · , · ): Fuel cost function of generation uniti a 0 i , a 1 i , a 2 i: Coefficients of the quadratic fuel cost function of generation uniti S C i ( · , · ): Startup cost function of generation uniti b 1 i , b 2 i , τ i: Coefficients of the startup cost function of generation uniti λ t L / λ t U: Lower/upper bound of the electricity price in hourt λ t: Random parameter representing the electricity price in hourt Γ: Budget lower bound of the sum of the electricity price λt within the scheduling horizon y i 0: Number of hours generation uniti being up or down at the end of hour 0 u i 0: Binary parameter to indicate the initial on/off status of generation uniti p i 0: Power generation level of generation uniti in hour 0 t i d: Last hour of the time periods during which the on/off statuses of generation uniti must be the same to its initial on/off status t i c: Earliest hour that generation uniti can be committed in the scheduling horizon γ: Weight for the artificial variable introduced in the algorithmDecision Variables u i t: Binary variable to indicate the on/off status of generation uniti in hourt v i t: Binary variable to indicate if generation uniti is started up in hourt y i t: State variable to indicate the number of hours generation uniti being up or down at the end of hourt p i t: Power generation level of generation uniti in hourt ω 0 , ω t: Dual variables z t: Introduced variable in the solution method q t: Artificial variable ξ 1 i t , ξ 2 i t , μ t , η t: Lagrangian multipliers


### 2.2. Modeling of Uncertain Electricity Prices

According to the robust optimization approach, we model uncertain electricity prices with the following uncertainty set:Λ ≔ λ 1 , … , λ T:(1a)λtL≤λt≤λtU,t=1,…,T(1b)∑t=1Tλt≥Γ.Constraints (1a) restrict the electricity price for each hour between a lower bound and an upper bound, which can be set to the 2.5% and 97.5% quantiles of the electricity price forecast, respectively. Constraint (1b) presents a budget lower bound Γ for the sum of the electricity price λt over the scheduling horizon. The budget value Γ is used to control the level of conservatism of the robust optimization approach. The smaller the value of Γ is, the more conservative the approach is. System operators can choose the budget value according to their requirements.Within the above uncertainty set, there are infinite electricity price scenarios. For a fixed generation self-schedule, the best-case electricity price scenario is the scenario corresponding to the maximum electricity sales revenue, while the worst-case electricity price scenario is the scenario corresponding to the minimum electricity sales revenue. It can be observed that the best-case electricity price scenario is atλt=λtU for all t for any fixed power generation level, while the worst-case electricity price scenario is related to the fixed power generation level and the budget value Γ [26].


### 2.3. The Robust Generation Self-Scheduling Model

The robust generation self-scheduling problem for a price-taking power producer, denoted by (RSS), is formulated as follows:(2)maxu,pminλ∈Λ∑t=1Tλt∑i=1npit-∑i=1n∑t=ticTFCipit,uit+SCiyi,t-1,vit(3)s.t.uit=ui0,i=1,…,n,t=1,…,tid(4)vit=0,i=1,…,n,t=1,…,tid(5)uit-ui,t-1≤uik,i=1,…,n,k=t+1,…,min⁡T,Ui+t-1,t=tid+1,…,T(6)ui,t-1-uit≤1-uik,i=1,…,n,k=t+1,…,min⁡T,Li+t-1,t=tid+1,…,T(7)uit-ui,t-1≤vit,i=1,…,n,t=tid+1,…,T(8)yit=uitmax⁡1,yi,t-1+1+1-uitmin⁡-1,yi,t-1-1,i=1,…,n,t=1,…,T(9)uitPiL≤pit≤uitPiU,i=1,…,n,t=1,…,T(10a)-RDiui,t+1-SDi1-ui,t+1≤pi,t+1-pit≤RUiuit+SUi1-uit,i=1,…,n,t=tic-1(10b)-RDiui,t+1-SDi1-ui,t+1≤pi,t+1-pit≤RUiuit+SUi1-uit,i=1,…,n,t=tic,…,T-1(11)∑i=1npit≤Dt,t=1,…,T(12)pit≥0,i=1,…,n,t=1,…,T(13)uit,vit∈0,1,i=1,…,n,t=1,…,T(14)yit:nonzerointeger,i=1,…,n,t=1,…,TIn the above formulation, the objective function (2) is to maximize the generation profit under the worst-case scenario within the uncertainty set of electricity prices. Under the max-min decision rule, the obtained generation schedule is robust against the variation of electricity prices within the uncertainty set. The generation profit is determined by the electricity sales revenue, the fuel cost, and the startup cost over the scheduling horizon where the fuel cost is FCi(pit,uit)=a0iuit+a1ipit+a2ipit2 and the startup cost is SCi(yi,t-1,vit)=[b1i(1-exp⁡(yi,t-1/τi))+b2i]vit according to [27]. Equations (3) and (4) show the impact of the initial statuses of units on decision-making. Inequalities (5) and (6) represent the minimum up-time and down-time requirements, respectively. Inequalities (7) represent the startup status. Equations (8) represent the relationship between the state variables in adjacent hours. Inequalities (9) represent the power generation capacity of units. Inequalities (10a) and (10b) represent the ramping-up and ramping-down rate limits. Inequalities (11) represent the electricity demand constraints. Constraints (12)-(14) show the value field of the decision variables.


## 2.1. Notations

Parameters n: Number of generation units T: Scheduling horizon R D i / R U i: Ramp-down/ramp-up rate of generation unit i S D i / S U i: Shutdown/startup ramp rate of generation unit i P i L / P i U: Minimum/maximum power output of generation unit i L i / U i: Minimum down-time/up-time of generation unit i D t: Electricity demand in hour t F C i ( · , · ): Fuel cost function of generation uniti a 0 i , a 1 i , a 2 i: Coefficients of the quadratic fuel cost function of generation uniti S C i ( · , · ): Startup cost function of generation uniti b 1 i , b 2 i , τ i: Coefficients of the startup cost function of generation uniti λ t L / λ t U: Lower/upper bound of the electricity price in hourt λ t: Random parameter representing the electricity price in hourt Γ: Budget lower bound of the sum of the electricity price λt within the scheduling horizon y i 0: Number of hours generation uniti being up or down at the end of hour 0 u i 0: Binary parameter to indicate the initial on/off status of generation uniti p i 0: Power generation level of generation uniti in hour 0 t i d: Last hour of the time periods during which the on/off statuses of generation uniti must be the same to its initial on/off status t i c: Earliest hour that generation uniti can be committed in the scheduling horizon γ: Weight for the artificial variable introduced in the algorithmDecision Variables u i t: Binary variable to indicate the on/off status of generation uniti in hourt v i t: Binary variable to indicate if generation uniti is started up in hourt y i t: State variable to indicate the number of hours generation uniti being up or down at the end of hourt p i t: Power generation level of generation uniti in hourt ω 0 , ω t: Dual variables z t: Introduced variable in the solution method q t: Artificial variable ξ 1 i t , ξ 2 i t , μ t , η t: Lagrangian multipliers


## 2.2. Modeling of Uncertain Electricity Prices

According to the robust optimization approach, we model uncertain electricity prices with the following uncertainty set:Λ ≔ λ 1 , … , λ T:(1a)λtL≤λt≤λtU,t=1,…,T(1b)∑t=1Tλt≥Γ.Constraints (1a) restrict the electricity price for each hour between a lower bound and an upper bound, which can be set to the 2.5% and 97.5% quantiles of the electricity price forecast, respectively. Constraint (1b) presents a budget lower bound Γ for the sum of the electricity price λt over the scheduling horizon. The budget value Γ is used to control the level of conservatism of the robust optimization approach. The smaller the value of Γ is, the more conservative the approach is. System operators can choose the budget value according to their requirements.Within the above uncertainty set, there are infinite electricity price scenarios. For a fixed generation self-schedule, the best-case electricity price scenario is the scenario corresponding to the maximum electricity sales revenue, while the worst-case electricity price scenario is the scenario corresponding to the minimum electricity sales revenue. It can be observed that the best-case electricity price scenario is atλt=λtU for all t for any fixed power generation level, while the worst-case electricity price scenario is related to the fixed power generation level and the budget value Γ [26].


## 2.3. The Robust Generation Self-Scheduling Model

The robust generation self-scheduling problem for a price-taking power producer, denoted by (RSS), is formulated as follows:(2)maxu,pminλ∈Λ∑t=1Tλt∑i=1npit-∑i=1n∑t=ticTFCipit,uit+SCiyi,t-1,vit(3)s.t.uit=ui0,i=1,…,n,t=1,…,tid(4)vit=0,i=1,…,n,t=1,…,tid(5)uit-ui,t-1≤uik,i=1,…,n,k=t+1,…,min⁡T,Ui+t-1,t=tid+1,…,T(6)ui,t-1-uit≤1-uik,i=1,…,n,k=t+1,…,min⁡T,Li+t-1,t=tid+1,…,T(7)uit-ui,t-1≤vit,i=1,…,n,t=tid+1,…,T(8)yit=uitmax⁡1,yi,t-1+1+1-uitmin⁡-1,yi,t-1-1,i=1,…,n,t=1,…,T(9)uitPiL≤pit≤uitPiU,i=1,…,n,t=1,…,T(10a)-RDiui,t+1-SDi1-ui,t+1≤pi,t+1-pit≤RUiuit+SUi1-uit,i=1,…,n,t=tic-1(10b)-RDiui,t+1-SDi1-ui,t+1≤pi,t+1-pit≤RUiuit+SUi1-uit,i=1,…,n,t=tic,…,T-1(11)∑i=1npit≤Dt,t=1,…,T(12)pit≥0,i=1,…,n,t=1,…,T(13)uit,vit∈0,1,i=1,…,n,t=1,…,T(14)yit:nonzerointeger,i=1,…,n,t=1,…,TIn the above formulation, the objective function (2) is to maximize the generation profit under the worst-case scenario within the uncertainty set of electricity prices. Under the max-min decision rule, the obtained generation schedule is robust against the variation of electricity prices within the uncertainty set. The generation profit is determined by the electricity sales revenue, the fuel cost, and the startup cost over the scheduling horizon where the fuel cost is FCi(pit,uit)=a0iuit+a1ipit+a2ipit2 and the startup cost is SCi(yi,t-1,vit)=[b1i(1-exp⁡(yi,t-1/τi))+b2i]vit according to [27]. Equations (3) and (4) show the impact of the initial statuses of units on decision-making. Inequalities (5) and (6) represent the minimum up-time and down-time requirements, respectively. Inequalities (7) represent the startup status. Equations (8) represent the relationship between the state variables in adjacent hours. Inequalities (9) represent the power generation capacity of units. Inequalities (10a) and (10b) represent the ramping-up and ramping-down rate limits. Inequalities (11) represent the electricity demand constraints. Constraints (12)-(14) show the value field of the decision variables.


## 3. Solution Methodology

As the objective function (2) is max-min-type and includes nonlinear cost functions, (RSS) is a max-min MINLP problem. To solve the problem, we first reformulate (RSS) as a minimization model using the idea provided in [26] and then develop an improved Lagrangian relaxation algorithm upon the reformulation.


### 3.1. Reformulation of (RSS)

According to constraints (1a), let(15)λt=λtU-λtU-λtLzt,t=1,…,Twhere zt,  t=1,…,T, are introduced variables satisfying(16)0≤zt≤1,t=1,…,T.Then we can reformulate minλ∈Λ⁡∑t=1Tλt∑i=1npit as (17)min∑t=1TλtU∑i=1npit-∑t=1TλtU-λtLzt∑i=1npit(18)s.t.∑t=1TλtU-λtLzt≤∑t=1TλtU-Γ.By dualizing constraints (16) and (18), we can transform (16)-(18) as follows:(19)max∑t=1TλtU∑i=1npit-∑t=1TλtU-Γω0-∑t=1Tωt(20)s.t.λtU-λtLω0+ωt≥λtU-λtL∑i=1npit,t=1,…,T(21)ωt≥0,t=1,…,T(22)ω0≥0where ωt, t=1,…,T, and ω0 are the dual variables for constraints (16) and (18), respectively.Based on the above transformation, (RSS) can be reformulated as(RSS1)(23)max∑t=1TλtU∑i=1npit-∑i=1n∑t=ticTFCipit,uit+SCiyi,t-1,vit-∑t=1TλtU-Γω0-∑t=1Tωts.t. constraints (3)-(14), and (20)-(22)or equivalently(RSS2)(24)min∑i=1n∑t=ticTFCipit,uit+SCiyi,t-1,vit-∑t=1TλtU∑i=1npit+∑t=1TλtU-Γω0+∑t=1Tωts.t. the same constraints in (RSS1).


### 3.2. The Lagrangian Relaxation Algorithm

(  1) Variable Splitting-Based Lagrangian Relaxation. When the Lagrangian relaxation algorithm is used to solve a generation scheduling problem, constraints that couple different units are typically relaxed to make the relaxed problem separable in units. In this paper, constraints that couple different units include constraints (11) and (20). If constraints (20) are relaxed, the resulting subproblem that contains dual variables ω0 and ωt,t=1,…,T, will be as follows:(25)min∑t=1TλtU-Γ-∑t=1TηtλtU-λtLω0+∑t=1T1-ηtωts.t. constraints (21) and (22)whereηt,  t=1,…,T, are nonnegative Lagrangian multipliers. Note that the above subproblem is unbounded when either ∑t=1TλtU-Γ-∑t=1TηtλtU-λtL or 1-ηt is negative. Consequently, using the typical Lagrangian relaxation algorithm may lead to an unbounded relaxed problem.To obtain a bounded relaxed problem, we reserve constraints (20) and introduce the variable splitting technique into the algorithm. The variable splitting technique is to duplicate certain decision variables by adding some artificial variables and variable copy constraints that link the decision variables and the added artificial variables into the problem [28]. This technique is usually used to obtain a stronger lower bound for a minimization problem [29–33], but we use it to construct a solvable relaxed problem in this paper. Based on the technique, we add artificial variables qt,  t=1,…,T, and the following variable copy constraints:(26)∑i=1npit=qt,t=1,…,Tinto (RSS2). We also introduce a parameter γ∈0,1 to represent the weight of the added artificial variables in the objective function. The resulting problem, denoted by (RSS3), is equivalent to (RSS2) and described as follows:(27)min∑i=1n∑t=ticTFCipit,uit+SCiyi,t-1,vit-∑t=1TλtU1-γ∑i=1npit+γqt+∑t=1TλtU-Γω0+∑t=1Tωts.t. constraints (3)-(10a), (10b), (12)-(14), (21), (22), (26)(28)λtU-λtLω0+ωt≥λtU-λtLqt,t=1,…,T(29)qt≤Dt,t=1,…,T(30)qt≥0,t=1,…,Twhere the value of γ will be discussed in Section 4.We relax constraints (10b) and (26) and incorporate them into the objective function (27) by introducing Lagrangian multipliers {ξ1it≥0,ξ2it≥0},  i=1,…,n,  t=tic,…,T-1,  πt≥0,  t=1,…,T,  and  μt∈R,  t=1,…,T. The obtained relaxed problem, denoted by (RP), is as follows:(31)min∑i=1n∑t=ticTFCipit,uit+SCiyi,t-1,vit-∑t=1TλtU1-γ∑i=1npit+γqt+∑t=1TλtU-Γω0+∑t=1Tωt+∑i=1n∑t=ticT-1ξ1itpi,t+1-pit-RUiuit-SUi1-uit+ξ2itpit-pi,t+1-RDiui,t+1-SDi1-ui,t+1+∑t=1Tμt∑i=1npit-qts.t. constraints (3)-(10a), (12)-(14), (21), (22), and (28)-(30).(  2) Solution of the Relaxed Problem. Given the Lagrangian multipliers, (RP) can be decomposed into n+1 independent subproblems which are denoted by (RP1) and (RP2i), i=1,…,n, and expressed as follows:(RP1)(32)min∑t=1TλtU-Γω0+∑t=1Tωt-∑t=1TγλtU+μtqts.t. constraints (21), (22), and (28)-(30).(RP2i)(33)min∑t=ticTFCipit,uit+μt-1-γλtUpit+SCiyi,t-1,vit+ξ2itic-ξ1iticpitic+ξ1iticSUi-RUiuitic+∑t=tic+1T-1ξ1i,t-1-ξ1it+ξ2it-ξ2i,t-1pit+ξ1itSUi-RUi+ξ2i,t-1SDi-RDiuit+ξ1i,T-1-ξ2i,T-1piT+ξ2i,T-1SDi-RDiuiTs.t. constraints (3)-(10a) and (12)-(14).(RP1) is a linear programming model and can be solved optimally by using a commercial solver. Specially, ifΓ=∑t=1TλtU, we can achieve the following optimal solution directly without calling the commercial solver: (34)qt=Dt,ifγλtU+μt≥0,0,otherwise,t=1,…,T,(35)ω0=anyvaluenolessthanmax1≤t≤Tqt,(36)ωt=0,t=1,…,T.Each (RP2i) corresponds a single-unit subproblem. According to the feature of the subproblem, we solve it using a two-stage backward dynamic programming procedure [34, 35]. In the first stage, define Cit(pit) as the generalized generation cost of unit i in hour t if unit i is on in hour t in (RP2i) and determine the corresponding optimal generation level (37)p^it=arg⁡min⁡Citpit:pi,t-1-RDi≤pit≤pi,t-1+RUiui,t-1+SUi1-ui,t-1,PiL≤pit≤PiU,ift=tic,arg⁡min⁡Citpit:PiL≤pit≤PiU,ift>tic.In the second stage, determine the on/off statuses of the unit using the backward dynamic programming. The state transition equations are (38)fi,T+1yiT=0,fityi,t-1=minuit,vit∈Ωyi,t-1SCiyi,t-1,vit+Citp^ituit+fi,t+1yit,t=T,T-1,…,tic,where yit is subject to constraint (8) and fit(yi,t-1) is the optimal value function.The obtained relaxation solution can provide a lower bound for the optimal objective function value of (RSS3).(  3) Construction of the Feasible Solution. Because constraints (10b) and (26) are relaxed, the optimal solution to the relaxed problem is generally infeasible for (RSS3). To obtain a feasible solution, we propose the following heuristic method based on the current relaxation solution.Step 0.
Initializet=1.Step 1.
For each uniti, adjust pit according to pi,t-1 to meet constraints (10b) corresponding to hourt without violating constraints (9).Step 2.
Decreasepit for some units or shut off some units to meet constraint (11) corresponding to hourt without violating constraints (5), (6), (9), and (10b).Step 3.
Ift≥T, go to Step 4. Otherwise, set t=t+1 and go to Step 1.Step 4.
Check if constraints (20) are satisfied. If constraints (20) are satisfied, stop. Otherwise, solve the following linear programming problem to determine ω0 and ωt, t=1,…,T.(39)min∑t=1TλtU-Γω0+∑t=1Tωts.t. constraints (20)-(22).The obtained feasible solution can provide an upper bound for the optimal objective function value of (RSS3).(  4) Updating of the Lagrangian Multipliers. We initialize the Lagrangian multipliers by zero and update them according to the subgradient algorithm [36]. The iteration is stopped when the maximum number of iterations is reached or the relative duality gap (ZU-ZL)/ZL×100% is smaller than a certain threshold, where ZL and ZU are the best lower bound and upper bound obtained so far for (RSS3), respectively.


### 3.3. Discussions on Special Cases

(RSS) can be simplified in the following two special cases:Γ=∑t=1TλtL and Γ=∑t=1TλtU. In the case of Γ=∑t=1TλtL, we have λt=λtL for all t according to constraints (1a) and (1b) and (RSS) can be simplified as a deterministic multi-unit generation self-scheduling problem. Similar result can be obtained in the case of Γ=∑t=1TλtU. Therefore, we propose the following two-stage algorithm for the two special cases:Stage 1.
Simplify (RSS) by determiningλt,  t=1,…,T, according to the budget value Γ. If Γ=∑t=1TλtL, then λt=λtL,  t=1,…,T. If Γ=∑t=1TλtU, then λt=λtU,  t=1,…,T.Stage 2.
Solve the simplified problem by using the typical Lagrangian relaxation algorithm where constraints (10b) and (11) are relaxed.Since both the variable splitting-based Lagrangian relaxation algorithm and the two-stage algorithm proposed above can solve (RSS) in the two special cases, we will choose the more effective one for each special case based on the numerical test in Section4.


## 3.1. Reformulation of (RSS)

According to constraints (1a), let(15)λt=λtU-λtU-λtLzt,t=1,…,Twhere zt,  t=1,…,T, are introduced variables satisfying(16)0≤zt≤1,t=1,…,T.Then we can reformulate minλ∈Λ⁡∑t=1Tλt∑i=1npit as (17)min∑t=1TλtU∑i=1npit-∑t=1TλtU-λtLzt∑i=1npit(18)s.t.∑t=1TλtU-λtLzt≤∑t=1TλtU-Γ.By dualizing constraints (16) and (18), we can transform (16)-(18) as follows:(19)max∑t=1TλtU∑i=1npit-∑t=1TλtU-Γω0-∑t=1Tωt(20)s.t.λtU-λtLω0+ωt≥λtU-λtL∑i=1npit,t=1,…,T(21)ωt≥0,t=1,…,T(22)ω0≥0where ωt, t=1,…,T, and ω0 are the dual variables for constraints (16) and (18), respectively.Based on the above transformation, (RSS) can be reformulated as(RSS1)(23)max∑t=1TλtU∑i=1npit-∑i=1n∑t=ticTFCipit,uit+SCiyi,t-1,vit-∑t=1TλtU-Γω0-∑t=1Tωts.t. constraints (3)-(14), and (20)-(22)or equivalently(RSS2)(24)min∑i=1n∑t=ticTFCipit,uit+SCiyi,t-1,vit-∑t=1TλtU∑i=1npit+∑t=1TλtU-Γω0+∑t=1Tωts.t. the same constraints in (RSS1).


## 3.2. The Lagrangian Relaxation Algorithm

(  1) Variable Splitting-Based Lagrangian Relaxation. When the Lagrangian relaxation algorithm is used to solve a generation scheduling problem, constraints that couple different units are typically relaxed to make the relaxed problem separable in units. In this paper, constraints that couple different units include constraints (11) and (20). If constraints (20) are relaxed, the resulting subproblem that contains dual variables ω0 and ωt,t=1,…,T, will be as follows:(25)min∑t=1TλtU-Γ-∑t=1TηtλtU-λtLω0+∑t=1T1-ηtωts.t. constraints (21) and (22)whereηt,  t=1,…,T, are nonnegative Lagrangian multipliers. Note that the above subproblem is unbounded when either ∑t=1TλtU-Γ-∑t=1TηtλtU-λtL or 1-ηt is negative. Consequently, using the typical Lagrangian relaxation algorithm may lead to an unbounded relaxed problem.To obtain a bounded relaxed problem, we reserve constraints (20) and introduce the variable splitting technique into the algorithm. The variable splitting technique is to duplicate certain decision variables by adding some artificial variables and variable copy constraints that link the decision variables and the added artificial variables into the problem [28]. This technique is usually used to obtain a stronger lower bound for a minimization problem [29–33], but we use it to construct a solvable relaxed problem in this paper. Based on the technique, we add artificial variables qt,  t=1,…,T, and the following variable copy constraints:(26)∑i=1npit=qt,t=1,…,Tinto (RSS2). We also introduce a parameter γ∈0,1 to represent the weight of the added artificial variables in the objective function. The resulting problem, denoted by (RSS3), is equivalent to (RSS2) and described as follows:(27)min∑i=1n∑t=ticTFCipit,uit+SCiyi,t-1,vit-∑t=1TλtU1-γ∑i=1npit+γqt+∑t=1TλtU-Γω0+∑t=1Tωts.t. constraints (3)-(10a), (10b), (12)-(14), (21), (22), (26)(28)λtU-λtLω0+ωt≥λtU-λtLqt,t=1,…,T(29)qt≤Dt,t=1,…,T(30)qt≥0,t=1,…,Twhere the value of γ will be discussed in Section 4.We relax constraints (10b) and (26) and incorporate them into the objective function (27) by introducing Lagrangian multipliers {ξ1it≥0,ξ2it≥0},  i=1,…,n,  t=tic,…,T-1,  πt≥0,  t=1,…,T,  and  μt∈R,  t=1,…,T. The obtained relaxed problem, denoted by (RP), is as follows:(31)min∑i=1n∑t=ticTFCipit,uit+SCiyi,t-1,vit-∑t=1TλtU1-γ∑i=1npit+γqt+∑t=1TλtU-Γω0+∑t=1Tωt+∑i=1n∑t=ticT-1ξ1itpi,t+1-pit-RUiuit-SUi1-uit+ξ2itpit-pi,t+1-RDiui,t+1-SDi1-ui,t+1+∑t=1Tμt∑i=1npit-qts.t. constraints (3)-(10a), (12)-(14), (21), (22), and (28)-(30).(  2) Solution of the Relaxed Problem. Given the Lagrangian multipliers, (RP) can be decomposed into n+1 independent subproblems which are denoted by (RP1) and (RP2i), i=1,…,n, and expressed as follows:(RP1)(32)min∑t=1TλtU-Γω0+∑t=1Tωt-∑t=1TγλtU+μtqts.t. constraints (21), (22), and (28)-(30).(RP2i)(33)min∑t=ticTFCipit,uit+μt-1-γλtUpit+SCiyi,t-1,vit+ξ2itic-ξ1iticpitic+ξ1iticSUi-RUiuitic+∑t=tic+1T-1ξ1i,t-1-ξ1it+ξ2it-ξ2i,t-1pit+ξ1itSUi-RUi+ξ2i,t-1SDi-RDiuit+ξ1i,T-1-ξ2i,T-1piT+ξ2i,T-1SDi-RDiuiTs.t. constraints (3)-(10a) and (12)-(14).(RP1) is a linear programming model and can be solved optimally by using a commercial solver. Specially, ifΓ=∑t=1TλtU, we can achieve the following optimal solution directly without calling the commercial solver: (34)qt=Dt,ifγλtU+μt≥0,0,otherwise,t=1,…,T,(35)ω0=anyvaluenolessthanmax1≤t≤Tqt,(36)ωt=0,t=1,…,T.Each (RP2i) corresponds a single-unit subproblem. According to the feature of the subproblem, we solve it using a two-stage backward dynamic programming procedure [34, 35]. In the first stage, define Cit(pit) as the generalized generation cost of unit i in hour t if unit i is on in hour t in (RP2i) and determine the corresponding optimal generation level (37)p^it=arg⁡min⁡Citpit:pi,t-1-RDi≤pit≤pi,t-1+RUiui,t-1+SUi1-ui,t-1,PiL≤pit≤PiU,ift=tic,arg⁡min⁡Citpit:PiL≤pit≤PiU,ift>tic.In the second stage, determine the on/off statuses of the unit using the backward dynamic programming. The state transition equations are (38)fi,T+1yiT=0,fityi,t-1=minuit,vit∈Ωyi,t-1SCiyi,t-1,vit+Citp^ituit+fi,t+1yit,t=T,T-1,…,tic,where yit is subject to constraint (8) and fit(yi,t-1) is the optimal value function.The obtained relaxation solution can provide a lower bound for the optimal objective function value of (RSS3).(  3) Construction of the Feasible Solution. Because constraints (10b) and (26) are relaxed, the optimal solution to the relaxed problem is generally infeasible for (RSS3). To obtain a feasible solution, we propose the following heuristic method based on the current relaxation solution.Step 0.
Initializet=1.Step 1.
For each uniti, adjust pit according to pi,t-1 to meet constraints (10b) corresponding to hourt without violating constraints (9).Step 2.
Decreasepit for some units or shut off some units to meet constraint (11) corresponding to hourt without violating constraints (5), (6), (9), and (10b).Step 3.
Ift≥T, go to Step 4. Otherwise, set t=t+1 and go to Step 1.Step 4.
Check if constraints (20) are satisfied. If constraints (20) are satisfied, stop. Otherwise, solve the following linear programming problem to determine ω0 and ωt, t=1,…,T.(39)min∑t=1TλtU-Γω0+∑t=1Tωts.t. constraints (20)-(22).The obtained feasible solution can provide an upper bound for the optimal objective function value of (RSS3).(  4) Updating of the Lagrangian Multipliers. We initialize the Lagrangian multipliers by zero and update them according to the subgradient algorithm [36]. The iteration is stopped when the maximum number of iterations is reached or the relative duality gap (ZU-ZL)/ZL×100% is smaller than a certain threshold, where ZL and ZU are the best lower bound and upper bound obtained so far for (RSS3), respectively.


## 3.3. Discussions on Special Cases

(RSS) can be simplified in the following two special cases:Γ=∑t=1TλtL and Γ=∑t=1TλtU. In the case of Γ=∑t=1TλtL, we have λt=λtL for all t according to constraints (1a) and (1b) and (RSS) can be simplified as a deterministic multi-unit generation self-scheduling problem. Similar result can be obtained in the case of Γ=∑t=1TλtU. Therefore, we propose the following two-stage algorithm for the two special cases:Stage 1.
Simplify (RSS) by determiningλt,  t=1,…,T, according to the budget value Γ. If Γ=∑t=1TλtL, then λt=λtL,  t=1,…,T. If Γ=∑t=1TλtU, then λt=λtU,  t=1,…,T.Stage 2.
Solve the simplified problem by using the typical Lagrangian relaxation algorithm where constraints (10b) and (11) are relaxed.Since both the variable splitting-based Lagrangian relaxation algorithm and the two-stage algorithm proposed above can solve (RSS) in the two special cases, we will choose the more effective one for each special case based on the numerical test in Section4.


## 4. Numerical Results

In this section, we consider test cases of different sizes to implement the numerical experiments. For convenience, the variable splitting-based Lagrangian relaxation algorithm is denoted by A1 and the two-stage algorithm for the two special cases is denoted by A2. The organization of the section is as follows. First, the generation of the test cases is described in detail. Second, we discuss the impact of the weightγ on the performance of algorithm A1 and determine the value of γ in the experiments. Third, we test the performance of algorithms A1 and A2, respectively. Finally, we discuss the effect of the budget value Γ on the self-scheduling.We use Visual C++ to implement the proposed algorithms on a PC with 2.83 GHz and 3.25-GB memory. The linear programming involved is solved by calling CPLEX 12.5.


### 4.1. Test Cases

Parameters for the test cases are presented as follows. The number of units is set to 10, 50, and 100, respectively. The scheduling horizon is set to 24, 96, and 168 hours, respectively. The combination of the two parameters forms nine problem sizes, in which10×24, 10×96, and 50×24 are corresponding to small-sized problems, 10×168, 50×96, and 100×24 are corresponding to medium-sized problems, and 50×168, 100×96, and 100×168 are corresponding to large-sized problems according to the number of decision variables and the number of constraints included in the problem. For each problem size, ten test cases are generated randomly and tested. Therefore, a total of ninety cases are tested in the experiments.For each test case, we letSDi=RDi+PiL, SUi=RUi+PiL, and RDi=RUi for convenience. Value ranges for parameters associated with units are partially based on those in [27] and shown in Table 2. We use Pennsylvania–New Jersey–Maryland (PJM) Interconnection Real Time data from 2005 to 2006 to forecast price ranges and set λtL and λtU to the endpoints of the following confidence interval at the 95% confidence level, respectively:(40)λ¯t-t0.975m-1mst,λ¯t+t0.975m-1mstIn (40), λ¯t is the sample mean, st is the sample standard deviation, m is the sample size, and t0.975(m-1) is the 97.5% quantile of Student’st-distribution with m-1 degrees of freedom. The forecasted price data for 168 hours are provided in Figure 1. The budget lower bound Γ is set to ∑t=1TλtL+budget∑t=1TλtU-λtL where budget is allowed to vary within the set {0,0.2,0.4,0.6,0.8,1}. Note that budget=0 and budget=1 are corresponding to the two special cases discussed in Section 3.3. The convergence threshold of the relative duality gap is set to 1%.Table 2
Range of values for parameters associated with units and demands.
Parameter Range of values P i L (MW) [ 40,400 ] P i U (≤1200) (MW) [ 3 P i L, 4PiL] R D i / R U i (MW) [ 0.3 P i U, 0.6PiU] U i / L i (h) P i U < 600 [ 1,3 ] P i U ≥ 600 [ 2,5 ] y i 0 (h) On (0.4 chance) [ 1,10 ] Off (0.6 chance) [ - 10 , - 1 ] x i 0 if ui0=1 (MW) [ P i L, PiU] D t (MW) [ 0.7 ∑ i = 1 n P i U, ∑i=1nPiU] a 0 i ($) [ 0.8 P i U, 1.2PiU] a 1 i ($/MW) P i U < 600 [ 16 . 0,17 . 5 ] P i U ≥ 600 [ 18 . 0,19 . 5 ] a 2 i ($/(MW)2) P i U < 400 [ 0 . 03,0 . 06 ] P i U ≥ 400 [ 0 . 01,0 . 03 ] b 1 i ($) P i U < 600 [ 4 P i U, 6PiU] P i U ≥ 600 [ 2 P i U, 4PiU] b 2 i ($) [ 0.4 b 1 i, 0.6b1i] τ i (h) P i U < 600 [ 1,4 ] P i U ≥ 600 [ 3,6 ]Figure 1
Price data for 168 hours ($/MWh).


### 4.2. Computational Results

(  1) Discussion on Weight γ. We allow γ to vary within the set {0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1} and discuss the impact of the weight γ on the performance of algorithm A1 as follows. First, we carry out the algorithm for each test case under various budget and γ settings. Then, corresponding to each γ setting, we calculate the average relative duality gap and computation time over all test cases and report the results in Table 3. From Table 3, we can make the following observations:(1)
The average relative duality gap varies between 0.90% and 0.99%. The results show that algorithm A1 can meet the convergence requirement under allγ settings.(2)
The average computation time varies between 0.61 s and 0.99 s and reaches the minimum whenγ is set to 0.4. The results indicate that algorithm A1 shows the fastest convergence behavior with γ setting to 0.4.Table 3
Average relative duality gap and computation time under variousγ settings.
γ 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Relative duality gap (%) 0.93 0.93 0.94 0.93 0.92 0.91 0.90 0.94 0.96 0.99 0.98 Computation time (s) 0.99 0.96 0.87 0.71 0.61 0.63 0.75 0.86 0.90 0.96 0.98Based on the above discussion, we setγ to 0.4 for algorithm A1 in the following experiments.(  2) Performance of the Algorithms. We use the relative duality gap, the computation time, and the number of cases, denoted by N, in which the algorithm stops before reaching the maximum number of iterations to measure the performance of the proposed algorithms. The numerical results of algorithm A1 for all problem sizes and budget settings are reported in Tables 4 and 5. The numerical results of algorithm A2 for the two special cases are reported in Table 6. The reported relative duality gap and computation time are both the average of ten test cases in the same problem size.Table 4
Relative duality gap of algorithm A1.
Size Relative duality gap (%) budget Avg. 0 0.2 0.4 0.6 0.8 1 10 × 24 0.92 0.94 0.92 0.96 0.92 0.88 0.92 10 × 96 0.92 0.95 0.94 0.92 0.93 0.93 0.93 10 × 168 0.95 0.96 0.96 0.91 0.92 0.93 0.94 50 × 24 0.92 0.86 0.94 0.91 0.88 0.90 0.90 50 × 96 0.92 0.92 0.94 0.95 0.92 0.90 0.92 50 × 168 0.93 0.90 0.95 0.96 0.93 0.90 0.93 100 × 24 0.92 0.86 0.92 0.92 0.89 0.86 0.90 100 × 96 0.93 0.86 0.91 0.92 0.84 0.91 0.90 100 × 168 0.95 0.87 0.91 0.94 0.89 0.91 0.91 Avg. 0.93 0.90 0.93 0.93 0.90 0.90 0.92Table 5
Numerical results of algorithm A1.
Size N Computation time (s) budget Avg. budget Avg. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 10 × 24 10 9 10 10 10 10 10 0.30 0.26 0.27 0.24 0.17 0.01 0.21 10 × 96 10 10 10 10 10 10 10 0.37 0.26 0.25 0.22 0.17 0.04 0.22 10 × 168 10 10 10 10 10 10 10 0.52 0.35 0.31 0.30 0.25 0.11 0.31 50 × 24 10 10 10 10 10 10 10 0.52 0.22 0.28 0.27 0.20 0.03 0.25 50 × 96 10 10 10 10 10 10 10 1.26 0.51 0.45 0.39 0.34 0.20 0.53 50 × 168 10 10 10 10 10 10 10 2.68 1.07 0.83 0.72 0.69 0.52 1.09 100 × 24 10 10 10 10 10 10 10 0.55 0.27 0.24 0.27 0.23 0.06 0.27 100 × 96 10 10 10 10 10 10 10 1.51 0.84 0.69 0.66 0.60 0.39 0.78 100 × 168 10 10 10 10 10 10 10 3.74 2.55 1.39 1.24 1.17 1.05 1.86 Avg. 10 10 10 10 10 10 10 1.27 0.70 0.52 0.48 0.43 0.27 0.61Table 6
Numerical results of algorithm A2 for the two special cases.
Size Relative duality gap (%) N Computation time (s) budget = 0 budget = 1 budget = 0 budget = 1 budget = 0 budget = 1 10 × 24 0.65 0.96 10 10 0.02 0.01 10 × 96 0.81 0.98 10 10 0.10 0.08 10 × 168 0.96 0.98 10 10 0.29 0.18 50 × 24 0.76 0.96 10 10 0.04 0.03 50 × 96 0.99 0.94 10 10 0.69 0.16 50 × 168 0.99 0.95 10 10 1.58 0.43 100 × 24 0.86 0.96 10 10 0.14 0.06 100 × 96 0.99 0.93 10 10 1.37 0.29 100 × 168 0.98 0.93 10 10 3.07 0.75 Avg. 0.89 0.95 10 10 0.81 0.22From Tables4 and 5, we have the following observations:(1)
The average relative duality gap is 0.92% and the maximum relative duality gap is 0.96%. The computation time increases linearly with the increase of the problem size. The average computation time is 0.61 s and the maximum computation time is 3.74 s. The results demonstrate that algorithm A1 can find a solution very close to the optimal one in a reasonable time for the cases in all sizes.(2)
The average number of cases in which algorithm A1 stops before reaching the iteration limit is 10. The result shows that algorithm A1 can meet the convergence requirement within the iteration limit for almost all cases and indicates good convergence behavior of the algorithm.(3)
The computation time shows a decreasing trend whenbudget increases from 0 to 1. This is because the conservatism of the solution decreases as the budget lower bound increases.From Table6, we can observe that algorithm A2 can meet the convergence requirement before reaching the iteration limit for all test cases. For budget=0, the average relative duality gap is 0.89% and the average computation time is 0.81 s. For budget=1, the average relative duality gap is 0.95% and the average computation time is 0.22 s. A comparison between the results of algorithms A1 and A2 indicates that algorithm A2 converges faster than algorithm A1 for the special cases. This is because (RSS) is simplified as a deterministic generation self-scheduling problem in algorithm A2.Based on the numerical results in Tables4–6, we suggest using algorithm A2 to solve (RSS) in the special cases and algorithm A1 to solve the problem in other cases.As a comparison, we also solve (RSS2) by using the MILP approach where the quadratic fuel function is approximated by a ten-piece piecewise linear function, the exponential startup cost function is linearized, and the resulting MILP model is solved by calling CPLEX MILP solver. For each test case, the time limit is set to 1800 s. The numerical results are reported in Table7. In Table 7, columns 2-7 report the numbers of solvable cases in which the MILP approach can reach the optimality within the time limit, denoted by NS, under various problem sizes and budget settings and columns 9-14 report the average computation times over the solvable cases. For the unsolved cases, the MILP approach either cannot be optimally solved within the time limit or runs out of memory. If all cases cannot be optimally solved, the average computation time is not given.Table 7
Number of solvable casesa and average computation time of the MILP approach.
Size N S Computation time (s) budget Avg. budget Avg. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 10 × 24 10 10 10 10 10 10 10 1.65 0.39 0.37 0.37 0.37 0.36 0.59 10 × 96 10 10 10 10 10 10 10 70.21 8.79 8.52 8.45 8.52 8.42 18.82 10 × 168 0 10 10 10 10 10 8 - 42.04 41.64 41.46 41.36 44.69 42.24 50 × 24 10 10 10 10 10 10 10 551.79 4.52 4.36 4.14 3.83 4.01 95.44 50 × 96 0 0 0 0 0 0 0 - - - - - - - 50 × 168 0 0 0 0 0 0 0 - - - - - - - 100 × 24 0 10 10 10 10 10 8 - 19.23 18.63 18.84 18.49 17.23 - 100 × 96 0 0 0 0 0 0 0 - - - - - - - 100 × 168 0 0 0 0 0 0 0 - - - - - - - Avg. 3 6 6 6 6 6 5 - - - - - - - a Corresponding to cases that can be optimally solved by using the MILP approach within the time limit.From Table7, we can have the following observations:(1)
For the problem in small size, all cases can be optimally solved. For the problem in medium size, only partial cases can be optimally solved. For the problem in large size, none of the cases can be optimally solved. The largest size of the cases that the MILP approach can solve within the acceptable time is100×24.(2)
For the solvable cases, the computation time of the MILP approach grows exponentially as the problem size increases.The observations show the computation limit of the MILP approach and imply the importance of proposing an effective solution algorithm. Because the proposed algorithm can find the near-optimal solutions within the acceptable time for the cases in all sizes, it is more effective than the MILP approach, especially for the problem in medium or large size.(  3) Effect of the Budget Value on the Self-Scheduling. To show the effect of the budget value on the self-scheduling, we compare the generation levels and the objective function values under various budget settings, respectively. The numerical results are reported in Tables 8 and 9. From Tables 8 and 9, we can observe that given the problem size, both the generation level and the generation profit increase as budget increases. This is because the problem becomes less conservative with the increase of budget.Table 8
Average generation level under variousbudget settings.
Size Generation level (MW) budget Avg. 0 0.2 0.4 0.6 0.8 1 10 × 24 81858 115014 121810 124072 126819 127648 116203 10 × 96 398015 511554 536327 548467 559971 565509 519974 10 × 168 703614 885812 912438 934035 951012 958648 890927 50 × 24 422552 593972 616302 635507 644835 652059 594205 50 × 96 2046938 2592646 2659243 2713494 2763548 2787158 2593838 50 × 168 3570388 4489274 4600006 4681392 4759769 4800147 4483496 100 × 24 858882 1209880 1248902 1283975 1310584 1324268 1206082 100 × 96 3945357 4989315 5127241 5223253 5321735 5341850 4991459 100 × 168 7137883 8959433 9198064 9372351 9531833 9606864 8967738 Avg. 2129498 2705211 2780037 2835172 2885567 2907128 2707102Table 9
Average objective function value under variousbudget settings.
Size Generation profit ($) budget Avg. 0 0.2 0.4 0.6 0.8 1 10 × 24 891716 1520921 2481063 3529387 4638719 5854005 3152635 10 × 96 5746889 7767879 11021736 14558938 18378060 22535395 13334816 10 × 168 10738997 13711115 18570527 23872184 29601911 35815388 22051687 50 × 24 4569506 7826961 12692856 17991655 23687784 29867163 16105988 50 × 96 29829539 40525389 56480048 73811914 92586379 112876879 67685024 50 × 168 54385184 69585348 93813824 120348245 149010436 180186711 111221625 100 × 24 9213887 15795275 25797137 36504533 47984444 60668555 32660638 100 × 96 57959284 78591223 109414984 142950970 179071718 218221036 131034869 100 × 168 108862841 139454061 188030093 241240640 298751036 360937693 222879394 Avg. 31355316 41642019 57589141 74978718 93745610 114106980 68902964


## 4.1. Test Cases

Parameters for the test cases are presented as follows. The number of units is set to 10, 50, and 100, respectively. The scheduling horizon is set to 24, 96, and 168 hours, respectively. The combination of the two parameters forms nine problem sizes, in which10×24, 10×96, and 50×24 are corresponding to small-sized problems, 10×168, 50×96, and 100×24 are corresponding to medium-sized problems, and 50×168, 100×96, and 100×168 are corresponding to large-sized problems according to the number of decision variables and the number of constraints included in the problem. For each problem size, ten test cases are generated randomly and tested. Therefore, a total of ninety cases are tested in the experiments.For each test case, we letSDi=RDi+PiL, SUi=RUi+PiL, and RDi=RUi for convenience. Value ranges for parameters associated with units are partially based on those in [27] and shown in Table 2. We use Pennsylvania–New Jersey–Maryland (PJM) Interconnection Real Time data from 2005 to 2006 to forecast price ranges and set λtL and λtU to the endpoints of the following confidence interval at the 95% confidence level, respectively:(40)λ¯t-t0.975m-1mst,λ¯t+t0.975m-1mstIn (40), λ¯t is the sample mean, st is the sample standard deviation, m is the sample size, and t0.975(m-1) is the 97.5% quantile of Student’st-distribution with m-1 degrees of freedom. The forecasted price data for 168 hours are provided in Figure 1. The budget lower bound Γ is set to ∑t=1TλtL+budget∑t=1TλtU-λtL where budget is allowed to vary within the set {0,0.2,0.4,0.6,0.8,1}. Note that budget=0 and budget=1 are corresponding to the two special cases discussed in Section 3.3. The convergence threshold of the relative duality gap is set to 1%.Table 2
Range of values for parameters associated with units and demands.
Parameter Range of values P i L (MW) [ 40,400 ] P i U (≤1200) (MW) [ 3 P i L, 4PiL] R D i / R U i (MW) [ 0.3 P i U, 0.6PiU] U i / L i (h) P i U < 600 [ 1,3 ] P i U ≥ 600 [ 2,5 ] y i 0 (h) On (0.4 chance) [ 1,10 ] Off (0.6 chance) [ - 10 , - 1 ] x i 0 if ui0=1 (MW) [ P i L, PiU] D t (MW) [ 0.7 ∑ i = 1 n P i U, ∑i=1nPiU] a 0 i ($) [ 0.8 P i U, 1.2PiU] a 1 i ($/MW) P i U < 600 [ 16 . 0,17 . 5 ] P i U ≥ 600 [ 18 . 0,19 . 5 ] a 2 i ($/(MW)2) P i U < 400 [ 0 . 03,0 . 06 ] P i U ≥ 400 [ 0 . 01,0 . 03 ] b 1 i ($) P i U < 600 [ 4 P i U, 6PiU] P i U ≥ 600 [ 2 P i U, 4PiU] b 2 i ($) [ 0.4 b 1 i, 0.6b1i] τ i (h) P i U < 600 [ 1,4 ] P i U ≥ 600 [ 3,6 ]Figure 1
Price data for 168 hours ($/MWh).


## 4.2. Computational Results

(  1) Discussion on Weight γ. We allow γ to vary within the set {0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1} and discuss the impact of the weight γ on the performance of algorithm A1 as follows. First, we carry out the algorithm for each test case under various budget and γ settings. Then, corresponding to each γ setting, we calculate the average relative duality gap and computation time over all test cases and report the results in Table 3. From Table 3, we can make the following observations:(1)
The average relative duality gap varies between 0.90% and 0.99%. The results show that algorithm A1 can meet the convergence requirement under allγ settings.(2)
The average computation time varies between 0.61 s and 0.99 s and reaches the minimum whenγ is set to 0.4. The results indicate that algorithm A1 shows the fastest convergence behavior with γ setting to 0.4.Table 3
Average relative duality gap and computation time under variousγ settings.
γ 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Relative duality gap (%) 0.93 0.93 0.94 0.93 0.92 0.91 0.90 0.94 0.96 0.99 0.98 Computation time (s) 0.99 0.96 0.87 0.71 0.61 0.63 0.75 0.86 0.90 0.96 0.98Based on the above discussion, we setγ to 0.4 for algorithm A1 in the following experiments.(  2) Performance of the Algorithms. We use the relative duality gap, the computation time, and the number of cases, denoted by N, in which the algorithm stops before reaching the maximum number of iterations to measure the performance of the proposed algorithms. The numerical results of algorithm A1 for all problem sizes and budget settings are reported in Tables 4 and 5. The numerical results of algorithm A2 for the two special cases are reported in Table 6. The reported relative duality gap and computation time are both the average of ten test cases in the same problem size.Table 4
Relative duality gap of algorithm A1.
Size Relative duality gap (%) budget Avg. 0 0.2 0.4 0.6 0.8 1 10 × 24 0.92 0.94 0.92 0.96 0.92 0.88 0.92 10 × 96 0.92 0.95 0.94 0.92 0.93 0.93 0.93 10 × 168 0.95 0.96 0.96 0.91 0.92 0.93 0.94 50 × 24 0.92 0.86 0.94 0.91 0.88 0.90 0.90 50 × 96 0.92 0.92 0.94 0.95 0.92 0.90 0.92 50 × 168 0.93 0.90 0.95 0.96 0.93 0.90 0.93 100 × 24 0.92 0.86 0.92 0.92 0.89 0.86 0.90 100 × 96 0.93 0.86 0.91 0.92 0.84 0.91 0.90 100 × 168 0.95 0.87 0.91 0.94 0.89 0.91 0.91 Avg. 0.93 0.90 0.93 0.93 0.90 0.90 0.92Table 5
Numerical results of algorithm A1.
Size N Computation time (s) budget Avg. budget Avg. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 10 × 24 10 9 10 10 10 10 10 0.30 0.26 0.27 0.24 0.17 0.01 0.21 10 × 96 10 10 10 10 10 10 10 0.37 0.26 0.25 0.22 0.17 0.04 0.22 10 × 168 10 10 10 10 10 10 10 0.52 0.35 0.31 0.30 0.25 0.11 0.31 50 × 24 10 10 10 10 10 10 10 0.52 0.22 0.28 0.27 0.20 0.03 0.25 50 × 96 10 10 10 10 10 10 10 1.26 0.51 0.45 0.39 0.34 0.20 0.53 50 × 168 10 10 10 10 10 10 10 2.68 1.07 0.83 0.72 0.69 0.52 1.09 100 × 24 10 10 10 10 10 10 10 0.55 0.27 0.24 0.27 0.23 0.06 0.27 100 × 96 10 10 10 10 10 10 10 1.51 0.84 0.69 0.66 0.60 0.39 0.78 100 × 168 10 10 10 10 10 10 10 3.74 2.55 1.39 1.24 1.17 1.05 1.86 Avg. 10 10 10 10 10 10 10 1.27 0.70 0.52 0.48 0.43 0.27 0.61Table 6
Numerical results of algorithm A2 for the two special cases.
Size Relative duality gap (%) N Computation time (s) budget = 0 budget = 1 budget = 0 budget = 1 budget = 0 budget = 1 10 × 24 0.65 0.96 10 10 0.02 0.01 10 × 96 0.81 0.98 10 10 0.10 0.08 10 × 168 0.96 0.98 10 10 0.29 0.18 50 × 24 0.76 0.96 10 10 0.04 0.03 50 × 96 0.99 0.94 10 10 0.69 0.16 50 × 168 0.99 0.95 10 10 1.58 0.43 100 × 24 0.86 0.96 10 10 0.14 0.06 100 × 96 0.99 0.93 10 10 1.37 0.29 100 × 168 0.98 0.93 10 10 3.07 0.75 Avg. 0.89 0.95 10 10 0.81 0.22From Tables4 and 5, we have the following observations:(1)
The average relative duality gap is 0.92% and the maximum relative duality gap is 0.96%. The computation time increases linearly with the increase of the problem size. The average computation time is 0.61 s and the maximum computation time is 3.74 s. The results demonstrate that algorithm A1 can find a solution very close to the optimal one in a reasonable time for the cases in all sizes.(2)
The average number of cases in which algorithm A1 stops before reaching the iteration limit is 10. The result shows that algorithm A1 can meet the convergence requirement within the iteration limit for almost all cases and indicates good convergence behavior of the algorithm.(3)
The computation time shows a decreasing trend whenbudget increases from 0 to 1. This is because the conservatism of the solution decreases as the budget lower bound increases.From Table6, we can observe that algorithm A2 can meet the convergence requirement before reaching the iteration limit for all test cases. For budget=0, the average relative duality gap is 0.89% and the average computation time is 0.81 s. For budget=1, the average relative duality gap is 0.95% and the average computation time is 0.22 s. A comparison between the results of algorithms A1 and A2 indicates that algorithm A2 converges faster than algorithm A1 for the special cases. This is because (RSS) is simplified as a deterministic generation self-scheduling problem in algorithm A2.Based on the numerical results in Tables4–6, we suggest using algorithm A2 to solve (RSS) in the special cases and algorithm A1 to solve the problem in other cases.As a comparison, we also solve (RSS2) by using the MILP approach where the quadratic fuel function is approximated by a ten-piece piecewise linear function, the exponential startup cost function is linearized, and the resulting MILP model is solved by calling CPLEX MILP solver. For each test case, the time limit is set to 1800 s. The numerical results are reported in Table7. In Table 7, columns 2-7 report the numbers of solvable cases in which the MILP approach can reach the optimality within the time limit, denoted by NS, under various problem sizes and budget settings and columns 9-14 report the average computation times over the solvable cases. For the unsolved cases, the MILP approach either cannot be optimally solved within the time limit or runs out of memory. If all cases cannot be optimally solved, the average computation time is not given.Table 7
Number of solvable casesa and average computation time of the MILP approach.
Size N S Computation time (s) budget Avg. budget Avg. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 10 × 24 10 10 10 10 10 10 10 1.65 0.39 0.37 0.37 0.37 0.36 0.59 10 × 96 10 10 10 10 10 10 10 70.21 8.79 8.52 8.45 8.52 8.42 18.82 10 × 168 0 10 10 10 10 10 8 - 42.04 41.64 41.46 41.36 44.69 42.24 50 × 24 10 10 10 10 10 10 10 551.79 4.52 4.36 4.14 3.83 4.01 95.44 50 × 96 0 0 0 0 0 0 0 - - - - - - - 50 × 168 0 0 0 0 0 0 0 - - - - - - - 100 × 24 0 10 10 10 10 10 8 - 19.23 18.63 18.84 18.49 17.23 - 100 × 96 0 0 0 0 0 0 0 - - - - - - - 100 × 168 0 0 0 0 0 0 0 - - - - - - - Avg. 3 6 6 6 6 6 5 - - - - - - - a Corresponding to cases that can be optimally solved by using the MILP approach within the time limit.From Table7, we can have the following observations:(1)
For the problem in small size, all cases can be optimally solved. For the problem in medium size, only partial cases can be optimally solved. For the problem in large size, none of the cases can be optimally solved. The largest size of the cases that the MILP approach can solve within the acceptable time is100×24.(2)
For the solvable cases, the computation time of the MILP approach grows exponentially as the problem size increases.The observations show the computation limit of the MILP approach and imply the importance of proposing an effective solution algorithm. Because the proposed algorithm can find the near-optimal solutions within the acceptable time for the cases in all sizes, it is more effective than the MILP approach, especially for the problem in medium or large size.(  3) Effect of the Budget Value on the Self-Scheduling. To show the effect of the budget value on the self-scheduling, we compare the generation levels and the objective function values under various budget settings, respectively. The numerical results are reported in Tables 8 and 9. From Tables 8 and 9, we can observe that given the problem size, both the generation level and the generation profit increase as budget increases. This is because the problem becomes less conservative with the increase of budget.Table 8
Average generation level under variousbudget settings.
Size Generation level (MW) budget Avg. 0 0.2 0.4 0.6 0.8 1 10 × 24 81858 115014 121810 124072 126819 127648 116203 10 × 96 398015 511554 536327 548467 559971 565509 519974 10 × 168 703614 885812 912438 934035 951012 958648 890927 50 × 24 422552 593972 616302 635507 644835 652059 594205 50 × 96 2046938 2592646 2659243 2713494 2763548 2787158 2593838 50 × 168 3570388 4489274 4600006 4681392 4759769 4800147 4483496 100 × 24 858882 1209880 1248902 1283975 1310584 1324268 1206082 100 × 96 3945357 4989315 5127241 5223253 5321735 5341850 4991459 100 × 168 7137883 8959433 9198064 9372351 9531833 9606864 8967738 Avg. 2129498 2705211 2780037 2835172 2885567 2907128 2707102Table 9
Average objective function value under variousbudget settings.
Size Generation profit ($) budget Avg. 0 0.2 0.4 0.6 0.8 1 10 × 24 891716 1520921 2481063 3529387 4638719 5854005 3152635 10 × 96 5746889 7767879 11021736 14558938 18378060 22535395 13334816 10 × 168 10738997 13711115 18570527 23872184 29601911 35815388 22051687 50 × 24 4569506 7826961 12692856 17991655 23687784 29867163 16105988 50 × 96 29829539 40525389 56480048 73811914 92586379 112876879 67685024 50 × 168 54385184 69585348 93813824 120348245 149010436 180186711 111221625 100 × 24 9213887 15795275 25797137 36504533 47984444 60668555 32660638 100 × 96 57959284 78591223 109414984 142950970 179071718 218221036 131034869 100 × 168 108862841 139454061 188030093 241240640 298751036 360937693 222879394 Avg. 31355316 41642019 57589141 74978718 93745610 114106980 68902964


## 5. Conclusions

In this paper, we propose an improved Lagrangian relaxation algorithm for the robust generation self-scheduling problem under electricity price uncertainty in the deregulated electricity market. The problem includes quadratic fuel cost, exponential startup cost, unit operation constraints, and electric demand constraints. To avoid obtaining an unbounded relaxed problem, variable splitting is introduced into the algorithm. For the special cases of the budget value, we also propose a two-stage algorithm. Numerical results demonstrate the good performance of the proposed algorithms. Future research can be focused on generalizing the proposed algorithm to solve other robust generation self-scheduling problems.



---

*Source: 6303596-2018-07-22.xml*