Analysis power shortage minimization methods in the modern processing software for adequacy assessment of electric power systems

. Analysis of domestic and foreign software systems for assessing the resource adequacy showed a variety of models and methods used in them. Many software systems use both linear and nonlinear models, these models are optimized according to various criteria to simulate the operation of the system. As tools for solving, software usually use commercial high-level modelling systems for mathematical optimization. However, in addition to the existing ready-made commercial solutions, the authors consider the effectiveness of optimization methods, as well as their parallelized versions, which can be independently implemented and applied as a solver for a specific problem. As a result, it was confirmed that these methods can be used to solve the problem, but they are less effective relative to a commercial solver. From the point of view of accuracy and resources spent on calculations, the most effective of the independently implemented methods turned out to be the parallelized method of differential evolution, which was confirmed by numerical experiments on small systems.


Introduction
The need to ensure a high level of reliability of electric power systems (EES) has always been relevant, however, modern consumers of electricity impose even higher requirements for ensuring the reliability of power supply. This is also associated with a high level of economic costs in the event of an interruption in the production and supply of electricity to consumers, and directly depends on the equipment and its failures. In view of constant changes towards the enlargement and complication of EPS, early planning, timely correction of changes and redundancy of its elements is one of the main directions of ensuring reliability. Thus, to minimize the number of cases of limiting the supply of electricity to consumers, it is necessary to implement in advance a set of technical and organizational measures aimed at increasing the reliability of the EPS. However, due to the fact that such activities are costly, an objective justification for their implementation requires a qualified assessment. For this, an assessment of the resource adequacy of prospective EPS schemes is carried out. The result of the assessment is reliability indicators that have an economic interpretation.
One of the stages of assessing the resource adequacy when applying the Monte Carlo method [1] is to determine the power deficits of the possible states of the EPS. The basis for calculating power deficits is the simulation of EPS, which includes a mathematical model of the EPS, as well as optimization methods that allow obtaining the value of the power deficit 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 a growing number of optimized parameters, depends on the optimization method used and the correctness of the mathematical model. The statement of the problem of minimizing the power deficit can be presented both in linear and nonlinear form [2]. The most adequate formulation is in a non-linear form, where the losses in power lines have a quadratic dependence on the transmitted power [3].
In the well-known domestic and foreign practice, various optimization methods are used to search for power shortages, for example, in the YANTAR software and computing complex [4][5], the method of internal points is used in conjunction with a linear and nonlinear model, in the ORION-M " [6] a dual simplex method and a linear model are used. In the software " Reliability "a high-level modeling system for mathematical optimization" GAMS "(CONOPT nonlinear optimization solver) is used, as well as various gradient and heuristic methods are investigated, linear and nonlinear models are considered [7 -8]. In turn, foreign software use both ready-made commercial solvers and their own implementations. For example, the GRARE complex uses the implemented methods of linear and quadratic programming to solve a linear problem [9], ROM uses the GAMS system and considers only nonlinear models [10], the ANTARES complex uses independently implemented linear models and optimization methods [11], in several in the PLEXOS complexes -UCo2, CCo1, EU 2030 [12], linear and nonlinear models are implemented, the Gurobi, CPLEX, XPRESS solvers are used, all of the above complexes are commercially available, closed projects.
Based on the proposed for consideration software computing complexes for assessing the resource adequacy, it can be concluded that the developers of the complexes use both independently implemented optimization methods and off-the-shelf products in the form of linear and nonlinear optimization problem solvers. The presented work analyzes the effectiveness of various independently developed optimization methods and their versions with embedded parallel technologies, as well as available commercial solutions.

