Application the differential evolution for solving the problem of minimizing the power shortage of electric power systems

Based on the experience of domestic and foreign researchers, it is known that there are various mathematical models, software systems, and optimization methods used to solve the set tasks for assessing the resource adequacy of electric power systems (EPS). However, the continuous development of EPS leads to the complication and integration of systems against the background of which it becomes necessary to take into account an increasing number of its elements and parameters in one task. Thus, for more effective (in terms of speed and accuracy) solutions of modified models, it is required to analyze and search the most appropriate set of optimization methods. In this connection, the purpose of this study is to analyze the applicability and the effectiveness of applying the method of differential evolution and gradient optimization methods for the model of minimizing the power deficit, which should be also compared. The article considers the analysis of the results of the work that optimization methods, studies were conducted on a test isolated EPS, with various tuning parameters. As a result, it was confirmed that these methods could be used to solve the problem. From the point of view of accuracy and resources expended on calculations, the most efficient among the implemented methods was the method of differential evolution, which was confirmed by numerical experiments on the small systems.


Introduction
Today's electricity consumers place high demands on ensuring the reliability of the electricity supply. This is due to the cost of a power failure with economic damage and situations dangerous to life and health of people that ensue. Interruption of electric power supply to consumers is due to failures of electric power system (EPS) equipment. In order to minimize the number of electricity constraints for consumers, it is necessary to implement a set of technical and organizational measures to improve the reliability of the EPS in advance. One of the main means of ensuring the reliability of the EPS is the early planning of the development of the system itself and the redundancy of its elements. Since maintaining the redundancy of generating capacities and the grid part of the EPS are costly measures, the justification of the redundancy of all types requires a qualified assessment. For this purpose, the resource adequacy of prospective EPS schemes is assessed. The result of the assessment is reliability indicators that lend themselves to an economic interpretation.
One of the stages of the resource adequacy assessment when applying the Monte Carlo method [1] is to determine the power shortages of possible states of the EPS. The basis for computing power shortages is a simulation of the EPS, which includes a mathematical model of the EPS, as well as optimization methods to obtain the power shortage amount for each of the considered states of the system. The quality of the results, including the speed and accuracy of the calculation, the ability to solve problems with an increasing number of optimized parameters, depends on the applied optimization method and the correctness of the mathematical model. The statement of the problem of minimizing the power shortage can be presented both in a linear and nonlinear form [2]. The most adequate is the statement in the nonlinear form, where losses in power lines have a quadratic dependence on the transmitted power [3].
In known practices at home and abroad, various optimization methods are used to solve this problem, so the "Amber" software and computer system (SCS) [4][5] makes use of the method of internal points, the "ORION-M" SCS [6] employs a dual simplex method. At present, in the USA several different software and computer systems are adopted, namely: GE MARS [7], GridView [8], MARELI (PROMOD IV) [9], SAM (Supply Adequacy Model) [10], N-Area Reliability Program (NARP) [11], as well as PLEXOS [12] are all commercially available, closed projects (models and methods are not disclosed). In European countries, they use the RTE Antares Simulator open-source SCS with several customizable linear mathematical models and optimization methods developed by RTE. The employed methods are characterized by the accuracy and speed of the power shortage calculation.
At present, in engineering applications, heuristic methods are increasingly used to solve optimization problems, and one of them is the method of differential evolution [13][14]. The method of differential evolution is used to find the global extremum of non-differentiable, non-linear, multi-modal (possibly having a large number of local extrema) functions of many variables. At present, there are 5 efficient modifications of the method available but besides that, a great increase in the method's performance can be achieved by adopting the technology of parallel and vector calculations.
As it is known, heuristic methods often have low performance in comparison to direct optimization methods, because they mainly use the stochastic approach and more evaluations of the objective function of each formed solution. Applicability and efficiency of the method often depend on the problem being solved; this article covers the method of differential evolution in its original form, and its comparative analysis together with the method of conjugate gradients with respect to the problem of minimization of the power shortage is carried out.

Problem statement
The problem of minimizing the power shortage is formulated as follows: to determine the optimal flow distribution in an EPS for known values of operable generating capacities, required levels of consumers' loads, transmission capacities of EPS connections and power loss coefficients in EPS connections [1], [4][5]. There are several types of models for minimizing the power shortage, and this paper will apply a model with nonlinear balance constraints, which takes into account the quadratic power losses. Mathematically, the problem is formulated as follows: when the balance constraints are respected: As well as constraints on optimized variables: 0 ≤ ≤ ̅ , = 1, … , , where: x_i -power used in zone i (MW), ( x) ̅ _iavailable power in zone i (MW), y_i -the load served in zone i (MW), y ̅ _i -the amount of load in zone i (MW), z_ij -power flow from zone i to zone j (MW), z ̅ _ijbandwidth of the power transmission line between nodes i and j (MW), a_ji -specified positive coefficients of specific power losses during its transfer from zone j to zone i, j≠i, i=1,…,n, j=1,…,n.

