An adaptive random search for short term generation scheduling with network constraints

This paper presents an adaptive random search approach to address a short term generation scheduling with network constraints, which determines the startup and shutdown schedules of thermal units over a given planning horizon. In this model, we consider the transmission network through capacity limits and line losses. The mathematical model is stated in the form of a Mixed Integer Non Linear Problem with binary variables. The proposed heuristic is a population-based method that generates a set of new potential solutions via a random search strategy. The random search is based on the Markov Chain Monte Carlo method. The main key of the proposed method is that the noise level of the random search is adaptively controlled in order to exploring and exploiting the entire search space. In order to improve the solutions, we consider coupling a local search into random search process. Several test systems are presented to evaluate the performance of the proposed heuristic. We use a commercial optimizer to compare the quality of the solutions provided by the proposed method. The solution of the proposed algorithm showed a significant reduction in computational effort with respect to the full-scale outer approximation commercial solver. Numerical results show the potential and robustness of our approach.


Introduction
A definition of Short Term Generation Scheduling (Unit commitment) used in Smart grid operations is a scheduling problem with two objectives: The power dispatch that requires distributing the system load to the generating units over a given time horizon and the start-up and shutdown schedules of the power output generators.
The planning horizon in the short term usually lasts from 24 to 168 hours (see for instance [1]). Additionally, the Short Term Generation Scheduling solves a very large-scale, time-varying, non-convex, mixed-integer optimization problem. For these reasons, this problem belongs to the set of NP-hard problems.
The Short Term Generation Scheduling (STGS) with network constraints in an NP-hard Mixed Integer Non Linear Problem. We readily know that exacts methods are inefficient for large power systems (see [2][3][4][5]). Because of this, we present an alternative for generating good solutions in short computing time based on an Adaptive Random Search strategy.
One advantage when the Lagrangian relaxation is used is to eliminate constraints that complicate the structure of the original problem. Because in our work we use Lagrangian relaxation, the model becomes drastically simplified and this allows to use the AGS algorithm in a natural way.
Lagrangian relaxation (LR) is the most widely used method to solve STGS by relaxing some constraints [6,7]. This method is based on the duality theory, and tries to find optimal dual variables that maximize the dual Lagrangian.
A way to effectively implement LR is to apply "divide and conquer principle", that is, to decompose the STGS into a master problem and several manageable subproblems to be solved separately. The sub-problems are linked by Lagrange multipliers that are added to the master problem to yield a dual problem. The dimensionality of such dual problem makes it easier to solve than the primal. The multipliers are updated through different techniques. A popular choice is the use of any subgradient method, but the mayor difficulty associated with this option is the feasibility of the initial solution. This phenomenon is due to the dual nature of the algorithm.
The commitment states obtained in LR generally are infeasible because STGS is a MINLP (Mixed Integer Non-Linear Programming) problem. The duality gap is an inherent disadvantage of this technique, i.e. the dual solution may be far away from the optimal solution (see [8]). The feasible commitment states are obtained after adjustment. Usually heuristic methods are needed to modify the dual solution obtained for LR into a feasible solution. The ways to obtain these feasible solutions may vary notably.
For instance, in [17] GA is used to identify the model parameters for an Ultracapacitor, based on time-domain data. In [18] GA is employed to extract the optimal model parameters based on the Hybrid Pulse Power Characterization (HPPC). Finally in [19] a GA is proposed for effectively achieving optimal component sizing of a hybrid energy storage system in an electric vehicle.
In [15] a new approach via multi-particle swarm optimization (MPSO) to solve the unit commitment (UC) problem is presented meanwhile in [13,14] a simulated annealing algorithm (SAA) is presented to solve UC.
However the obtained results by EP, GA, SA, PSO and TS required a considerable amount of computational time especially for a large system size. This kind of techniques are not suitable for STGS problem due its NP-hard nature.
Because the main objective of an effective method for solving UC with LR is how to obtain feasible solutions. We present an alternative for generating feasible solutions in short computing time. Our proposal is based on an Adaptive Gibbs Sampling (AGS) algorithm.
The random search is one of the pillars of most heuristic methods to solve engineering optimization problems. The success of these methods to find good solutions (near-optimal) in an optimization problem is mainly achieved by tuning parameters. The perturbation of a parameter makes it possible to explore large regions of a landscape and potentially escape from local optima, resulting in the exploration of different local optima. The GA introduces a stochastic perturbation in their mutations [20], while SA presents these perturbations through their temperature levels and the cooling task [21]. For instance, a bad selection of the perturbation parameters in these methods could be result in a large risk of getting trapped in local regions. This fact points out of a need for an accurate selection of the step size parameters which dictate the amount of noise in the random search. On the other hand, the optimal scale of this perturbation, in order to achieve a good balance between exploration and exploitation, depends on the shape of the search landscape associated with the optimization problem. Therefore, this dependence makes the parameter selection a major issue in the design of heuristic algorithms.
In this study, we present an adaptive random search approach based on the Markov Chain Monte Carlo method to address a short term generation scheduling with network constraints.
The key to the proposed method is that the noise level of the random search can be adaptively controlled according to the landscape. Thus, the noise intensity allows exploring the entire search space, and noise-reduction allows exploitation in the promising regions where local optima exist. The algorithm is well described and compared in unconstrained global optimization problems by [22]. The effectiveness and robustness of this method was proved in several complex problems in order to find reasonable quality solutions for global optimization problems. Motivated by the success of the performance of the method over complex problems, we decided to apply our method in this engineering optimization problem. Besides, we use a fullscale outer approximation commercial solver to compare the quality of the solutions provided by the proposed method.
This paper is organized as follows: in the following section we present the problem formulation of STGS. In section "Solution Methodology", we introduce a novel optimization method which we call AGS algorithm based on infeasible solutions calculated for LR algorithm. We present our computational experience in the section this called, additionally we show our results on three test systems. Finally we draw our conclusions and summarize future work in the last section.