Interior point method
Interior point (IP) methods (also called barrier methods) are a specific class of algorithms for solving linear and nonlinear convex optimization problems [13].
The key idea of the interior point algorithms is to eliminate constraints -inequalities from the problem -by introducing a quadratic or logarithmic penalty in the objective function for approaching the boundaries of the admissible region. According to the methods, the starting point for the search can only be selected within the valid area. The choice of the starting point of the search is carried out depending on the formulation of the problem. In the absence of constraints or converting them to penalty functions with an outer point, the starting point is chosen arbitrarily. When constrained or converted to penalty functions with an interior point, the starting point is selected within the valid range. In this case, the set of points is divided into acceptable and unacceptable, depending on the restrictions. In turn, the set of admissible points, depending on the constraints, is also divided into boundary and internal. A function ∶  → ℝ is called a barrier function for a set  if ( ) → +∞ with → , where  is the interior space of  and  is the boundary of . Instead of the original problem, it is proposed to solve the problem: min ( , ) = ( ) + ( ). (2.1) and are given only in , the barrier property guarantees that by minimum of is exist, at the same time, the larger provide the greater influence of . Under reasonable conditions, it can be achieved that if would tends to infinity, then the minimum of converges to the solution of the original problem.
If the set is given as a set of inequalities ( ) ≤ 0, 1 ≤ ≤ , then the standard choice of the barrier function is the logarithmic barrier: (2. 2) The minimum points * ( ) of the function ( , ) for different forms a curve, which is usually called the central path, and the interior point method tries to follow this path.

Conjugate gradient method
The conjugate gradient method (CG) is an iterative numerical method of unconstrained optimization in a multidimensional space [14][15][16]. The main idea of the CG method, like other gradient methods, is to descend in the direction of decreasing the gradient, but the search for the desired step and trajectory is carried out using conjugate directions.
The definition of conjugacy is formulated as follows: two vectors x and y are called A-conjugate (or conjugate with respect to the matrix A) or A-orthogonal if the scalar product of x and Ay is equal to zero, that is: Conjugacy can be considered a generalization of the concept of orthogonality. Indeed, when the matrix A is the identity matrix, in accordance with equality (3), the vectors x and y are orthogonal. One possible way to calculate conjugate directions is to use the iterative method to calculate conjugate directions, Fletcher-Reeves: Expression (4) means that the new conjugate direction is obtained by adding the antigradient at the turning point and the previous direction of motion, multiplied by the coefficient calculated by (5). The directions calculated by (4) turn out to be conjugate if the function to be minimized is given in the form: That is, for quadratic functions, the conjugate gradient method finds the minimum in n steps (n is the dimension of the search space). For general functions, the algorithm ceases to be finite and becomes iterative. At the same time, Fletcher and Reeves propose to restart the algorithmic procedure every n + 1 steps. A restart is necessary in order to forget the last direction of the search and start the algorithm again in the direction of the fastest descent.
In this paper, in addition to the work of the original algorithm, a version with the use of parallelization technology is also considered. As mentioned above, for the algorithm to work, it is necessary to zero the last search direction every n + 1 steps, i.e. restart of the procedure. However, this parameter is often selected manually and the rate of convergence of the method, as well as obtaining the result, can vary greatly depending on the chosen restart moment. Thus, not always, a once selected iteration number for zeroing makes it possible to find the optimal solution efficiently in terms of resource and time use. As an experiment, 5 different times of direction zeroing were chosen, and indeed the convergence rate of the method varied greatly, including if the numbers of zeroing iterations changed during the execution of the algorithm. Proceeding from the fact that the rate of convergence of the method depends on the moment of zeroing the direction, it is proposed to investigate the operation of the conjugate gradient method with the condition of parallel launching of two or more calculations with different times of zeroing, as well as synchronized stopping, evaluating the obtained solution and choosing the solution as the starting point, with the smallest value of the gradient norm.