Minimizing power shortage
Since model (1-6) is a non-linear programming problem, one can use various methods of conditional and unconditional optimization to solve it. However, this problem cannot be solved by standard methods of unconditional optimization due to various constraints that are of the type of equations and inequalities, for this purpose it is necessary to transform the objective function and all constraints into those of the single objective function type or apply optimization methods, where these constraints act as parameters of the optimization method.
In the present studies, the penalty function method is applied to the transformation of a conditional optimization problem to an unconditional one. The penalty function method can be applied to optimization problems with various types of constraints. The method enables us to transform an initial problem with constraints into a problem, the solution of which can be obtained by methods of unconditional optimization. Such transformation allows not only to use various methods of unconditional optimization but also to increase the accuracy of calculations given the correct selection of setup parameters. The main changes are made to the objective function that has constraints added to it in the form of penalty functions. Thus, changes in the system can result in triggering the penalty function, the value of which will begin to increase dramatically. In this case, the response to the penalty will be regulated by the optimization method and, ultimately, the function will be directed to the desired solution.
The studies of the power shortage search that is efficient in terms of time and effort were conducted within the framework of the following sets of methods: a set of methods of penalty functions and the method of conjugate gradients with the Fletcher-Reeves coefficient, where the step length value is calculated as the Armijo rule, as well as a set of methods of penalty functions and differential evolution. The complete study included a larger set of methods that were also implemented programmatically, but due to the inefficiency of some of the methods, they passed only the first stage of the studies with their results presented in the part of this article that deals with experimental studies. The complete list of implemented methods included: the gradient descent method, steepest descent method (with and without step normalization), conjugate gradient method (the variants that come with Fletcher-Reeves and Polak-Ribière coefficients) and differential evolution method. The following algorithms have been implemented as one-dimensional optimization and line search methods to calculate the step length value in the steepest descent method and conjugate gradient method: 2 different algorithms of the Golden-section method, the combined Brent's method, the Powell's method, methods based on the Armijo rule conditions, the strong Wolfe conditions, the Armijo-Goldstein conditions, the parabolic method. Due to the instability of the obtained results, most of the above described onedimensional optimization methods were not used. The main method to be tested was the line search by the Armijo rule. The conjugate gradient method is an iterative numerical method (of the first order) for solving optimization problems, which allows us to determine the extremum (minimum or maximum) of the objective function [15][16][17]. The conjugate gradient method is a further development of the fastest descent method, which combines two concepts: the gradient of the objective function and the conjugate direction of vectors. In general, the process of finding the minimum of a function is an iterative procedure, the algorithm of which can be described by the following set of steps:

Algorithm of the employed conjugate gradient method
Step 1: Analytical expressions (in their symbolic form) are defined to calculate the gradient of function ∇f(x_1,x_2,…,x_n ) using formula (7) Step 2: The initial approximation is set = { 1 , … , } Then the iterative process is performed.
Step 3: The necessity to restart the algorithmic procedure for zeroing the last direction of search is determined. As a result of the restart, the search is carried out anew in the direction of the steepest descent.
Step 4: The coordinates of unit vector ( 1 , 2 , … , ) are calculated using the formula obtained in Step 1, and the coordinates of the new point are determined when moving in the direction of the unit vector as a function of the calculation step.
calculation of the weight coefficient and unit vector of conjugate directions at the current calculation step (the Fletcher-Reeves formula): -for the first step of the calculation ( = 0), the weighting coefficient is not calculated (the same applies to the case of the algorithm restart), and the unit vector of the conjugate directions is determined as follows: -for the following calculation steps (k=1,2,…), the weighting coefficient and unit vector of the conjugate directions are calculated on the basis of the following ratios: In this case, the coefficient calculated using the Fletcher-Reeves formula is presented as a formula below formula (9).
Step 5: we determine the calculation step length based on the condition of the extremum search for the following function = { ± • ( )} (the solution of the onedimensional optimization problem).
⟹ { ± • ( )} → Step 6: New values of the function arguments are defined after the k-th step of the calculation: where the "+" sign is used to find the maximum of a function and the "-" sign is used to find the minimum of a function; Step 7: Checking the stopping criteria of the iterative process. The calculational process ends when the point at which the gradient estimate is zero (response function coefficients become insignificant) is reached. Otherwise, there is a return to Step 3 and the iterative calculation continues.
To find the step length value, one has to solve the problem using one-dimensional optimization methods. However, in practice, a complete solution of the problem is either not achieved due to the complexity of the function or it takes a large amount of time and internal iterations to find a solution. A different approach can be used to reduce the number of operations: the values of the calculation step length are selected so that they meet the condition presented below.
The condition (the Armijo rule) is an adaptive method of searching for the value of the calculation step length, which indicates that function { ± • ( )} should not exceed the value of some decreasing linear function equal to ( ) at the zero point: where coefficient σ∈(0,1) and the calculation step length λ are determined iteratively by multiplying the initial step length λ_0 by coefficient β∈(0,1)until the condition is met.
The algorithm for determining the optimization problem calculation step length as per the Armijo rule can be represented by the following procedure: Step 1. Set coefficient σ within the range from 0 to 1 and the initial step length valueλ_0.
Search procedure (verifying that the the Armijo rule condition is respected) Step 2. If the Armijo rule condition is not met, then it is necessary to adjust the calculation step length λ_k= λ_0 •β^k, where variable β can take any value from 0 to 1. By default, variable β is assigned a value of 0.5, and k is the current iteration number of the search.
Step 3. If the Armicho rule condition is met, then the calculation step length can be assumed to be λ=λ_k and the search procedure is completed.
This rule requires a single calculation of the gradient, after which a small number of iterations are spent on selecting the appropriate step length. Each of these nested iterations, in turn, requires the value of the objective function to be calculated without a gradient, i.e. the tests performed are relatively lightweight. It should be noted that this condition is satisfied for all sufficiently small λ.
It should be noted that in the course of the studies the corrective values were selected independently and were determined as β=0.85, with coefficient σ= 10e-4,, under the conditions of experiments, the number of iterations spent on the search did not exceed 14, and in general, this rule had a high rate of convergence and also provided sufficient accuracy of the calculation step length for optimization methods.

The algorithm of the adopted method of differential evolution
Differential evolution is a method of multi-dimensional mathematical optimization that belongs to the class of stochastic optimization algorithms and uses some of the ideas of genetic algorithms but does not require working with variables in the binary code [18][19][20].
This method requires only the possibility to calculate the values of the objective functions, but not those of its derivatives, so it is a direct "method". Differential evolution is intended to find a global minimum (or maximum) of non-differentiable, non-linear, multi-modal (a large number of local extrema) functions of many variables. The method is easy to implement and use (it contains few control parameters that require their selection) and can be parallelized.
The algorithm of the method of differential evolution can be represented as follows: Step Step 3. Cross-breeding: a) set the mutation probability coefficient CR = [0, 1]; b) form a vector of random numbers P ∈[0, 1] with the number of dimensions equal to that of R1 c) form a "child vector" CH, the ordinal number in P is greater than CR, the gene from is inherited from R1 or, otherwise, MV. d) evaluate the objective function for CH vector values; Step 4. Selection: (a) Compare vectors R1 and CH; b) introduce a vector with a lower value of the objective function into the new population.
Go to Step 2 as part of the cycle.
Step 5. Check whether the limit on the number of generations has been reached.

