Solving Unit Commitment Problem Using Modified Subgradient Method Combined with Simulated Annealing Algorithm

Abstract

This paper presents the solving unit commitment (UC) problem using Modified Subgradient Method (MSG) method combined with Simulated Annealing (SA) algorithm. UC problem is one of the important power system engineering hard-solving problems. The Lagrangian relaxation (LR) based methods are commonly used to solve the UC problem. The main disadvantage of this group of methods is the difference between the dual and the primal solution which gives some significant problems on the quality of the feasible solution. In this paper, MSG method which does not require any convexity and differentiability assumptions is used for solving the UC problem. MSG method depending on the initial value reaches zero duality gap. SA algorithm is used in order to assign the appropriate initial value for MSG method. The major advantage of the proposed approach is that it guarantees the zero duality gap independently from the size of the problem. In order to show the advantages of this proposed approach, the four-unit Tuncbilek thermal plant and ten-unit thermal plant which is usually used in literature are chosen as test systems. Penalty function (PF) method is also used to compare with our proposed method in terms of total cost and UC schedule.

1. Introduction

UC is very important problem for power system engineering. The problem can be described as a nonlinear, mixed-integer, and nonconvex and is considered to be a nondeterministic polynomial-time hard (NP-hard) problem [1]. The real difficulty in solving the problem is the high dimension of the possible solution space. Meta-heuristic and mathematical based methods have been developed for solving the thermal and hydrothermal UC problem in the literature. The most used metaheuristic methods are SA [2–5], expert systems [6], tabu search [7–9], evolutionary programming [10, 11], genetic algorithms [12–15], memetic algorithm [16], particle swarm optimization [17], interior point method [18], neural network [19, 20], and greedy randomized adaptive search procedure [21]. The mathematical-based methods depend on the duality theory [22–28]. The other methods used for solving the UC problem are dynamic programming [29, 30], mixed-integer programming [31], benders decomposition [32], and some hybrid methods [33–35].

Meta-heuristic methods are used for solving difficult combinatorial optimization problems. To solve UC problem by using these methods, priority list is precalculated and all the constraints are heuristically incorporated.

The LR method subsequently tries to find the values of the Lagrange multipliers that maximize the dual objective function based on the duality theory. The duality gap, a major problem in the nonlinear programming, has been long recognized as an inherent disadvantage of these methods. If the LR-based methods are used for solving the UC problem, dual solution may be far away from the optimal solution. The duality gap for the problem of UC is an important measure of the quality of the solution. When the gap is smaller, the solution is better [36]. In [37], different mathematical-based methods LR, penalty function, and augmented Lagrangian penalty function are compared to each other according to feasible cost, dual cost, duality gap, number of iterations, and duration time. According to [37], it is seen that there are differences between primal value and dual value. MSG method does not require any convexity and differentiability assumptions. In the nonlinear NP, the duality gap has been investigated and the theoretical tools for zero duality gap condition have been improved extensively in [38–42].

In this paper, one of the methods based on dual optimization technique, MSG method, which has the best performance in eliminating the duality gap in the literature, is used for solving the UC problem. A dual problem with respect to the sharp augmented Lagrangian is constructed for UC problem. The disadvantage of the MSG method is that the zero-duality gap value depends on the initial value of the upper dual value. This dual value is found by using SA algorithm. In this proposal approach constructs the dual problem and solves it without any duality gap for large class of nonconvex constrained problems. This proposed approach is compared to PF method because the MSG method removes some of the problems occurred in this method.

The remaining sections are outlined as follows. Section 2 provides a description of the UC problem formulation. This section includes an objective function and the constraints of the problem. In Section 3, MSG method with SA algorithm Explained in detailed. In Section 4, applications and numerical results are presented and discussed. Finally the conclusion is given.

2. Unit Commitment Problem Formulation