Differential evolution method
Differential evolution method (DE) -is a metaheuristic method of multidimensional mathematical optimization, belongs to the class of stochastic optimization algorithms and uses some ideas of genetic algorithms, but does not require working with variables in a binary code [17][18][19]. The method (DE) uses the generation of a certain set of chromosomes or vectors of parameters = { 1 , … , }, called a generation, the size of which does not change during the execution of all operations. At each iteration, these vectors are updated by generating a new generation obtained from the previous one using a special procedure that combines the operations of mutation and crossover.
During the existence of the method, many of its modifications were invented, in particular, they concerned the procedure for the formation of new generations. In the implementation under consideration, the essence of this procedure is as follows: in the previous generation, the target vector R1 is determined, as well as 3 different, random vectors (R2, R3, R4), then the mutation coefficient F is set (a real number in the interval [0,1]), after that the vector "mutant" is formed according to the following expression: (2.7) Then the target vector and the mutant vector are crossed, depending on the probability of mutation.
This method only requires the ability to calculate the values of the objective functions, without taking into account the upper and lower constraints on the variables, i.e. these restrictions are already use in the method. DE is designed to find the global minimum (or maximum) of non-differentiable, nonlinear, multimodal (a large number of local extrema) functions of many variables. The method is simple to implement and use (it contains few control parameters that require selection) and in reality is an iterative process of creating new, improved generations with mutations, where the end is finding the best set of chromosomes. Despite this, the differential evolution method can be considered one of the most susceptible to parallelization. The implementation of this mechanism consists in the parallel formation of new descendants, since here they are not directly dependent on each other; the mechanism can have the character of data parallelization. Thus, the process of forming each new vector in a generation can be taken out into a separate problem and solved independently of the others, at each iteration. At the end of the procedure for forming the entire population, a process of stream synchronization is required to start the next iteration or stop calculations.

A set of optimization methods for the CONOPT package
Many foreign and some domestic computational systems are based on commercial solvers, it is allowed to move away from the direct implementation of optimization methods, but such complexes are no less interesting than their own developments. This paper discusses the use of a high-level modeling system for mathematical optimization "GAMS" [20], this complex is aimed at solving various mathematical problems, including optimization ones. GAMS automatically connects the CONOPT4 solver as a nonlinear optimization problem solver. This solver has a hided code and a complex system for transforming the problem into a simplified form, and also divides the solution of the problem into several stages where several optimization methods interact.
The first step begins with a "starting point setting" procedure, changing individual variables one at a time. The procedure is based on Newton's method with some heuristic modifications. Next, there is a transition to the next stages, where the model is solved by the method of sequential linear programming (SLP), after which, the descent by the method of conjugate gradients occurs, then the solution is started by the method of sequential quadratic programming (SQP).
Successive Linear Programming (SLP), also known as Sequential Linear Programming, is an optimization technique for approximately solving nonlinear optimization problems.
Starting at some estimate of the optimal solution, the method is based on solving a sequence of first-order approximations (i.e. linearizations) of the model. The linearizations are linear programming problems, which can be solved efficiently. As the linearizations need not be bounded, trust regions or similar techniques are needed to ensure convergence in theory.
Sequential quadratic programming (SQP) is one of the most common and effective general-purpose optimization algorithms, the main idea of which is the sequential solution of quadratic programming problems that approximate a given optimization problem. For optimization problems without constraints, the SQP algorithm is transformed into Newton's method of finding a point at which the gradient of the objective function vanishes. To solve the original problem with equality constraints, the SQP method is transformed into a special implementation of Newtonian methods for solving the Lagrange system.

Mathematical formulation of the problem of minimizing power deficit
The problem of minimizing the power deficit is formulated as follows: for known values of operable generating capacities, required levels of consumer loads, transmission capacities of EPS links and power loss coefficients in EPS links, it is necessary to determine the optimal flow distribution in the EPS [1], [4][5]. There are several types of models for minimizing power shortages; in this article, the applied models will be considered. The following is a linear formulation of the problem: Mathematically, the problem is formulated as follows: when the balance constraints are respected: -specified positive coefficients of specific power losses during its transfer from zone j to zone i, j≠i, = 1, … , , = 1, … , .
The considered model (3.1-3.6) is a common flow distribution model in the field of assessing balance reliability, the solution of which is carried out by minimizing the power deficit. Model (3.1-3.6) is a transport problem. To solve the presented optimization problem, in view of its relative simplicity, the simplex method and the dual simplex method in their various variations are mainly used.
Based on the linear model (3.1-3.6), which is a flow distribution model, it is possible to apply a modification of the balance constraint, so in [4][5] there is a reasonable conclusion that a model where power losses depend on the square of the transmitted power is a more adequate model, close in physical sense to a real model. For this, the model  Thus, the set task can be presented in two formslinear and nonlinear programming problems. The type of the problem strictly depends on the used balance constraints of formulas (3.1), (3.3-3.7). These models are actively used in the "Reliability" computer complex, and were also previously used in the "YANTAR" software.