Experimental studies of optimization methods to solve the problem of power shortage minimization
Within the scope of the studies and the software implementation of algorithms and mathematical models a personal computer with the suite of software products and technical specifications indicated in Table 1 was used.  [17] At the first stage of the studies, the correctness of the implementation of gradient methods and the method of differential evolution was tested. For this purpose, widely known special functions with known minimum values were used. First of all, the testing was carried out on the Rosenbrock function [21], a non-convex function used to evaluate the performance of optimization algorithms proposed by Howard Rosenbrock in 1960. The Rosenbrock function for two variables is defined as: It has a global minimum at point ( , )= ( 1 , 1 ), where ( , ) = 0.
Launching each of the software implemented methods has shown that all of them find a global minimum when setting different levels of accuracy. Thus, the following accuracy parameters were used for gradient descent methods: 1e-6. As can be seen from the results, the required accuracy of calculations is achieved, the result is unambiguous. The comparison of the number of iterations cannot provide an unbiased assessment of the amount of performed calculations, because the method of the steepest descent and that of conjugate gradients, in addition to the basic iterations, perform a one-dimensional step length search, which also includes the evaluation of the function, while the method of differential evolution on the first iteration forms a population of 20 random vectors, after which there is a change of generations reaching up to 300 of them. Thus, the most unbiased way to compare these methods is to estimate the number of calls to the calculations of the objective function value. However, one should also take into account that besides evaluating the function, gradient methods also calculate the gradient (that of the derivative of each variable), which also makes the calculations heavier and greatly affects the speed of the algorithms. At the same time, the method of differential evolution calculates only the values of the objective function.
At the second stage of the study, the efficiency of the software implemented methods was assessed as applied to the problem of minimizing the power shortage. The studies were conducted based on mathematical model (1)(2)(3)(4)(5)(6). The 3-zone system was chosen as the system to be tested test, which is an isolated EPS with the "ring" topology (see Figure 1), made up of three zones of reliability and three inter-zone links, presented in the figure as nodes of the graph and its edges, respectively. = 10 was chosen as the initial parameter of the penalty function method with its subsequent increase up to 1,000 with the 10 times increase step. However, as evidenced in practice, the change of coefficient with its subsequent increase is required only for the conjugate gradient method that is gradually approaching the desired solution, while in the case of the method of differential evolution it is enough to specify the coefficient once to obtain the final solution.
The results of the performance of the conjugate gradient method in one of the numerous experiments are shown in Table 3.  Table 3, the solution matches the required correct values within the permissible error of 1% that occurs due to the numerical instability of the penalty function method and the error of representation of real numbers in the computer memory.
The following values were used as setup parameters for the method of conjugate gradients: the accuracy of the calculations: 10e -6, the restart point was set at 18 iterations (chosen experimentally), the maximum number of iterations was 50,000. Then, an experiment was carried out to calculate the solution to the problem of minimizing the power shortage with various starting points. The minimum number of method iterations was 1,422; the maximum number of method iterations was 6,793; on average, 3,776 iterations were required to obtain a proper solution, the average number of calls to the evaluation of the objective function was 46,332 times, and the average gradient was calculated 3,776 times, which corresponds to the average number of iterations.
After that, the method of differential evolution was tested with the penalty coefficient = 10 set as setup parameters, the number of populations was 96, the mutation coefficient = 0.5 (this coefficient was chosen experimentally and provides the best results of the convergence rate of the method), the coefficient of crossover speed was 0.9, and the maximum number of generations amounted to 1,500.
As a result of the application of the method of differential evolution, the solution was obtained using only one penalty value = 10. The final solution is presented in Table 4. 0,00005 0,00005 As can be seen from Table 4, the solution also matches the required correct values within the permissible error of 1%, which occurred due to the error of representation of real numbers in the computer memory and the numerical instability of the penalty functions method. The minimum number of method iterations was 680, the maximum number of method iterations was 1,500, on average 943 iterations were required to get a "proper" solution, the average number of calls to the evaluation of the objective function was 56,688 times, but such operations are deemed to be lightweight and do not require heavy additional calculations.

Conclusion
The speed and accuracy of solving the problem of minimizing the power shortage have an effect on obtaining adequate values of the EPS reliability indicators and on further solving of the subsequent problems, for example, the justification of the EPS redundant generating capacity. The problem of minimizing the power shortage is non-linear and non-convex, so it is required to apply appropriate methods to solve it. The paper treats the issue of applicability and efficiency of the solution by methods of conjugate gradients and differential evolution. Experimental studies have shown that both methods enable us to solve the problem with a given accuracy, but the effort for solving this problem is different. Since it is not sound to compare methods by the number of iterations spent on arriving at the solution, because the effort taken by iterations of each method varies greatly, it was suggested to compare the number of evaluations of the objective function. Backed by the analysis of the performance of the methods, distinguishing features of the method of differential evolution may be highlighted. Namely, there is no need to take into account the constraints of the maximum and minimum values of the level of generation, loads, and bandwidth at the level of penalty functions, since these constraints are addressed at the level of the method itself. Furthermore, there is no need to increase the value of the parameter of the penalty with the subsequent calculation, because the method works correctly with the value of the penalty equal to 10. Despite the fact that the method of differential evolution evaluates the objective function 10,356 more times than the conjugate gradient method, it is necessary to take into account that such calculations are more lightweight, because the objective function is free of some of the penalties, and there is no need to calculate the function gradient at each iteration.
Thus, the method of differential evolution was determined as the most efficient among the implemented methods in terms of accuracy and resources spent on calculations, which is confirmed by numerical experiments performed for small systems.