UC has been used to plan over a given time horizon the most economical schedule of committing and dispatching generating units to meet forecasted demand levels and spinning reserve requirements while all generating unit constraints are satisfied. The objective function can be represented mathematically as𝐹(𝑃^𝑡_𝑖,𝑈^𝑡_𝑖)=^𝑇𝑁_𝑡=1𝑖=1[𝐹_𝑖(𝑃^𝑡_𝑖)+SC^𝑡_𝑖(1−𝑈_𝑖^𝑡−1)]𝑈_𝑖^𝑡+^𝑇𝑁_𝑡=1𝑖=1SD^𝑡_𝑖×(1−𝑈_𝑖^𝑡−1)𝑈_𝑖^𝑡−1. (2.1) In (2.1), 𝑇 is period, 𝑁 is the number of generators, 𝑃^𝑡_𝑖 is the generation power output of unit 𝑖 at hour 𝑡, 𝐹_𝑖(𝑃^𝑡_𝑖) is quadratic fuel cost function of generating unit 𝑖 (𝐹_𝑖(𝑃^𝑡_𝑖)=𝑎_𝑖+𝑏_𝑖𝑃^𝑡_𝑖+𝑐_𝑖(𝑃^𝑡_𝑖)²), 𝑈^𝑡_𝑖 is status of unit 𝑖 at hour 𝑡 (on = 1, off = 0), 𝑆𝐶^𝑡_𝑖 is startup cost of unit 𝑖 at hour 𝑡, 𝑆𝐷^𝑡_𝑖 is shutdown cost of unit 𝑖 at hour 𝑡 [43].

The minimization of the objective function is provided to the following constraints.

System Constraints
Power Balance Constraints. For satisfying the load balance in each stage, the forecasted load demand should be equal to the total power of the generated power for feasible combination ^𝑁_𝑖=1𝑈^𝑡_𝑖𝑃^𝑡_𝑖=𝑃^𝑡_𝑑. (2.2) In (2.2), 𝑃^𝑡_𝑑 is system load demand at hour 𝑡.

Unit Constraints
Generation Limits. Each unit must satisfy the generation range 𝑃_𝑖,min𝑈^𝑡_𝑖≤𝑃^𝑡_𝑖≤𝑃_𝑖,max𝑈^𝑡_𝑖,𝑖=1,2,…,𝑁. (2.3) In (2.3), 𝑃_𝑖,min is minimum power output of unit 𝑖, and 𝑃_𝑖,max is maximum power output of unit i.
Ramp Up and Ramp Down Constraints. For each unit, output is limited by ramp up/down rate at each hour as follows: 𝑃_𝑖^𝑡−1−𝑃^𝑡_𝑖≤RD_𝑖,if𝑈^𝑡_𝑖=1&𝑈_𝑖^𝑡−1𝑃=1,_𝑡𝑖−𝑃_𝑖^𝑡−1≤RU_𝑖,if𝑈^𝑡_𝑖=1&𝑈_𝑖^𝑡−1=1, (2.4)RD_𝑖 is ramp down rate limit of unit 𝑖, and RU_𝑖 is ramp up rate limit of unit 𝑖.

The problem is nonconvex because it is structured binary variables. These variables cause a great deal of trouble and difficulty in solving the UC. Load balance is coupling constraint for the UC problem. The coupling constraints across the unit so that one unit affects what will happen on other units if the coupling constraints are met.

3. Modified Subgradient Method Combined with Simulated Annealing Algorithm

3.1. Simulated Annealing Algorithm for Unit Commitment

SA algorithm is a strong technique for solving hard combinatorial optimization problems without specific structure. This method has the ability of escaping local minima by performing uphill moves. The main advantage of SA algorithm is that it does not need large computer memory. SA is based on the iterative method, originally proposed by Metropolis et al. [44], which simulates the transition of atoms in equilibrium at a given temperature. The main disadvantage of this method is very greedy regarding computation time requirements, due to the large number of iterations needed for the convergence of the method. All the generating and acceptance distributions depend on the temperature [3].

One of the measure steps in the SA is the cooling schedule. This schedule has three components: the initial temperature setting, the temperature decreasing scheme, and a finite number of initial value of temperature should be large enough to allow all transients to be accepted for unconstrained optimization problem. On the other hand, the temperature is lowered based on multiplying the temperature in which typical values lie between 0.8 and 0.99. If the acceptance ratio is low or the sampled mean and variance of cost values at current temperature have big drops, then the factor is adjusted to a higher value in order to avoid getting stuck at a local optimal configuration. Otherwise, the factor is adjusted to a lower value to increase the convergence speed. The number of moves at each temperature is based on the requirement that at each temperature quasiequilibrium is to be restored. The final temperature is obtained if, at five consecutive temperatures, either the sampled mean values of cost function do not change or the acceptance ratio is small enough [2, 3].