Analysis of optimization methods used in software systems
Experimental calculations were carried out for systems of various configurations, however, in this work, a comparative analysis of optimization methods and the GAMS complex used in different software, applied to a computational test system with the same initial information, is considered, then a test system (TS). Previously, this test system has already been considered for calculating and checking the correctness of models and applied optimization methods in the YANTAR software [4][5], as well as for their research [5], this in turn will add to the comparison methods of internal points implemented by the authors of this complex. Also, for calculations, manually and automatically, a model was implemented, formed from this scheme (Fig. 1) within the framework of the GAMS complex. To test independently implemented optimization methods of the "Reliability" software, this system was created automatically from the descriptive part of this scheme. The main characteristics of the reliability zones of the test system are shown in Table 1. It is worth noting that some reliability zones have excess generating capacity (4,6), and some are in short supply (1,2,3,5,7). In total, the system has 7 links, the throughput of each link in the directions was the same, in the forward and backward directions. Also, each link had an individual loss rate.
For an indicative comparison of the effectiveness of certain methods, a computational test system was used (Fig. 1), with the characteristics described in Tables 1 and  2. Information about the system and its characteristics, about the work of the method of interior points using quadratic approximations (IPQA) and based on linearization (IPBL) with different calculation accuracy was taken from the source [5]. The calculations obtained using the high-level modeling system for mathematical optimization "GAMS", namely using the CONOPT solver (GAMS_C), as well as the methods of conjugate gradients (CG), differential evolution (DE) and their parallelized versions (PCG, PDE) were carried out independently the authors of the article. Despite the fact that the calculations were carried out on the presented scheme with identical characteristics, the most relevant are the calculations carried out by the authors themselves. A PC with the following set of characteristics was used as the testing equipment: CPU Intel (R) Core i7-8700K @ 3.70GHz, boost 4.50GHz, RAM DDR4 48.0 GB, 15/15/15/36, 2133 MHz. In this case, the only objective and available parameter for comparing these methods is the time to solve the problem. The obtained calculation results show the high efficiency of the methods used in the CONOPT solver. The high speed of the solution is achieved by using various approaches and technologies, including simplifying the resulting model, reducing the number of required parameters by expressing them through other variables. The speed of work of independently implemented methods is inferior to a commercial solver. However, due to the fact that the mechanisms used in a commercial solver are not available for open study, it is difficult to estimate the total amount of work done during the period of solving the problem, including taking into account the accuracy of the solution. Unfortunately, the comparison of these methods by the criterion of the number of iterations is incorrect, due to the fact that iterations of different methods can have a different amount of computation. For example, the SG method in one iteration calculates the value of the function, additionally calculates the value of the gradient and several times the value of the function during onedimensional optimization, while the DE method calculates in one iteration only the value of the function for each element of the population.
Despite the fact that the calculations by the DE method last about the same as the calculations by the SG method, the parallelized version of the PDE using 6 streams received an almost 6-fold reduction in the calculation time. This speaks of the efficiency of using parallelization technologies, and in turn allows using more efficient systems, with a large number of threads, and performing calculations faster.

Conclusions
The speed and accuracy of solving the problem of minimizing the power deficit affect the obtaining of adequate values of the EPS reliability indicators and are a necessary condition for solving the problem of substantiating the value of the power reserve. This paper considers various methods and systems used in modern complexes for assessing the resource adequacy of EPS, as well as the application of parallelization technologies to the methods of conjugate gradients and differential evolution. As an experiment, the calculation of a test 7node system was carried out using the method of internal points, methods of conjugate gradients and differential evolution, as well as using the CONOPT solver from the GAMS package, which is popular in foreign complexes for assessing reliability. Experimental studies have shown that the CONOPT solver from the GAMS package turned out to be the most efficient in terms of the speed of solving the problem, however, when using powerful multiprocessor systems, the differential evolution method in its classical design can be no less effective. Thus, it should be noted that the parallelization mechanisms have shown their efficiency in this situation can significantly reduce computation time and equipment downtime.