Evolutionary techniques applied to the optimal short-term scheduling of the electrical energy production

Evolutionary techniques applied to the optimal short-term scheduling of the electrical energy production Alicia Troncoso , Jose ´ C. Riquelme , Jesu ´s S. Aguilar-Ruiz , Jesu´s M. Riquelme Santos Abstract This paper presents an evolutionary technique applied to the optimal short-term scheduling (24 h) of the electric energy production. The equations that define the problem lead to a non-convex non-linear programming problem with a high number of continuous and discrete variables. Consequently, the resolution of the problem based on combinatorial methods is rather hard. The required heuristics, introduced to assure the feasibility of the constraints, are analyzed, along with a brief description of the proposed genetic algorithm (GA). The GA is used to compute the optimal on/off status of thermal units and the fitness function is obtained by solving a quadratic programming problem by means of a standard non-linear Interior Point (IP) method. The results from real-world cases based on the Spanish power system are reported, which show the good performance of the proposed algorithm, taking into account the complexity and dimensionality of the problem. Finally, an IP algorithm is adapted to deal with discrete variables that appear in this problem and the obtained results are compared with that of the proposed GA. Keywords: Genetic algorithms; Scheduling; Optimization; Feasibility; Interior point algorithms 1. Introduction The optimal short-term scheduling of the electrical energy production [1] aims at determining which generating units should be online and the corresponding optimal generation of thermal and hydro units along the scheduling period, usually 24 h, in order to minimize the expected total cost satisfying the forecasted system load. The scheduling task leads to a non-linear mixed-integer programming problem. Moreover, this problem is coupled in time by the maximum speed that generating units, specially thermal units, are able to change the produced energy (known as up and down ramps), and also by the topology of the hydroelectric power plants, with a delay in hours between the water of a reservoir being used and the availability of that water in the * Corresponding author. E-mail addresses: [email protected] (A. Troncoso), riquelme@ lsi.us.es (J.C. Riquelme), [email protected] (J.S. Aguilar-Ruiz), [email protected] (J.M. Riquelme Santos). reservoirs downstream. A really large number of variables, both continuous and discrete variables, is needed to properly model this problem. Many approaches have been proposed for the resolution of this optimization problem, ranging from Dynamic Programming to Linear Mixed-Integer Programming or Lagrangian Relaxation [2], the latter being the most widely used optimization method in commercial programs. Genetic Algorithms (GAs) [3,4], a general-purpose stochastic search method based on the mechanics of natural selection, have also been successfully applied to the electrical energy scheduling problem since the adaptation is quite straightforward due to the combinatorial nature of this problem. In the last few years, the Interior-Point (IP) or Logarithmic-Barrier class of methods has become the preferred numerical approach to solve non-linear optimization problems, as a consequence of its enhanced capability to deal with inequality constraints [5,6]. Since the early developments, intended for Linear Programming problems, many improvements have been proposed in order to extend the application of IP methods to convex and non-convex non-linear problems. These refinements have to do with stepsize control [7], line-search techniques to force convergence to a local minimum from arbitrary starting points [8,9], use of trust regions [10], among others. However, despite these improvements, its performance on mixed-integer problems is frequently disappointing because the solution is far from the global optimum [11]. This has raised the interest in computationally expensive algorithms taken from the artificial intelligence area such as evolutionary methods [12], tabu search [13], particle swarm optimization [14], simulated annealing [15] or ant colony optimization [16]. Although convergence to the global optimum cannot be theoretically guaranteed, the ability to escape from local minima makes them an attractive choice for many applications [17,18]. In this paper, a GA applied to the optimal shortterm (24 h) electrical energy production scheduling is presented. Some heuristics are included in order to assure the feasibility of the constraints that appear in the problem. Results from real-world cases based on Spanish power system are reported. The electric energy production scheduling presents a large number of variables and constraints, a non-convex non-linear objective function and integer and continuous variables. Thus, the main difficulty is to find feasible solutions and the novelty of the paper is addressed to improve the feasibility. The main contributions of the encoded GA can be stated as: a procedure to generate the initial population taking into account the ramp constraints; dynamic constraints of minimum limits on the hourly energy production of the thermal units to consider the starting and stopping periods as an alternative to the inclusion of other binary variables to model these states; a crossover operator adapted to the features of this problem leading to an adequate percentage of feasible individuals; and sparsity techniques and optimal ordering used to reduce the computational overhead. In spite of the sparsity techniques, a potential limitation of the GA is rather related to the CPU time when a complex topology of hydroelectric power plants is analyzed. The paper is organized as follows: Section 2pre- sents the equations used to model the scheduling problem, leading to a non-linear mixed-integer programming problem with a large number of both continuous and discrete variables. In Section 3a briefly description of the primal-dual IP algorithm is made. Section 4introduces the proposed GA, and several implementation issues that are crucial to obtain feasible solutions are discussed. Finally, Section 5reports some results obtained from realistic cases based on the Spanish power system, and the main conclusions of the paper are outlined. 2. Formulation of the problem The objective of the scheduling problem is to determine the on/off state and the energy production of thermal and hydro units at each hour of the scheduling period, in order to minimize the total cost of the system satisfying the forecasted hourly demand and the technical constraints of thermal and hydro power plants. The standard notation used for the scheduling of the electrical energy production problem is summarized in Table 1. This notation describes fixed parameters of the thermal and hydro units, indexes, number of elements and variables. 2.1. Objective function The total energy production cost of the scheduling period is defined by CT¼X nt t¼1X ng i¼1 ½CiðPi;tÞþSUiUi;tð1Ui;t1Þ þSDið1Ui;tÞUi;t1;ð1Þ where n t is the number of hours of the scheduling period, n g the number of thermal units, each having a quadratic cost function, C i (P i,t ), of the energy production, P i,t ;SU i and SD i are, respectively, the startup and shut-down cost of the thermal generator i, and U i,t is a binary variable representing the on/off state of the thermal generator iat hour t. It can be observed that the total production cost is a sum of quadratic functions of the energy of each thermal generator if the state of each generator was previously stated by the GA. This is the case of the proposed technique because the on/off states are managed by the GA. Notice that the production cost is only due to the production of thermal generators P i,t , i.e., generators that produce energy by burning a fuel or by atomic means. Hydro units provide free-of-charge energy PH h,t that is only subject to the availability of water in the corresponding reservoirs. 2.2. Constraints The minimization of the objective function is subject to technical constraints, water balance in hydroelectric power plants and the associated reservoirs, and to the system energy demand and reserve balances: •Maximum and minimum limits on the hourly energy production of the thermal and hydro generators, Pm i6Pi;t6PM i;i¼1;...;ng;t¼1;...;nt;ð2Þ PHm h6PHh;t6PHM h;h¼1;...;nh;t¼1;...;nt; ð3Þ Table 1 Definition of the data and variables of the problem Data or fixed parameters SU i The start-up cost of the thermal unit i(€) SD i The shut-down cost of the thermal unit i(€) C i (Æ) Quadratic cost function of the thermal unit i(€) Pm iLower bound of the hourly energy production of the thermal unit i(MWh) PM iUpper bound of the hourly energy production of the thermal unit i(MWh) PHm hLower bound of the hourly energy production of the hydro plant h(MWh) PHM hUpper bound of the hourly energy production of the hydro plant h(MWh) VHm hLower bound of the water level of the reservoir hin terms of energy (MWh) VHM hUpper bound of the water level of the reservoir hin terms of energy (MWh) UR i Upper bound of the up rate of the thermal unit i(MWh/h) DR i Lower bound of the down rate of the thermal unit i(MWh/h) W h Inflow of the reservoir hin terms of energy (MWh) D t Energy demand at hour t(MWh) R t Generating capacity in reserve at hour t(MWh) DT i Number of hours that the unit imust be shut-down after stopping UT i Number of hours that the unit imust be functioning after starting d(k) Water delay time between reservoir kand the next reservoir downstream (in h) n(k) Next reservoir downstream regarding the reservoir k Indexes and number of elements iThermal unit index hHydro plant index tHour index n g Number of thermal units n h Number of hydro plants n t Number of hours of the scheduling period Variables P i,t Energy production of the thermal unit iat hour t(MWh) U i,t On/off state of the thermal generator iat hour t PH h,t Energy production of the hydro plant hat hour t(MWh) VH h,t Stored energy of the reservoir hat hour t(MWh) where n h is the number of hydro plants, PH h,t the energy production of hydro plant hat hour t, and Pm i,PM i,PHm hand PHM hare the limits on the hourly energy production of the thermal unit i and hydro plant h, respectively. Eq. (2) cannot be fulfilled when thermal generators are either starting or stopping, as starting and stopping periods begin, respectively, when the corresponding state changes to ON or OFF. In order to avoid this problem, this equation is modified for thermal units that are either being started-up or shut-down, 06Pi;t6PM i;i¼1;...;ng;t¼1;...;nt:ð4Þ Moreover, the energy produced by thermal units during periods of shutting-down (U i,t = 0) is out of the optimal scheduling. Consequently, penalty terms proportional to this energy are added to the objective function as follows: C0 T¼CTþX nt t¼1X ng i¼1 CpPi;tð1Ui;tÞ:ð5Þ •Maximum up and down ramps of thermal units. The thermal units can not increase or decrease the production of energy at consecutive hours by more than a given maximum rate, DRi6Pi;tPi;t16URi;i¼1;...;ng; t¼1;...;nt;ð6Þ where UR i yDR i are, respectively, the maximum up and down rates of the thermal generator i, usually known as ramp limits. •Limits on the available water. The hydro units use water to generate electrical energy and water is a limited resource. Thus, the energy produced by a hydro unit is limited by the volume of available water in the associated reservoir. In consequence, reservoir levels are subject to capacity limits, VHm h6VHh;t6VHM h;h¼1;...;nh; t¼1;...;nt;ð7Þ where VH h,t is the stored energy of reservoir hat hour t, corresponding to the hydro unit h;VHm h and VHM hare respectively the minimum and maximum limits on the stored energy imposed by the maximum and minimum possible water level of reservoir h. •Hydraulic coupling between reservoirs. Time coupling exits due to cascaded reservoirs, since the water used to produce energy in a hydro unit will be available later to the next hydraulic unit downstream with a certain delay, obviously when the water has arrived to the corresponding reservoir. VHh;t¼VHh;t1PHh;tþX nðkÞ¼h PHk;tdðkÞþWh; ð8Þ where d(k) is the water delay time in hours between reservoir kand the next reservoir downstream, n(k), that is supposed to be reservoir h, and W h is the natural inflow of reservoir h. •The total hourly energy production must be equal the total energy demand at that hour, D t , which has been previously forecasted. X ng i¼1 Pi;tUi;tþX nh h¼1 PHh;t¼Dt;t¼1;...;nt: ð9Þ •The total energy that can be produced at each hour must exceed the forecasted demand by a specified amount, R t , i.e., the generating capacity in reserve to be used if an unexpected event such as the failure of a plant or a large error on the forecasted demand happens. X ng i¼1 PM iUi;tþX nh h¼1 PHM hPDtþRt; t¼1;...;nt:ð10Þ •Minimum up and down times of thermal units. The minimum up time, UT i , is the minimum number of hours that the unit imust be functioning after starting. Besides, the minimum down time, DT i , is the minimum number of hours that the unit imust be shut-down after stopping. X DT i1 k¼0 ð1Ui;tþkÞPDT i if unit iis shut-down at hour tð11Þ X UT i1 k¼0 Ui;tþkPUT i if unit iis started at hour t:ð12Þ Start-up and shut-down costs of realistic cases tend to reduce the number of shut-downs and start-ups to a minimum, making the minimum-time constraints useless in most cases. Moreover, the inclusion of hydraulic generation facilitates the fulfillment of the thermal unit constraints because the hydro units are faster in response and produce energy at no cost, i.e., the hydraulic energy will be strategically distributed among the hours of the scheduling horizon in order to avoid the starting of more thermal units than the strictly required. As an example, Table 2 shows the number of constraints, binary and continuous variables of the above problem for a test system comprising 49 thermal units, two hydro units and the scheduling horizon embracing 24 h. 3. Primal-dual IP algorithm Among the distinctive features of the above optimization problem, the most important are: large number of variables and constraints in practical cases; non-convexity of the objective function; and presence of integer and continuous variables. The IP methods only can be applied when all the variables of the problem are continuous. An alternative frequently used in practice consists in relaxing the discrete nature of U i,t and imposing instead the next constraint: 06Ui;t61:ð13Þ This leads to a simplified model without discrete variables which requires that a heuristic procedure be applied during or at the end of the iterative process in order to determine the best integer value for every U i,t . Thus a mixed-integer programming problem is considered from a continuous perspective. In addition, every discrete variable can be handled as a continuous variable provided that the following quadratic constraint is added: Ui;tð1Ui;tÞ¼0:ð14Þ However, this procedure increases the non-linearity and non-convexity of the initial problem. Based on the above comments and practical experience, the following procedure has been chosen to solve the optimal short-term scheduling of the electrical energy production. (1) The problem is solved adding the constraint (Eq. (13)). (2) Using the obtained solution in the before step to initialize the IP algorithm, solve the problem including the constraint (Eq. (14)). A brief description of a classical IP algorithm is provided next. This method introduces auxiliary positive slack variables in order to turn inequality restrictions into equality constraints: xj6xj,xjþsj¼xjsjP0;ð15Þ where x j represents any variable subject to a limit, and s j is the corresponding positive slack variable. In order to guarantee the positiveness of the slack variables, logarithmic penalty terms are included in the objective function by means of a penalty factor lthat is progressively reduced throughout the iterative process [5]. f0ðxj;sj;lÞ¼fðxjÞlX j ln sj:ð16Þ The main steps of the IP algorithm are the following: (1) Initialize the variables so that the slack variables are positive. (2) Initialize the penalty factor lso as to make the logarithmic terms dominate over the original objective function. (3) The minimization of the corresponding Lagrangian function is performed by solving the non-linear optimality equations using an one-step Newton’s algorithm, and the optimal increment of primal and dual variables is computed. (4) The step-length ais reduced, if necessary, so that the slack variables remain positive. Lagrange multipliers associated to equality constraints arising from the introduction of auxiliary slack variables must also remain positive because optimality conditions lead to equations of the form: sjzj¼l;ð17Þ where z j is the Lagrange multiplier associated with the slack variable s j . (5) Update primal and dual variables taking into account the necessary step-length limitation. (6) Reduce the penalty factor l. The proportional relationship between the penalty factor and the duality gap (dugap) defined by Eq. (17) Table 2 Dimension of the problem for a test system Number of constraints Number of variables Binary Continuous (2 Æn g +3Æn h +2)Æn t +2Æn g n g Æn t (n g +2Æn h )Æn t 2642 1176 1272 provides the most common approach to reduce this penalty factor: l¼cPnl j¼1sjzj nl ;ð18Þ where c61andn l is the number of inequality constraints of the original problem. Steps 3–6 are iteratively repeated until optimality conditions are satisfied and the penalty coefficient l, and consequently the average duality gap, is small enough. More sophisticated versions of steps 3 and 4, including line searches, modified Hessians, etc. [8–10] could be needed in the non-convex case, in order to avoid divergence or convergence to unacceptable points. However, such refinements have not been actually implemented because the behaviour of the IP method has proven good enough for the application tested. The application of Newton’s method to solve the non-linear optimality equations yields a very large, sparse linear system, specially when ramp and hydraulic couplings are considered. Consequently, sparsity techniques and optimal ordering [19] must be used to reduce the computational overhead. Fig. 1 shows the fill-ins generated when solving a small example (five hydro plants, five thermal plants and a 5-hour scheduling horizon), with and without optimal ordering. It can be noted that the fill-ins are reduced a 50% approximately when an optimal ordering is made. Table 3 shows relative execution times of the IP algorithm for two realistic problems (73 thermal plants, 24 h, 8 and 30 reservoirs, respectively). As can be noticed, execution times grow considerably with the number of variables, specially when standard ordering is performed. Finally, Fig. 2 shows the relative overhead of the different processes comprising the IP algorithm. 4. The proposed genetic algorithm As presented in the previous section, the optimal scheduling of the electric energy production is a non-linear, non-convex, combinatorial, mixed-inte- ger and very large problem. Hence, there is no technique that would always lead to the optimal solution of the problem for realistic cases. In the last years, techniques based on heuristics, dynamic programming, linear mixed-integer programming and lagrangian relaxation have been applied to this particular problem. Techniques based on heuristics rely on simple rules that depends on the knowledge of power plant operators. Constraints of realistic problems are not properly modelled by dynamic programming approaches, and the number of required states increases exponentially, thus leading to excessive computation times. Linear programming approaches cannot properly model neither the non-linear objective function nor the non-linear constraints, and crude approximations are required. Finally, the use of heuristic techniques is required by lagrangian relaxation approaches to calculate feasible solutions, deteriorating the quality of the obtained solutions. Consequently, new methods are still needed to obtain more optimal solutions to realistic problems. In this paper, a GA [20,21] has been used to solve the scheduling problem due to its ability to deal with non-linear functions and integer variables. The proposed GA algorithm is used to compute the optimal on/off states of thermal units, i.e., the binary variables, while the optimal continuous variables, i.e., the hourly energy production of hydro and committed thermal units, are calculated solving a typical quadratic programming problem by a classical IP optimization algorithm in which the on/off states of thermal units are known. Convergence characteristics of GA depend on several key implementation issues that are discussed in the rest of this section. 4.1. Codification of the individuals Each individual is represented by the on/off states of thermal generators during the scheduling period. Thus, individuals are represented by 0/1 matrices, with columns corresponding to time scheduling intervals and rows associated with thermal units. If the element (i,j) is equal to one, the state of thermal unit iduring time interval jis on. Similarly, if the element (i,j) is equal to zero, the state of thermal unit iduring time interval jis off. Fig. 3 shows the representation of a certain individual of the population. It can be observed that the thermal unit 1 is on from 1am to 4am and the rest of hours is off; the thermal unit 2 is off from 1am to 6am and the rest of hours is on; the thermal unit 3 is on during all scheduling horizon; the thermal unit 4 is only on from 11am to 4pm, etc. 4.2. Initial population Up and down ramp constraints of thermal units (Eq. (6)) are a key factor in the convergence of the GA: if the initial population is strictly randomly selected, ramp constraints lead to many infeasible individuals in the initial generation, which makes successive generations suffer from poor diversity, and the GA may converge prematurely. To assure that the initial population contains an adequate 0 20 40 60 80 100 0 20 40 60 80 100 Columns Rows 0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100 Columns Rows Fig. 1. Structure of the linear system including fill-ins. percentage of feasible individuals, initial on/off schedulings are randomly selected but modified to account for the minimum start-up and shut-down times imposed by ramp constraints. For example, if generator g, with a maximum down ramp equal to 100 MWh, is on at hour 3 producing an energy of 400 MWh, this generator would require 4 hours to shut-down and, consequently, the generator at hours 4, 5 and 6 should be on. The state U g,3 is strictly randomly generated but the states for the following hours, U g,4 ,U g,5 and U g,6 , are given by Ug;3¼1)Ug;4¼Ug;5¼Ug;6¼1:ð19Þ 4.3. Fitness function The fitness function evaluates the quality of an individual of the population. In this case, the function is the inverse of the total production cost of the individual. The total production cost is obtained solving a quadratic programming problem by using a non-linear Interior Point method [22,23].An extra-high-cost fictitious generator is included to satisfy the system demand (Eq. (9)). This fictitious generator generates the necessary energy that the rest of generators cannot produce to satisfy the demand of the customers. A penalty term proportional to the deficit in reserve requirements is added in the cost function aiming at satisfying the reserve constraint. Penalty terms only apply to infeasible individuals, which are consequently eliminated throughout the evolutionary process. 4.4. Selection operator To produce a new generation, parents are randomly selected using a roulette wheel selection technique that selects the best individuals for reproduction. The probability of a particular individual being selected is in proportion to its fitness function, taking into account that the total generation cost, including possible penalizations, is being minimized. The individuals chosen to be parents are included in the following generation. 4.5. Crossover operator Offspring is obtained by adding the binary strings that results from random partitions of each row, as shown in Fig. 4a. A column-partitioning procedure may also be applied (Fig. 4b). This crossover operator is a particular case of the multi-point crossover operator where the number of points is equal to the number of rows or columns, respectively. As rows are associated with the thermal units, the first approach yields mainly the infeasibility of new individuals in terms of minimum up and down times (Eqs. (11) and (12)), while the second approach has an effect bigger on the constraint of the demand (Eq. (9)) and the reserve (Eq. (10)). This above statement is shown in the next example. Let be a test Table 3 Relative execution time and iteration number. Number of variables Natural ordering Optimal ordering Iterations Time Iterations Time 2376 33 4.00 34 1.00 3960 27 140.98 34 2.45 Others 23% Optimal ordering 6% Linear system building 10% Linear system solution 61% Fig. 2. Relative overhead of the different processes. Hours of the Schedulin g Horizon Thermal Units 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 ..... 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 Fig. 3. Representation of an individual of the population. system comprising four thermal units and a 4-hour scheduling horizon. The minimum up and down time are 2 hours for all thermal units. Two individuals selected to be parents are represented in Fig. 5. Fig. 6 shows the children obtained using a crossover operator by rows from the random partition (1,2,1,2). It can be observed that the children do not satisfy the Eqs. (11) and (12) (generator three) while the Eqs. (9) and (10) can be feasible. Fig. 7 shows the children obtained using a crossover operator by columns from the same partition. It can be observed that the individual on the right does not satisfy the Eqs. (9) and (10) since all generators are off at hours 3 and 4. However, in this case the Eqs. (11) and (12) are fulfilled. The crossover probability has been set to one, i.e., two individuals that have been selected to be parents are always combined to obtain a new individual. In the final version of the GA, the crossover by rows has been chosen because start-up and shutdown costs of realistic cases, along with the inclusion of hydraulic generation, tend to reduce the number of shut-downs and start-ups to a minimum, making the minimum-time constraints useless in most cases. All the rows are always combined to obtain a new individual, though probabilities might have been used to determine which rows should be combined. 4.6. Mutation operator After the crossover process, the individuals of the population are mutated to introduce some new genetic material according to a pre-defined mutation probability p. Consequently, the percentage of mutated individuals of a generation is equal to 100p%. The mutation of an individual means the mutation of an only gene. The gene to be mutated 1 2 3 t 1 2 3 t 11 22 33 .. .. .. gg 1 2 3 t 1 2 3 t 11 22 33 .. .. .. gg PARENTS CHILDREN 123 t 123 t 11 22 ... ... gg 123 t 123 t 11 22 ... ... gg CHILDREN PARENTS ab Fig. 4. Crossover Operator: (a) random partitions of rows and (b) random partitions of columns. h1 h2 h3 h4 h1 h2 h3 h4 g1 1 1 1 1 g1 1 1 0 0 g2 1 1 0 0 g2 1 1 0 0 g3 0 0 0 0 g3 1 1 1 1 g4 1 1 0 0 g4 0 0 1 1 Fig. 5. Parents. h1 h2 h3 h4 h1 h2 h3 h4 g1 1 1 0 0 g1 1 1 1 1 g2 1 1 0 0 g2 1 1 0 0 g3 0 1 1 1 g3 1 0 0 0 g4 1 1 1 1 g4 0 0 0 0 Fig. 6. Children by partitions of rows. h1 h2 h3 h4 h1 h2 h3 h4 g1 1 1 1 1 g1 1 1 0 0 g2 1 1 0 0 g2 1 1 0 0 g3 1 1 1 1 g3 0 0 0 0 g4 0 0 1 1 g4 1 1 0 0 Fig. 7. Children by partitions of columns.