In solving the UC problem by using SA algorithm, two types of variables need to be determined: the unit status (binary) variables U and the units output power (continuous) variables P. Then, this problem can be considered into two subproblems, a combinatorial optimization problem in U and a nonlinear optimization problem in P. The flow diagram of SA algorithm for UC problem is given in Figure 1 [3]. The economic dispatch problem is solved by using lambda-iteration method in the SA.

Figure 1 

SA algorithm for UC problem.

3.2. Modified Subgradient Method for Unit Commitment

The MSG method proposed by Gasimov has a notable performance in having zero duality gaps for large class of nonconvex problems [38, 39]. In the standard form the nonlinear programming problem can be expressed as follows: min_𝐾𝑓(𝐾) (3.1)subjecttoℎ(𝐾)≤0,𝐾∈Ω, (3.2)

where ℎ(𝐾)=[ℎ₁(𝐾)ℎ₂(𝐾)⋯ℎ_𝑛(𝐾)] is the constraint vector.

The primal problem is (3.1). Sharp augmented Lagrangian is defined as follows:𝐿(𝐾,𝑣,𝑐)=𝑓(𝐾)+𝑐‖ℎ(𝐾)‖−𝑣^𝑇ℎ(𝐾), (3.3) where 𝑣∈𝑅³ and 𝑐∈𝑅₊. Define the dual function as 𝐻(𝑣,𝑐)=min_𝐾∈Ω𝐿(𝐾,𝑣,𝑐). (3.4)

The dual problem (P *) ismax_{(𝑣,𝑐)∈𝑅3×𝑅+}𝐻(𝑣,𝑐). (3.5) Dual function H(v,c) is not a definite function; however it is a convex function for continuous f(K), h(K) and compact set S. Inequality constraints can be converted to standard form by adding nonnegative slack variables to left hand sides of inequalities; one of the advantages of the MSG method is that the fact that the solution convergence is proved. Thus, in each iteration a better solution can be found. This method differs from the classical LR functions, scans the solution with conics, and hence obtains zero dual gap value. In other words, there is no difference between primal and dual problems. Thus it gives an optimal solution. Another advantage of this method is that it does not require any convexity and differentiability assumptions.

Using the definitions, the MSG method is as follows.

Initialization
Choose a pair (𝑣₁,𝑐₁) with 𝑣₁∈𝑅³and𝑐∈𝑅₊, 𝑐₁≥0and let𝑗=1, and go to Step 1.

Step 1. Given (𝑐_𝑗,𝑣_𝑗), solve the following subproblem: min_𝐾𝑓(𝐾)+𝑐_𝑗‖ℎ(𝐾)‖−𝑣_𝑗ℎ(𝐾)=𝐻(𝑣_𝑗,𝑐_𝑗) (3.6) subject to 𝑘∈Ω.

In (3.6) (𝑐_𝑗,𝑣_𝑗) is multiplier of sharp augmented Lagrangian, 𝑓(𝐾) is objective function, ‖ℎ(𝐾)‖ is norm of the constraint vector, and 𝐻(𝑣_𝑗,𝑐_𝑗) is dual function of the problem. Let 𝐾_𝑗 be a solution of (3.6). If ℎ(𝐾_𝑗)=0, then stop; (𝑣_𝑗,𝑐_𝑗) is an optimal solution to the dual problem and 𝐾_𝑗 is a solution to (3.1); so 𝑓(𝐾_𝑗) is the optimal value of problem (3.1). Otherwise, go to Step 2.

Step 2. Update (𝑣_𝑗,𝑐_𝑗) by 𝑣_𝑗+1=𝑣_𝑗−𝑧_𝑗ℎ(𝐾_𝑗𝑐),_𝑗+1=𝑐_𝑗+(𝑧_𝑗+𝜀_𝑗)‖ℎ(𝐾_𝑗)‖, (3.7) where 𝑧_𝑗 and 𝜀_𝑗 are positive scalar step sizes defined below. Replace 𝑗 by 𝑗+1 and go to Step 1.