Problem formulation
The STGS with network constraints problem consists in determining the mix of generators and their estimated output level to meet the expected demand of electricity over a given time horizon (a day or a week), while satisfying the load demand, spinning reserve requirement and transmission network constraints.
In this work, we address a STGS based on the notation presented in [23], where network constraints are represented through a DC model (see [1]) and we consider a multi-period time horizon. The objective is to minimize a function that includes fixed costs, start-up costs and operating costs. A second order polynomial describes the variable costs as a function of the electric power. The following notation is used in the mathematical model: nr Reference node with angle zero.

Decision variables:
If plant j is started up at the beginning of period k; 0 otherwise ( Objetive function: Constraints: Spinning reserve X Generation limit Transmission capacity limits Start-up and shut-down of power units Angular limit voltage The objective Eq (1) is to minimize the start up cost A j y jk and the operating cost of each plant. The operating cost of each plant j includes a fixed cost F j v jk and a variable cost E j (t jk ). There is a power balance Constraint (2) per node and time period. In each period, the production has to satisfy the demand and losses in each node. Power line losses are modeled through cosine approximation and it is assumed that the demand for electric energy is known, and is discretized into t periods. There are many approximations to model power line losses, some of them are linear, and some are non-linear. Further details of the cosine approximation can be found in [24]. Spinning reserve requirements are modeled in Eq (3). In each period the running units have to be able to satisfy the demand and the pre-specified spinning reserve. In Eq (4), each unit has a technical lower and upper bounds for the power production. The transmission capacity limits of lines in Eq (5) serve the purpose of avoiding problems in the dynamic stability of the system. The Constraint (6) describes how the units start-up, run and shut-down (a running unit cannot be started-up). Finally, the angle in all buses has lower and upper bounds given by Eq (7).

Solution methodology
In this paper, we extend the AGS algorithm (see [22]) for continuous optimization to tackle mixed-variable optimization problems. First, the LR was used just for bounding purposes, and then the AGS algorithm was used to construct feasible solutions in reasonable computational times.

Lagrangian relaxation framework
Because STGS is an NP-Hard problem, the global optimal solution cannot be obtained for large scale power systems. For this reason, we use LR to calculate a lower bound of the optimal solution. The main disadvantage of this method is the difference between primal and the dual solutions, which is the duality gap. This situation determines that solutions obtained by LR are infeasible for the original problem.
LR is based on the duality theory, which tries to find optimal dual variables that maximize the dual Lagrangian. These dual variables (Lagrangian multipliers) need to be updated in order to improve the lower bound. In this work, we use subgradient method to update the lagrangian multipliers. The parameters of the subgradient and other especifications are shown in [25].
For LR we decomposed the STGS problem into n subproblems, one per generation node [25]. DICOPT solver (see [26]) was used to maximize the Lagrangian bound. The LR algorithm was used just for bounding purposes. Since the main objective of an effective method for solving STGS with LR is to obtain feasible solutions, we present an alternative to constructing feasible solutions in short computing time based on the AGS algorithm.
In this paper, the objective function used by the AGS algorithm does not consider an explicit mechanism for handling constraints. For this reason we use the LR in the original problem. All of the system constraints are dualized in original objective function. This manipulation is necessary due to fact that AGS needs a function without constraints.
Applying Lagrange duality to the Constraints (2), (3) and (5) in STGS yields: Dual Function: Dual Subproblem: where λ nk is the Lagrangian multiplier associated to a power balance constraint of node n in period k; μ k is the Lagrangian multiplier associated to a spinning reserve requirement in period k; γ nk , β nk are the Lagrangian multipliers associated to a transmission capacity limits of node n in period k. The above model is subject to Constraints (4), (6) and (7), called box constraints. Dualizing these constraints produces a dual subproblem that is less expensive to solve and speed up the solution of this subproblem. The algorithm initialized with a set of Lagrange multipliers, in this case, we defined these multipliers heuristically as a result of the knowledge of the original problem. Then, we use subgradient method to improve Lagrangian multipliers.
As the Lagrangian function constituted by the dualization of the complicating constraints is concave and non-differentiable, the AGS algorithm allows optimization of the dual function, since it is able to optimize such functions. The dual function can be decomposed into subproblems, one for each generation unit, however this procedure is not explored in this paper.
A brief summary of the subgradient method is given in Algorithm 1 below. Determine the subgradient direction d k of the function L Z DS at α k 7 Determine step size t k λ k (UB − L Z DS (α k ))/||d k || 2 8 Update multiplier vector by using α k+1 max{0, α k +t k Á d k } The adaptive gibbs sampling algorithm The AGS algorithm is based on the Markov Chain Monte Carlo (MCMC) method combined with the one-dimensional Metropolis-Hastings algorithm and the multi-dimensional Gibbs sampler. This MCMC algorithm is called Metropolis-within-Gibbs (MWG) algorithm, and was suggested in [27]. The proposed optimization heuristic is a population-based method that generates a set of new potential solutions through a random process given by the MWG algorithm. The global information about the landscape is extracted through the population of solutions generated at each iteration. With the global information generated in the step described above, the random process identifies the most promising regions of the search space and then uses this information to generate another set of new potential solutions.
The main key of the proposed method is that the noise level of the random search can be adaptively controlled according to the landscape. Thus, the noise intensity allows exploring the entire search space, and the noise reduction allows the exploitation in the regions where local optima exist. In order to improve the obtained solutions, we consider coupling a local search into a random search process. The algorithm is well described and implement for unconstrained global optimization problems in [22]. A brief summary of the AGS method is given in Algorithm 2 below. Selection(X ) 7c Update Intensification(x) 12 if (x better than S best ) then 13 S best x 14 until (termination criteria are not met); 15 return S best General form of AGS algorithm.
Step 0: Initialization: Randomly select an initial solutionx ¼ ðx 1 ; x 2 ; :::; x n ; :::; x N Þ within the feasible region and go to Step 1. Note that in this step we could also provide an initial solution obtained by other method. In this paper we use the solution that provides the subgradient method.
Step 1: Sampling: Generate candidate points for each variable by where Z is a standard normal random variable and c n is a scale parameter. The candidate point will be accepted as the next value x tþ1 n ¼ x Ã n with probability P ¼ min 1; pðx 1 ; x 2 ; :::; x Ã n ; :::; x N Þ pðx 1 ; x 2 ; :::; x t n ; :::; If the candidate point is not accepted, then the current value of x is retained: x tþ1 n ¼ x t n . Simulating one value in turn for each individual variable is called one cycle of Gibbs sampling where a newx vector solution is built. We can draw a population of M solutions by performing M Gibbs cycles. The output of the sampling step is a populationX and a vector of acceptance ratesỹ. Finally, go to Step 2.
Step 2: Selection: Estimate a mode solutionx for each variable in the populationX and go to Step 3.
Step 3: Update: Adjust the scale parametersc by the following rule where c o n is a constant chosen so that initially the acceptance rates are close to zero. The actual iteration number, τ, is initialized at the beginning with τ = 1 and will be increased iteratively by τ = τ + 1. Go to Step 4.
Step 4: Mutation: Replace the variable value by a random value within the search space according to the following rule: If y n > ε; then x n ¼ rand; c n ¼ c o n and t ¼ 1; where rand is an uniform random variable within the feasible region. If hỹi > b, where hỹi is an average of acceptance rates vector over all variables, go to Step 5; Otherwise return to Step 1.
Step 5: Intensification: Improve the solutionx via a local search strategy. We can use an arbitrary local search method. In this paper, we use the Nelder-Mead method as a local search strategy [28]. Finally, return to Step 1.

Parameter settings:
The AGS parameters used in this paper are chosen such as to be a robust setting and therefore, in our experience, applicable to a wide range of optimization problems. The parameters used are: population size M = 100, initial scale parameters c o = (0.1, . . ., 0.1), β = 0.7, ε = 0.95 and λ = 2.
The algorithm stops when the number of function evaluations is reached, or when the neighbor solution was not improved after the time period elapses. In Step 2, the global information generated by the population allows identifying the most promising (or more likely) regions of the search space. Therefore, the starting point for the next iteration will belong to such promising region. Note that in the Steps 3 and 4, the noise level of the random search can be adaptively controlled through acceptance rates. The acceptance rates provide information about the landscape. Values close to 1 on the acceptance rates allow to identify local optima and exploiting on them via local search. In addition, the mutation mechanism allows the escape and avoids being trapped in local optima, exploring other regions of the search space.
In this way, the noise intensity allows exploring the entire search space, and the noise reduction allows exploitation in the promising regions where local optima exist.
AGS algorithm for STGS problem. Before using the AGS algorithm to solve the STGS problem, we must define the representation of a solution. A solutionx is composed of two parts. The first part contains continuous variables, and the second contains integer variables, that is,x ¼x c [x b . In the STGS problem,x c ¼ ft jk ; d nk g contains the variables of power output of plant j in period k and the angle of node n in period k, respectively.x b ¼ fv jk ; y nk g contains the variables v jk that represents if plant j is committed in period k and y nk represents if plant j is started up at the beginning of period k.
In the initialization process, Step 3 in Algorithm 2, an initial solution is provided by the method RL (see Section "Lagrangian Relaxation Framework"). Next, in Step 4, the integer variables are fixed, and continuous variables are modified by the Gibbs cycles. The integer variables can be modified after Step 9 in Algorithm 2, by the following rule: After modifying the y nk value, v jk will make a copy of the value, that is, v nk = y nk . The Algorithm 3 shows the exchange of information between the Lagrangian relaxation scheme and AGS algorithm. Selection(X ) 7c Update(c o , λ) 8x Mutation(ỹ; ε) 9 Generate k at random, k = 1, 2, . . .K 10 if (y nk−1 = 1) then 11 y nk 0 12 v nk y nk 13 if (y nk−1 = 0) then 14 r rand(0, 1) 15 Intensification(x) 19 if (x better than S best ) then 20 S best x 21 until (termination criteria are not met); 22 return S best

Computational experience
The AGS algorithm was developed in C++ and simulated on a desktop computer with an AMD Phenom TM II N970 Quad-Core with a 2.2 GHz processor and 8 GB RAM. For Lagrangian Relaxation, the mathematical model of STGS was implemented in the modeling environment GAMS using the solver DICOPT (see [26] and [29]) for solving the MINLP problems (dual subproblems and full STGS scale model).
Three test systems are presented to evaluate the performance of the proposed AGS algorithm in the engineering optimization problems. The test cases used in the experimentation are power systems of 24, 104 and 118 units with a planning horizon of 24 h.
In order to illustrate the structure of test systems, we choose the basic information of the IEEE 24-bus test system. The single-line diagram of the IEEE 24-bus test system with 24 nodes, 24 thermal units and 38 transmission lines is shown in Fig 1. The data used in the IEEE-24 bus system were based on the original system which is available at [30]. The data of the minimum and maximum capacities generating units are presented in Table 1. Table 2 lists the costs, initial state and power output of each generating unit at time 0. In Table 3 and Fig 2 the load profile is illustrated. The node location of the loads, as well as load at each node as a percentage of the total system demand are presented in Table 4. The transmission lines data is given in Table 5. The lines are characterized by the nodes that are connected, as well as the reactance and the capacity of each line.
The 104-bus system data were extracted from [24] and correspond to the energy system of mainland Spain. The data used in the IEEE 118-bus system were extracted from the original  system which is available at [31]. This data contains information about reactance and capacitance of the transmission lines, the demand profile, and costs of the generation.

Results
The computational complexity of the problem is shown in Table 6. This complexity depends on the number of thermal units connected in each test system. This table shows the number of variables and constraints for each test system. As the AGS algorithm is a stochastic approach, 20 runs are executed for AGS on each test case, and the average costs of the 20 runs are determined. The comparison results of AGS algorithm and full-scale outer approximation commercial solver (DICOPT) are shown in Table 7. These results show that AGS solutions are very close to the GAMS solutions z Ã . The GAP in all cases is less than 0.05%. Fig 3 shows the average Gap for the three test systems using AGS. Standard deviation is also shown in Table 7. The GAP between the best solution through GAMS solver and the AGS is calculated by Additionally, it can be seen in Table 7, that the CPU time of proposed algorithm is much less than DICOPT. Fig 4 shows

Conclusions and future work
In this paper, we present a novel optimization method which we call Adaptive Gibbs Sampling (AGS) algorithm to address a Short Term Generation Scheduling with network constraints, which determines the startup and shutdown schedules of thermal units over a given planning  horizon. The proposed heuristic is a population-based method that generates a set of new potential solutions via random search strategy. The random search is based on a Markov Chain Monte Carlo method. The key to the proposed method is that the noise level of the random search is adaptively controlled in order to exploring and exploiting the entire search space. In order to improve the solutions, we consider coupling a local search into a random search process. This paper proposes an enhanced method that combines Lagrangian relaxation and AGS to solve STGS. The AGS algorithm is based on infeasible solutions calculated for LR algorithm.
We evaluated the performance of AGS algorithm against a full-scale outer approximation commercial solver in order to compare the quality of the solutions provided by the proposed method. Several groups of instances are tested to evaluate the performance of the proposed heuristic. The experimental results show that AGS algorithm is robust, as it is capable of finding reasonable quality solutions for this engineering optimization problem. Our AGS method converges to the near-optimal solution at a faster rate than the direct solution obtained by the solver DICOPT. In addition, AGS produces much tighter bounds on the optimal solution values than standard Lagrangian Relaxation.
In the future, we will extend the experimentation in other test systems with different structures in order to evaluate the performance of the proposed algorithm, and compare their results against other heuristic methods.