Step Size Calculation
Let us consider the pair (𝑣_𝑗,𝑐_𝑗) and calculate 𝐻(𝑣_𝑗,𝑐_𝑗)=min_𝐾∈Ω𝐿(𝐾,𝑣,𝑐), (3.8) and let 𝐻(𝐾_𝑗)≠0 for the corresponding 𝐾_𝑗, which means that 𝐾_𝑗 is not optimal. Then the step size parameter 𝑧_𝑗 can be calculated as 0<𝑧_𝑗<2(𝐻−𝐻(𝑣_𝑗,𝑐_𝑗))5‖ℎ(𝐾_𝑗)‖²,0<𝜀_𝑗<𝑧_𝑗, (3.9) where 𝐻 is an upper bound for the dual function. Considering the dual function formed by using the sharp Lagrangian, its value at any feasible point is not larger than primal problems’ objective function value. The equality occurs at a point when both primal and dual problems achieve their optimal values. It has been proven in [38] that if ℎ(𝐾)=0 for any K obtained from (3.6), then it is the solution of the primal problem. If ℎ(𝐾)≠0, then the value of 𝐻 calculated from (3.6) is strictly less than the optimal value of (3.1). In this case we update dual variables using Step 2, which leads to an increase in the value of the dual function. Solution of (3.6) corresponding to the updated (𝑣,𝑐) is always greater than the value obtained in the previous step. Note that this property is not guaranteed by the multiplier and penalty methods [43, 45]. According to MSG method for solving the UC problem, the objective function and constraints are defined as follows.

Objective Function
One has min𝐹(𝑃^𝑡_𝑖,𝑈^𝑡_𝑖)=^𝑇𝑁_𝑡=1𝑖=1[𝐹_𝑖(𝑃^𝑡_𝑖)+SC^𝑡_𝑖(1−𝑈_𝑖^𝑡−1)]𝑈^𝑡_𝑖+^𝑇𝑁_𝑡=1𝑖=1SD^𝑡_𝑖×(1−𝑈_𝑖^𝑡)𝑈_𝑖^𝑡−1. (3.10) Slack variables are added to inequality constraints. Equality and inequality constraints can be defined as follows.

Generation Limit Constraints
One has ℎ₁(𝑁,𝑇)=𝑈^𝑡_𝑖𝑃_𝑖,min−𝑃^𝑡_𝑖+ℎ^𝑡_𝑖ℎ=0,₂(𝑁,𝑇)=𝑃^𝑡_𝑖−𝑈^𝑡_𝑖𝑃_𝑖,max+ℎ^𝑡_𝑖=0. (3.11)

Load Balance Constraints
One has ℎ₃(𝑇)=^𝑁_𝑖=1𝑈^𝑡_𝑖𝑃^𝑡_𝑖−𝑃^𝑡_𝑑=0. (3.12)

Ramp Up and Ramp Down Constraints
One has ℎ₄(𝑁,𝑇)=𝑃_𝑖^𝑡−1−𝑃^𝑡_𝑖−RD_𝑖,+ℎ^𝑡_𝑖,if𝑈^𝑡_𝑖=1&𝑈_𝑖^𝑡−1ℎ=1,₅(𝑁,𝑇)=𝑃^𝑡_𝑖−𝑃_𝑖^𝑡−1−𝑅𝑈^𝑡_𝑖+ℎ^𝑡_𝑖,if𝑈^𝑡_𝑖=1&𝑈_𝑖^𝑡−1=1. (3.13)ℎ^𝑡_𝑖 shows slack variables from (3.11) to (3.13).
H = [ℎ₁(𝑁,𝑇)ℎ₂(𝑁,𝑇)ℎ₃(𝑇)ℎ₄(𝑁,𝑇)ℎ₅(𝑁,𝑇)]. According to (3.3), sharp augmented Lagrangian 𝐿_SAL is defined as follows:

𝐿_SAL=^𝑇𝑁_𝑡=1𝑖=1[𝐹_𝑖(𝑃^𝑡_𝑖)+SC^𝑡_𝑖(1−𝑈_𝑖^𝑡−1)]𝑈^𝑡_𝑖+^𝑇𝑁_𝑡=1𝑖=1SD^𝑡_𝑖×(1−𝑈^t_𝑖)𝑈_𝑖^𝑡−1+𝑐⋅[sqrt((𝑈^𝑡_𝑖𝑃_𝑖,min−𝑃^𝑡_𝑖+ℎ^𝑡_𝑖)²+(𝑃^𝑡_𝑖−𝑈^𝑡_𝑖𝑃_𝑖,max+ℎ^𝑡_𝑖)²+(^𝑁_𝑖=1𝑈^𝑡_𝑖⋅𝑃^𝑡_𝑖−𝑃^𝑡_𝑑)²+(𝑃_𝑖^𝑡−1−𝑃^𝑡_𝑖−𝑅𝐷^𝑡_𝑖+ℎ^𝑡_𝑖)²+(𝑃^𝑡_𝑖−𝑃_𝑖^𝑡−1−RU^𝑡_𝑖+ℎ^𝑡_𝑖)²)]−𝑣^𝑡_𝑖[(𝑈^𝑡_𝑖𝑃_𝑖,min−𝑃^𝑡_𝑖+ℎ^𝑡_𝑖)+(𝑃^𝑡_𝑖−𝑈^𝑡_𝑖𝑃_𝑖,max+ℎ^𝑡_𝑖)+(^𝑁_𝑖=1𝑈^𝑡_𝑖⋅𝑃^𝑡_𝑖−𝑃^𝑡_𝑑)+(𝑃_𝑖^𝑡−1−𝑃^𝑡_𝑖−RD_𝑖+ℎ^𝑡_𝑖)+(𝑃_𝑖^𝑡−𝑃_𝑖^𝑡−1−𝑅𝑈_𝑖^𝑡+ℎ^𝑡_𝑖)]. (3.14)

The epsilon value is chosen as 0.95z_j in the UC problem solution. After (3.14) is constructed using the objective function and constraints, all the steps of the MSG method are applied to the UC problem. In this study, SA algorithm is used in order to assign the appropriate initial value 𝐻 for MSG method. Then, the UC problem is solved using this proposed method for finding the optimal solution.

4. Applications and Numerical Results

The UC problem for four-unit Tuncbilek thermal plant is solved by using the MSG method combined with SA algorithm. First the SA algorithm is applied to the UC problem and then the cost function value found from the SA algorithm is used as upper dual value in MSG method for four-unit Tuncbilek thermal plant and ten-unit thermal plant. SA algorithm is coded in MATLAB PF and MSG methods are coded in GAMS which is a high-level modeling system for mathematical programming problems [46]. The data for four-unit Tuncbilek thermal plant is taken from Turkish Electric Power Company and Electricity Generation Company. The unit characteristics for four-unit Tuncbilek thermal plant are given in Table 1. In this study, a 24-hour day is subdivided into 8 discrete stages for four-unit Tuncbilek thermal plant. The load demands for the stages are given in Table 2. The unit characteristics for ten-unit thermal plant are given in Table 3. The load demands for each period are given in Table 4 for ten unit thermal plant. UC schedule for the SA algorithm for four-unit Tuncbilek thermal plant is given in Table 5. UC schedule for the SA algorithm for ten-unit thermal plant is given in Table 6.

Table 1 

Unit characteristics for four-unit Tuncbilek thermal plant.

Table 2 

Load data for four-unit Tuncbilek thermal plant (MW).

Table 3 

Unit characteristics for ten-unit thermal plant.

Table 4 

Load data for ten-unit thermal plant (MW).

Table 5 

UC schedule for SA algorithm for four-unit Tuncbilek thermal plant.

Table 6 

UC schedule for SA algorithm for ten-unit thermal plant.

The cost value $56249.9 found from the SA algorithm is used for upper dual value in MSG method for four-unit Tuncbilek thermal plant. The cost value $565907.9 found from the SA algorithm for ten-unit thermal plant. Optimal solution of the dual model set up by the augmented Lagrange function, which is proposed by Azimov and Gasimov for nonconvex problems, is equal to the primal solution of the system [38–40]. CONOPT2 is used as a GAMS solver for MSG and PF methods. The MSG method is run until the norm is equal to 9.713717×10⁻⁷ for four-unit Tuncbilek thermal plant; the norm is equal to 6.61489×10⁻⁴ for ten-unit thermal plant. Then all the constraints reduce to zero and feasible solution is obtained. Primal value is equal to the dual value yielding zero duality gap value for both of these systems. The value for the 𝑐 parameter is 9957.327 for four-unit Tuncbilek thermal plant; 𝑐 parameter is 13468.451 for ten-unit thermal plant. The UC schedule for MSG method combined with SA algorithm is given for four-unit Tuncbilek thermal plant and ten-unit thermal plant in Tables 7 and 8, respectively.

Table 7 

UC schedule for MSG method combined with SA algorithm for four-unit Tuncbilek thermal plant.

Table 8 

UC schedule for MSG method combined with SA algorithm for ten-unit thermal plant.

In Table 7, the total cost value is found $55581.683 and in Table 8 total cost value is found $563986.5 for MSG method combined with SA algorithm. It can be seen from Tables 7 and 8 that MSG reaches to zero duality gap value (primal value = dual value) for UC problem.

To show the advantages of this method, PF method is used. UC schedule, primal-dual values for PF method are given in Tables 9 and 10 for four-unit Tuncbilek thermal plant and for ten-unit thermal plant, respectively.

Table 9 

UC schedule for PF method for four-unit Tuncbilek thermal plant.

Table 10 

UC schedule for PF method for ten-unit thermal plant.

It is seen that there are differences between primal value and dual value in the PF method. The quality of the solution is improved when the duality gap is decreased. In this paper, zero duality gap is achieved and feasible solution is attained by using total cost value of the SA algorithm as an upper dual value for the MSG method.

5. Conclusion

The quality of the solution of UC problem is improved when the duality gap is decreased. In this paper zero duality gap is achieved and feasible solution is attained by applying the novel proposed method, MSG method combined with SA algorithm, to solving the UC problem. The cost function found from the SA algorithm is used for upper dual value in MSG method. The most attractive feature of the proposed approach is that the duality gap value of MSG method over the scheduled time horizon is zero. Note that there is a general acceptance that whenever the system size gets smaller, the duality gap value gets bigger for UC problem. However, it is shown that this is not the case with the MSG method for a small size system. The results of the proposed method for solving the UC problem are very impressive, and the quality of feasible solution is significantly improved. This approach can be applied to UC problem for any size of systems to obtain the feasible schedule.

List of Symbols

𝑐:	Multiplier of sharp augmented Lagrangian
𝑓(𝐾):	Objective function
𝐹𝑖(𝑃𝑡𝑖):	Generator fuel cost function in a quadratic form, 𝐹𝑖(𝑃𝑡𝑖)=𝑎𝑖+𝑏𝑖𝑃𝑡𝑖+𝑐𝑖(𝑃𝑡𝑖)2 ($/h)
𝐻:	Upper bound of the dual function
𝐻(𝑣,𝑐):	Dual function
ℎ(𝑥):	Equality constraint
𝑔(𝑥):	Inequality constraint
LR:	Lagrangian Relaxation
𝑁:	Number of generating units
𝑃:	Primal problem
𝑃∗:	Dual problem
PF:	Penalty function method
𝑃𝑡𝑑:	Nominal demand at hour 𝑡 (MW)
𝑃𝑡𝑖:	Generation output of unit 𝑖 at hour 𝑡 (MW)
𝑃𝑡max,𝑖:	Maximum available capacity of unit 𝑖 at hour 𝑡 (MW)
𝑃𝑡min,𝑖:	Minimum available capacity of unit 𝑖 at hour 𝑡 (MW)
RU𝑖:	Ramp up rate of unit 𝑖 (MW/h)
RD𝑖:	Ramp down rate of unit 𝑖 (MW/h)
SC𝑖:	Start up cost of unit 𝑖 ($)
SD𝑖:	Shut down cost of unit 𝑖 ($)
UC:	Unit commitment
𝑈𝑡𝑖:	Status value of unit 𝑖 at time 𝑡
𝑇:	Time horizon for UC (h)
TC:	Total cost
𝑣:	Multiplier of sharp augmented Lagrangian.

Acknowledgment

The support and guidance on MSG method by Professor Rafail N. Gasimov is gratefully acknowledged.