Simulated Annealing, Differential Evolution and Directed Search Methods for Generator Maintenance Scheduling

: Generator maintenance scheduling presents many engineering issues that provide power system personnel with a variety of challenges, and one can hardly afford to neglect these engineering issues in the future. Additionally, there is vital need for further development of the repair planning task complexity in order to take into account the vast majority of power flow constraints. At present, the question still remains as to which approach is the simplest and most effective, as well as appropriate for further application in the power flow-oriented statement of the repair planning problem. This research compared directed search, differential evolution, and very fast simulated annealing methods based on a number of numerical calculations and made conclusions about their prospective utilization in terms of a more complicated mathematical formulation of the repair planning task. A comparison of results shows that the effectiveness of directed search methods should not be underestimated, and that the pure differential evolution and very fast simulated annealing approaches are not essentially reliable for repair planning. The experimental results demonstrate the perspectivity of unifying single-procedure methods in order to net out risk associated with specific features of these approaches. and V.P.O.;


Introduction
Currently, there is a significant national and industrial push for a reduction in emissions, as well as for reductions in reliance on fossil fuels as a whole and conventional generation in particular. According to this trend, the development of new power system control approaches on the one hand, and wide integration of renewable generation on the other, have shown an unprecedented pace in recent years. Modern demand response strategy, electric vehicle incorporation, increasing penetration of renewable generation, and energy storage systems only provide higher equipment and grid efficiency, environmental friendliness, and electrical power system flexibility. Nevertheless, the same novelties also pose new significantly complicated problems for reliable and economic control of power systems. In fact, the uncertain characteristics of primarily renewable generation and load profiles equipped with energy storage devices may develop into serious and, at times, disastrous situations, causing, first and foremost, economic losses for electrical industry subjects.
Flexible and optimal control of the generating equipment state and output power is the primary effective means of power system state control for system operators. As a result, the problems mentioned above increase the importance of every generating unit operating in power systems. The concern regarding the controllability of generating equipment underlines, in turn, the significance of generator maintenance scheduling procedures. These measures can be completed using generating adequacy approaches, being responsible for the assessment of the balance between generated and consumed energy [1][2][3]. For instance, the criterion being widely used in order to estimate the generating adequacy of a power system is expected an energy not supplied (EENS) coefficient [4].
The scientific community recognized the risk and seriousness of inappropriate repair planning long before this research was completed. Despite this fact, there is still no general approach to major and midlife maintenance planning for generating equipment. The primary barrier for a successful solution is the complicated integer entity of the optimization problem, being multiextremal and requiring additional measures to find the global optimum or a close quasi-optimum point. Then, all of the possible approaches to the problem could be subdivided into deterministic and heuristic ones. Different researchers advocate the effectiveness of the certain methods corresponding to these groups. As a result, the development and implementation of specialized generating equipment maintenance scheduling methods are needed, which are discussed below in detail.
Currently, there is a rich body of specialized literature related to the application of integer programming methods for generator maintenance scheduling [5][6][7][8]. There is certain disagreement over how adequate these types of methods are, since they were one of the first approaches used to solve the problem under consideration. Nevertheless, the idea does not correspond to the current situation, as, for example, researchers in Ref. [6] and Ref. [7] propose to apply integer programming methods for generator maintenance scheduling, taking into account demand response strategies and AC/DC high voltage directed current (HVDC) link constraints, respectively. However, many researchers, as it was mentioned above, advocate the use of heuristic approaches to solve the problem-in particular, increasing attention today is being paid to ant colony algorithms [9][10][11], fuzzy and artificial neural network approaches [12][13][14], evolution [15,16] and genetic methods [17][18][19], simulated annealing methodologies [20][21][22][23], particle swarm techniques [24][25][26], taboo search procedures [27][28][29], etc. In addition, one should not forget about directed search methods [30], which combine deterministic nature with the simplest implementation, while ensuring appropriate planning results.
Taking into account the diversity of the presented possible approaches, a question arises as to what method is the most appropriate for further application in a more complicated scheduling statement. Currently, it is an acceptable practice to consider only active power line loading constraints [6,7,21]. Nevertheless, in order to plan repairs or outages, for example, for either static reactive compensators or energy storage devices, it is necessary to consider more complex constraints, such as allowable voltage levels and small-signal stability limitations of a power system. The need becomes dramatically increased when equipment to be scheduled operates in the weak networks that still exist and operate widely around the world. From this point of view, there is a significant need for comparison of the different currently applied methods in order to develop existing or create new methods that are simple enough, reliable, and fast in their calculations.
It is proposed, in terms of research, to compare differential evolution and very fast simulated annealing methods as modifications of the genetic algorithm and convenient simulated annealing methods, respectively. The directed search method is to be used as a reference method in terms of this comparison. Despite their heuristic characteristics, annealing and evolution methods belong to two different approaches to the solution searching process: the simulated annealing method is very similar to the probabilistically modified directed search method, which means that on every step of the iterative process, there is always one current solution point and a certain potential variant to be assessed. Meanwhile, the differential evolution approach relies on the diversity and significant quantity of current solutions unified in a population, and looks somewhat similar, from this point of view, to the particle swarm approach. As a result, a comparison of these methods would provide not only profitable information about their own effectiveness, but could also uncover the most appropriate methods for generator maintenance scheduling.
A comparison is based on the resultant optimization function values for a number of test systems, the simplicity of the implementation of the method, the calculation time needed to complete scheduling, the possibility to take into account power flow constraints, and the number of objective function calculations. Analysis of the results and further perspectives of the use of these methods for complex power system equipment maintenance scheduling procedures is not merely a matter of academic interest, but is vital in industrial application. This is reasonable, as the resultant collected data can vary significantly depending upon the dimensionality of the problem, the proportion of generators to be repaired in a power system, the generation and consumption profile, and so on.
In addition, currently there is a very limited number of publications considering the effectiveness of the initialization and initial parameter settings for the heuristic methods of maintenance scheduling procedures. In particular, such research as in Ref. [15] for the differential evolution method and Ref. [20] for the simulated annealing procedure estimate primarily the possibility to apply these methods but not behavior of the results in different scenarios. In this research, the initialization procedure for the simulated annealing method is presented. The generation profile diversity and the power system's reserve effect on the maintenance schedule are also assessed for simulated annealing and differential evolution methods. Moreover, the influence of the different parameters of these methods on the resultant maintenance plans is assessed in order to justify a certain approach for their initialization, as well as to analyze the behavior of the simulation results.
The rest of this paper is organized as follows: Section 2 presents the mathematical formulation of the maintenance planning problem statement considered and a description of the proposed implementation of the directed search, differential evolution, and very fast simulated annealing methods, as well as providing descriptions of the experimental procedure and the test model. Section 3 represents the discussion of the comparison results, while Section 4 concludes this paper.

The Methods under Consideration and a Description of the Test Materials
This section provides insight into the mathematical formulation of the generator scheduling problem, the implementation of the directed search, differential evolution, and very fast simulated annealing generator maintenance scheduling methods, and the test data used in terms of the experimental part of this research.

Mathematical Formulation of the Repair Planning Problem
First of all, attention should be paid to the mathematical statement of the generator maintenance planning problem.
Task. For all generating units to be repaired within the time interval being considered, set moments of the beginning of the repairs, taking into account the optimization criterion.
Mathematical method. In this case, the very fast simulated annealing, differential evolution, and directed search methods are applied.
Variables. Variables are discontinuous moments of the beginning of generating unit maintenance ∈ N , which can be represented as vector ̅ = , , … , . Planning horizon. The time period within which all repairs should be completed is stated to be equal to 1 year, while the simulation step is 1 week-= 52. It is suggested that the generation profile remains constant during the week.
Optimization criterion and objective function. The optimization criterion is generating adequacy oriented in character. It is based on assessment of the EENS coefficient. The idea behind the criterion can be described as follows: as a system for which generation maintenance planning is to be completed is assumed to be strongly interconnected, the power state constraints and topology of the system could be omitted in the first approximation. Therefore, an undersupply of energy can arise only when the electrical load in the system is greater than the total available generation for a certain moment. Then, let us observe an equation for EENS for one hour, taking into account the probability characteristics of the system consumption: where ( ) is the EENS for hour t of week w, MWh; is the expected electrical energy consumption for hour t, MWh; , is the standard deviation of consumption for hour t, MWh; is the total available generation for hour t, MWh; and , , , and ( , , , ) are the consumption cumulative distribution and the probability density function for hour t, respectively. Obviously, taking into account that the calculation step is equal to 1 week and the fact that the generation profile remains constant during this interval, the EENS coefficient calculated using Equation (1) reaches its maximum value when is the peak load value for week w. As a result, the generating adequacy of the system for week w can be assessed based on the value only.
From the above reasoning, it can be concluded that the generating adequacy optimization criterion can be stated as the minimum of the sum of the EENS values for the peak hours of every week w within the time period under investigation. Then, the objective function for the period of 1 year can be written as follows: where is the EENS for peak load hour of week w.
Constraints. In order to integrate the constraints into the problem statement, it is necessary to introduce new auxiliary variables , for all generators to be repaired i and in which weeks w. If generator i is undergoing maintenance in week w, then , = 0; otherwise, , = 1. Using these parameters, it is stated that the objective Equation (2) is subjected to the following constraints:

•
The maintenance of generator i to be repaired is to be continuous over repair: • The maintenance of generator i to be repaired is to be completed within period T, while the duration of repair procedure is equal to : where is the week in which the maintenance of generating unit i is started and is the duration in weeks of the maintenance action of generation unit i; • The maintenance of generator i to be repaired is to be started after initial week , and is to be finished before final week : where is the week after which the maintenance of generating unit i is to be started and is the week before which the maintenance of generating unit i is to be finished.

•
As there is a limited capability of maintenance crews in a system, there cannot be more than simultaneous maintenance actions in week w: where is the capability of the maintenance crews in week w (the number of crews that are able to provide maintenance actions individually).
Assumptions. The following assumptions are made in this research. Power state constraints and the topology of a power system are not considered, as the system is supposed to be a strongly interconnected one. It is supposed that maintenance crews are able to complete maintenance actions on all types of units. The economical relationships between the subjects of the electricity market and the economic characteristics of generators are neglected.
In the following sections, on the basis of this mathematical formulation of the problem, a description of the methods under investigation is provided.

Directed Search Method for Generating Equipment Maintenance Planning
The directed search methods proposed first by Oboskalov [30] are methods that rest on a twostep approach: an initial sort of the list of generators to be repaired, and then the local enumeration procedure. The research in Ref. [30] was primarily related to deep investigation into the effectiveness of the directed search method and showed the probabilistic characteristics of the power system generation reserve. In order to provide insight into the nature of the directed search method, it would be profitable to present briefly its principal algorithm in the current research. The algorithm of the second-order directed search method can be described as follows: (1) For every unit i in the list of generators to be repaired, criterial coefficient is calculated according to the following formula: where is the available active power of generating unit i, MW, and is the probability of failure of the generating unit i, p.u.
(2) The list of generators is sorted in descending order of . As a result, after the sorting procedure, the generators at the top of the list are the most reliable, powerful, and long since repaired units.
In particular, these machines will be repaired first in the following stages. It is of note that research in Ref. [30] mentioned that the above proves that the initial sorting procedure of the generator list ensures that generator maintenance scheduling using this method finds one of the best solutions in the top 2% of all possible combinations of moments to begin the repairs. Depending upon the quantity of the generators enumerated in the fourth stage, the order of the method is defined. For example, if there are three generators enumerated in every iteration, then this method is named the third-order directed search method.

Differential Evolution Method for Generating Equipment Maintenance Planning
The differential evolution approach, first proposed by Storn [31] as mentioned above, is a heuristic method closely connected to genetic algorithms. According to the methodology, the solution space is considered a linear one. As an almost immediate corollary, this feature allows application of vector space mathematical tools in order to complete crossing-over operations and to find appropriate solutions. The maintenance planning algorithm on the basis of this method can be represented as follows: (1) In order to obtain an initial estimate for further optimization processing, the first-order directed search method is applied to the list of generators to be repaired. It is fair to ask why the firstorder method is used here; the fact is that the strongest bases on which the differential evolution method rests are the diversity and dispersion of the initial data. Consequently, in this step, there is no need to find the best solution, because the results of the initial optimization are needed only for further localization of the initial maintenance plan vectors in the following steps. This method returns the vector of the repair moments for generators ̅ . (2) A population matrix is formed, comprising column vectors ̅ : where n is the number of generators to be repaired; s is the population size (i.e., the number of repair plan vectors included in the population); and , is the repair moment for generating unit i according to repair plan vector j.
(3) The population matrix (Equation (8)) is changed in order to obtain an initial population matrix according to the following equations: where ∆ is the pseudorandom departure of the repair moment of generating unit i of column vector j, drawn according to a further probability distribution: where = • is the repair moment standard deviation for generating unit j. Since, after operation (Equation (10)), the repair moments do not belong to a set of natural numbers ∉ , the resultant elements of the population matrix must be rounded off to the nearest integer values. The idea behind the initialization of the standard deviation value for the presented probability distribution is that the pseudorandom departure of the repair moment for unit i cannot be greater than half of the duration of the maintenance period for this unit (i.e., -three standard deviations for a normal probability distribution, Equation (10)). This approach to the initialization of repair moment departures, which directly affects the diversity of the initial set of repair plans, can be changed depending upon the user's assessment.
(4) Every iteration of the search process k for population can be subdivided into the following steps: (a) For every column vector ̅ of the population matrix, two more column vectors, ̅ and ̅ , are chosen randomly, while . Using these three elements, a trial vector is calculated as follows: where is the crossing-over coefficient, defined by the user. The coefficient can be represented not only as a certain number, but also as a diagonal matrix.
(b) The child vector is defined based on the trial vector and the parent vector ̅ according to the idea illustrated graphically by Figure 1: two integer numbers, h and g, are chosen randomly (ℎ , ℎ ∈ 1, , ≤ − ℎ) , while h is the starting index of the interval to be exchanged by the elements of the same interval of , and g is the number of elements that are to be exchanged.  Despite the resultant vector ̅ being the purpose of the whole procedure, the final population matrix consists of many suboptimal repair plans and, great number of such vectors is considered during the whole iteration process. This is a distinctive feature of the differential evolution method as, for example, these elements are valuable as possible starting points for further applying the methods strongly influenced by the initial conditions. Finally, this method seems perfect for the first step of a bee algorithm.
To conclude, a question remains as to how to choose the parameters of the method, particularly the population size s and the crossing-over coefficient . Moreover, here one may observe another significant advantage of the method. First, although these parameters influence the calculation results, it is intuitively obvious that a greater population size has a positive effect on the resultant generator maintenance schedules. Meanwhile, from purely a logical point of view, the appropriate crossing-over coefficient values are localized in a comparatively narrow range of possible values from 0.1 to 1.0. s, speaking only about the crossing-over coefficient, it should be noted that initial diversity of the population guarantees that the method will provide at least a number of acceptable maintenance plan variants, regardless of the coefficient chosen. As a result, there is no need for a special approach to define the certain optimal setting of the method. Thus, its parameters could be set intuitively or on the basis of the tentative results of preliminary experiments.

Very Fast Simulated Annealing Method for Generating Equipment Maintenance Planning
The convenient simulated annealing method was invented by Metropolis [32] and developed by Kirkpatrick [33]. The ideas involved in the method are related to a mathematical representation of the physical metal annealing process. According to the conventional simulated annealing method, the temperature value is established for a certain system of elements. On every iteration of the optimization process, a new system state is found and its energy , which represents an objective function, is calculated. If a new system state is better, then it is chosen as a new basic state for the following process and the temperature is decreased; otherwise, depending upon the temperature value, the system could potentially move to a worse solution point. The last mechanism allows the avoidance of local extremum points. The very fast simulated annealing method, first proposed by Ingberg [34], is an improved version of this algorithm, providing faster convergence of the iterative process.
The proposed algorithm for generating equipment maintenance planning by application of the very fast simulated annealing method, with the first-order directed search method as a means to obtain an initial estimate of the vector of repair plan moments, can be represented by the following steps: (1) By applying the first-order directed search method, an initial estimate of the vector of repair plan moments ̅ is completed. (2) The objective function value ( ̅ ) is calculated according to (2) in order to estimate the initial repair plan ̅ . The initial vector of the repair moments ̅ is considered as the base one; therefore, for the first iteration, ̅ = ̅ . (3) Every iteration k of the searching process can be subdivided into the following steps: (a) The vector of the repair moment departures ∆ ̅ is calculated according the following equation: where is the independent random variable, corresponding to generation unit i, uniformly distributed on the interval [0,1]; is the temperature, corresponding to generating unit i; d is the divider parameter needed to reduce the dispersion of the repair moments in terms of a solution search procedure; , is the last week when it is allowable to start repairing the generating unit i; and , is the very first week when it is allowable to start repairing the generating unit i. The fundamental idea behind divider is that the methodology proposed by Ingberg does not completely correspond to the specific characteristics of the maintenance planning problem. In particular, the use of the pure difference , − , without a divider leads to significant repair moment departures ∆ ̅ , meaning that further local searching becomes meaningless. However, the basic idea on which the approach lies is to complete a local search of the possible appropriate variants. As a result, the divider parameter is used to localize the solution search process in the neighborhood of a starting point by limiting the random repair moment departures in every iteration.
(b) A new vector of the repair moments ̅ is defined using the departure vector ∆ ̅ : (c) The objective function value ( ̅ ) is calculated for the new vector of the repair moments ̅ according to Equation (2).
, the base vector is equated to ̅ : and for all dimensions of the repair moment vector, the temperature is updated according to the following equation: where is the initial temperature for dimension i of the solution state space; is the temperature coefficient; is the iteration number; and is the number of generators to be repaired.
(e) If ( ̅ ) ≥ ( ̅ ), the independent random variable , uniformly distributed on interval [0; 1], is drawn for every generator i to be repaired and in the case when the repair moment and temperature for the generating unit i are redefined according to Equations (14) and (15), respectively. Otherwise, these parameters remain unchanged.
(f) The number of iterations k is increased by 1 and the calculation process continues from step 3a. (4) The iteration process is stopped if the user-defined minimal temperature or maximum number of iterations is reached. Nevertheless, it is fairly difficult to assess precisely the minimal temperature without initial experimental test results following application of the method. Consequently, it is reasonable to follow the iteration maximum number criteria in order to simplify the research matter. (5) The resultant vector of the repair moments ̅ is chosen among all system states considered within the iteration process, as a vector with the smallest objective function value ( ̅ ).
The disadvantage of the very fast simulated annealing method, widely referred to in the specialized literature, is related to the significant number of settings needed to operate with this tool (e.g., for the case when generator maintenance scheduling for 30 units is required, there are more than 60 method parameters to be established before its application). The initialization procedure, as a result, requires a long time initial calculations. In order to avoid this complication of the calculation process, it is assumed that the temperature and its coefficient are the same for all dimensions of the problem. This step essentially allows simplification of the simulation process, since only two parameters (i.e., general temperature and its coefficient ) are to be defined, disregarding the dimensionality of the problem.
Nevertheless, even when using the simplified simulated annealing approach, a question remains as to how the method parameters are to be initialized. This problem is caused by the fact that, in contrast to the differential evolution method, the effect of these coefficients on the simulation results is not completely evident. Neither of them can be initialized or considered independently. Moreover, there is one more important distinction between the differential evolution and the simulated annealing methods, which is that inappropriate coefficient values for the latter leads to a lack of appropriate optimization results. Then, it becomes clear that an initialization methodology is needed to operate the very fast simulated annealing method in practice. The methodology is proposed below.
The main idea behind the initialization procedure is as follows. If one considers the dependence of the repair moment departures (Equation (12)) upon the temperature value being equal for all dimensions, it is apparent that the parameter does not affect the maximum and minimum values of said dimensions. On the contrary, it affects only the trace form of the probability distribution. As a result, in order to set the parameters of the method, one could operate with only Equation (17) to determine the probability of accepting a worse solution equation : In order to ensure the effectiveness of the simulated annealing algorithm, it is necessary to set such an initial temperature that, for the mean objective function deviation (∆ ) , the probability of accepting a worse solution is adequate during the simulation process. This means that it is to be relatively high at the very beginning of the searching process and then gradually decreases to significantly low values by the end. The temperature coefficient is established in the second stage using the Fibonacci search method [35]. Let us consider briefly this algorithm: (1) The approximate mean objective function deviation (∆ ) during the process is assessed in order to define an order and a range of its possible deviations. (2) Based on the ranges of possible values of ∆ , and the dimensionality of the problem's dependence of on (∆ ) and are worked out. Using the resultant surface illustrated by Figure 2a, it is necessary to choose such a temperature value that the (∆ ) trace corresponding to it lays on the imaginary border between rapidly and slowly decreasing ones. In fact, further definition of the temperature coefficient partly provides some protection against possible mistakes in this step. Hence, it is just necessary to avoid too low a temperature being established. For instance, in Figure 2a, the acceptable initial temperature value is = 1.5. (3) Using the set of traces for dependence of on the temperature coefficient for several objective function deviation ∆ values, while is equal to the value established in the previous step and the iteration step is equal to 1 ( = 1), the range of possible temperature coefficient values is established. A graphical representation of these traces is provided in Figure 2b. A possible temperature coefficient minimum border is to correspond to the maximum values of the probability of accepting a worse solution , while the maximum border is to correspond to minimum but non-zero probability values. Such a range of values guarantees that at the very beginning of the search process, the probability is high enough for the initially flexible and diversiform development of the simulation process. For instance, in the case of Figure 2b, a further temperature coefficient search would be quite correct if localized within [0.1; 3.0] borders. (4) Using the Fibonacci search method, a local search of the optimal temperature coefficient is organized within the range established. In every step of the search process, the number of simulated annealing generator maintenance scheduling procedures are completed for chosen in step 2 and the current temperature coefficient value . In order to overcome an obstacle connected to the random character of the planning results calculated using the simulated annealing method, the selection criterion is proposed to be the minimum value of the mean objective function for a certain parameter . The Fibonacci search is stopped when the length of the range of possible values reaches the user accuracy of calculation .
After these computational and analytical procedures, one could find appropriate initial temperature and temperature coefficient values for further generator maintenance scheduling using the very fast simulated annealing method.

The Experiment Methodology and Test Data Used
Let us briefly consider the test data and then, on the basis of this information, move to features of the experimental methodology of this research.
As a trivial basis for further comparison and analysis, it is proposed that three variations of the reliability test system (RTS) RTS-96 represented in Ref. [36] are used. The main difference between these three systems is in their scales: The first test system (1TS), graphically illustrated in Figure 3, is the smallest one with a short list of generators to be repaired and a relatively low consumption, in contrast to the second and third test systems (2TS and 3TS, respectively). In order to organize as much comprehensive analysis as possible, it seems necessary to compare the methods presented on the basis of the different dimensionalities of the problem. The test systems are considered as strongly connected ones, and the interest is related only to the test system generation profiles and generators to be repaired, the duration of these repairs, and the annual load profiles. Former data are presented for the three RTS systems in Tables 1-3. Table 4 contains the load profile data in p.u. in respect of the maximum annual load for the weeks of the year under consideration. The peak load values for the three test systems are 2850 MW, 5700 MW, and 8550 MW.   The available active generation power for the second test system is 6108 MW. Table 3. The generation profile of the third test system and quantity of generators subjected to major and midlife repair.

Quantity of Generators Subjected to Major Repair (8 Weeks)
Quantity of Generators Subjected to Midlife Repair (4 Weeks) 12 15 5  3  20  12  3  3  50  18  3  6  76  12  3  3  80  9  3  2  100  9  3  2  155  12  3  3  350  3  1  1  400  6  2  1 1 The available active generation power for the third test system is 9162 MW. In relation to the experimental section of the research, it should be noted that the objectives of the research were twofold. The first was to show the dependence of the maintenance planning results on the parameters of the heuristic methods considered, and the second was to compare them based on the results of the generator maintenance scheduling procedures completed using these methods.
In order to achieve the first aim, sequential maintenance planning procedures for 2TS were organized using the differential evolution and very fast simulated annealing approaches. The experimental dependence of the resultant objective function values collected laid the foundation for logical conclusions.
The second problem was solved through generator maintenance scheduling procedures for the three test cases with both user-defined parameters and the best parameters of the methods. Thus, conclusions could then be drawn directly from the resultant objective function values and their final deviation from the initial values.
Finally, in order to complete the computational experiments, all of the methods considered in the research were implemented in the C# software programming language using the Microsoft Visual Studio 2019 version 16.4.1 programming support environment. A graphical representation of collected results is provided, in turn, in Python language using the seaborn and matplotlib libraries described in much depth and represented in Ref. [37].

Influence of the Method Parameters on the Optimization Results
First of all, it is useful to consider, as far as possible, the state space for the best objective function values collected by applying the differential evolution and simulated annealing methods. In order to complete this task, the generator maintenance scheduling procedure was completed while the parameters of the methods were generated as pseudorandom uniformly distributed values. The results of these sequential simulations are represented in Figure 4a,b for the differential evolution and very fast simulated annealing methods, respectively.  The initial objective function value for 2TS calculated using the first-directed search method is 124.842 MW. Figure 4a shows that the use of the differential evolution method is connected with a significant probability of obtaining inappropriate results post-simulation. It can be seen that the optimization criteria values for the majority of the states are greater than the initial objective function value. This means that it is more profitable to use the directed search method. Nevertheless, it should be noted that for certain crossing-over coefficient values, in particular for the case shown in Figure  4a in the neighborhood of 0.2, there is an appreciable probability of obtaining better scheduling results.
In relation to the very fast simulated annealing method application results presented graphically in Figure 4b, one can see that there is a plateau of the state points on the initial objective function value level. This phenomenon could be interpreted as representing unsatisfactory planning results. The right-hand group of state points in the diagram corresponds to certain local minimum values of the state space, which has become an invisible barrier for the method on the way to achieving a global optimum value. Moreover, even for the points not included in either the massive plateau or the local minimum group, there is also the probability that the iteration process will stop at a local minimum.
For a more comprehensive analysis of influence of the parameters of the heuristic methods, the further subsections prove the analytical observation of the dependence of the resultant objective function values upon the individual parameters of the differential evolution and very fast simulated annealing methods.

Differential Evolution Method
In order to investigate how the population matrix size and the crossing-over coefficient independently affect the generator maintenance planning results, two experiments were conducted under the following conditions.

•
One hundred generator maintenance scheduling procedures were conducted for population sizes within the interval of [10; 100] with steps of 10, while the crossing-over coefficient = 0.4 remained constant.
• One hundred generator maintenance scheduling procedures were conducted for the crossingover coefficient within the interval of [0.2; 1.0] with steps of 0.1, while the population size = 50.0 remained constant.
The results of these consequent calculation experiments are represented in Figures 5 and 6 as sets of boxplots for every stage of the calculations.   Figure 5 shows that an increase in the population size results in a gradual decrease in the dispersion of the resultant objective function values, as well as the EENS. However, this positive effect is limited, because using a greater population size in the problem solution gives rise to the calculational efforts needed to complete one whole simulation. The last feature becomes more and more critical when taking into account the random character of the application results of the methods.
If one considers Figure 6, it may seem that there is no strong correlation between the crossingover coefficient and the generator maintenance scheduling results. In turn, an alert reader will notice that for low (e.g., = 0.2) and high (e.g., = 0.9) parameter values, the dispersion of the objective function is significantly greater. The mean EENS value is also different for all crossing-over coefficients, and the lowest ones are concentrated in the central part of the interval under investigation, except for the = 0.5 case.
The practical implication of the coefficient analysis by the differential evolution method is twofold. On the one hand, a greater size of population improves the quality of further simulation results but leads to a significant number of calculations. On the other hand, the optimal value of the crossing-over coefficient is highly likely to be located within the [0.3; 0.6] interval, but the latter conclusion depends on the length of the allowed maintenance periods during the time interval under consideration. Nevertheless, by following the differential evolution procedure, one may obtain inappropriate results even when the coefficients seem to be optimal. In particular, keeping in mind that the initial objective function value during the experiment was ( ̅ ) = 124.842 MW, almost all of the state points can be considered acceptable. Nevertheless, if the second-order directed search method is used as a reference method, then it appears that half of the points provide a worse solution for the task.
Finally, in addition to the assessment of the influence of the independent parameters on the maintenance schedule, it would be profitable to consider how the resultant objective function value changes because of the diversity of the generators to be repaired, as well as the reserves of the existing generation − . In order to represent the diversity of the generating equipment, the available power of all of the generators to be repaired was divided by 4, 2, 1, and 0.5, and simultaneously, the quantity of these generating units was multiplied by the same numbers. The value of the total generating power to be repaired in the power system remained the same, but the number of generating units to be scheduled was increased or decreased depending upon the provision of a new level of maintained generation profile diversity. Then, in order to consider the different existing generation reserves based on the 2TS data, the peak load for every week was increased in such a manner as to achieve a fraction of the power system generation reserve equal to 0.7, 0.8, 0.9, and 1.0 in relation to the initial value. Figure 7 illustrates the probability distributions of the resultant objective function values relative to the values calculated for the cases using the second-order directed search method. The first and obvious conclusion is that in all cases, the differential evolution method provides suboptimal schedules as a result of the optimization procedures. Nevertheless, one can see that the maintenance result becomes relatively worse when the diversity of the generating units decreases and the available active power of the units increases. In addition, Figure 7 underlines the fact that at first sight, for a number of cases, the differential evolution method provides inappropriate solutions. Nevertheless, deviation from the second-order directed search method results does not exceed 4% in any of the cases. Considering further application of the heuristic method for more a complicated problem statement, such a difference is not significant and rather underlines the closeness to the optimal solutions. Figure 7. The probability distributions of the resultant objective function values for a number of fractions of the power system generation reserve (70%, 80%, 90%, and 100% of the initial reserve) and the degree of diversity of the list of generators to be repaired (128 generators to be repaired (red), 64 generators to be repaired (blue), 32 generators to be repaired (orange), and 16 generators to be repaired (green)).

Very Fast Simulated Annealing Method
In relation to the very fast simulated annealing method, the effects of the temperature and temperature coefficient on the maintenance scheduling results should be considered. The following calculational procedures were completed in order to investigate these.

•
Five hundred generator maintenance scheduling procedures were provided for the pseudorandom temperature coefficient uniformly distributed within the interval of [0.10; 3.50], while the temperature = 1.75 remained constant.
• Five hundred generator maintenance scheduling procedures were provided for the pseudorandom temperature uniformly distributed within the interval of [0.01; 5.00], while the temperature coefficient = 1.50 remained constant. Figure 8 demonstrates dot-plot dependence of optimal solution objective function value on temperature coefficient. As it was in general shown by Figure 4a, influence of the last parameter of maintenance planning results is considerably complicated by configuration of extremum points in state space. In particular, right hand cluster corresponds to a number of local minimum solution points, as well as left hand maximum point group illustrates inability of the algorithm for such scenarios to find better solution then initial one. The behavior of the final objective function value when the temperature coefficient is constant is no less intricate, but the temperature itself varies in turn, which is illustrated in the dot plot in Figure 9. Nevertheless, even in this case, the points localized on the right side of the diagram can be seen to form a certain trend in contrast to the group of state points on the left, which again illustrates the unsuccessful simulation results. At this point, consideration of generation profile diversity and power system generation reserve effects on the simulation results is of interest in connection to the simulated annealing method. For the same test cases as were used previously for the investigation relating to the differential evolution method, a number of scheduling procedures were provided. The probability distributions for the test cases are represented in Figure 10. Figure 10. The probability distributions of the resultant objective function values for a number of fractions of the power system generation reserve (70%, 80%, 90%, and 100% of the initial reserve) and the degree of diversity of the list of generators to be repaired (64 generators to be repaired (blue), 32 generators to be repaired (orange), 16 generators to be repaired (green)).
As shown in Figure 10, the main difference between the application of the differential evolution and simulated annealing methods is that the latter does not guarantee finding a certain suboptimal solution in terms of maintenance planning. In particular, for the case of 16 generators and a 70% fraction, the scheduling results appear to be almost 10% worse than those found using the secondorder directed search method. However, despite the high inhomogeneity of the solution space, the cases of 32 and 64 generators present great planning results, being at least close to the suboptimal solution.
Finally, it can be concluded that the complexity of the solution space for the very fast simulated annealing method leads to difficulty in the initialization of the method parameters, as well as to a high probability of obtaining an inappropriate solution to the problem.

Comparison of the Generator Maintenance Scheduling Results
The comparison of the generator maintenance scheduling results is divided into two steps. First, a direct comparison is made of the calculational efforts needed and the resultant objective function values for the three test cases described above using the methods under discussion. Second, a graphical comparison is made of the possible deviation from the initial objective function values for the same cases using the heuristic methods.

Comparison of the Generator Maintenance Scheduling Results
Tables 5-7 represent the calculational effort estimations, as well as the resultant objective function values ( ̅ ). In order to assess the calculational efforts, the number of objective function calculations was counted. The justification for this approach to the estimation of former characteristics is that the perspective, especially in an optimization criteria assessment, is made more complicated when incorporating power flow constraints into the mathematical formulation. As a result, the number of objective function computations during the process provides an opportunity to analyze how efficient a certain method is. Let us consider the results presented. Tables 5-7 contain estimations of the implementation complexity of the calculation and the possibility to take into account complicated constraints, such as power flow constraints. In addition, the tables present the average calculation time for each method. Here, it should be noted that this parameter strongly depends upon the programming implementation of the calculation procedures. As a result, in order to avoid possible misunderstanding or under-or overestimation of the efficiency of the methods, the number of objective function calculations introduced above should be taken into account when considering the data.
Based on the results for 1TS represented in Table 5, one can come to the following important conclusions. The best method in this case is the second-order directed search method, which allows the lowest resultant objective function value to be obtained with insignificant calculational efforts in comparison to other ones. Nevertheless, what is more interesting here is that the very fast simulated annealing method provides results that are insufficiently better than the first-order directed search method. On this occasion, it seems that the initial maintenance plan variant used in the experiment was unsatisfactory, resulting in search process localization close within the neighborhood of the starting point. Additional calculational experiments related to similar planning procedures completed using the very fast simulated annealing method with the other parameters did not show any better results. Therefore, a question obviously remains as to how to appropriately initialize the simulated annealing method to ensure its effectiveness for all cases.
In contrast, the results in Table 6 illustrate a different situation-the very fast simulated annealing method provides one of the best results, which is also certainly true for the differential evolution method.
Similarity to the last test case is seen in Table 7, providing generator maintenance scheduling data for 3TS. In particular, the very fast simulated annealing method provides appropriate simulation results. However, the best optimal maintenance plan was calculated using the differential evolution method.
What was not expected was that the directed search method proved the effectiveness of the three methods, especially taking into account the comparatively insufficient computational efforts needed to use them. This demonstration allows application of these methods as a reference in further scientific and practical generator maintenance scheduling operations.
Despite it being reasonable to use the directed search methods, as proven by the results, their application is limited only to relatively simple maintenance planning problem formulation (i.e., without power flow constraints and only when generator scheduling is to be planned). Then, based only on the results presented in the tables above, it can be concluded that the differential evolution method is the more reliable heuristic approach to the problem in contrast to the simulated annealing method. The latter seems to be strongly affected not only by the dimensionality of the problem, but also by the initial starting point of the search process. At the same time, one should remember that the differential evolution method can hardly be referred to as a computationally effective approach, as the number of objective function calculations it requires is one of the greatest among all of the experiments. It would potentially be profitable to combine these methods in terms of one maintenance scheduling method in order to increase their applicability for the problem mentioned.

Dispersion of the Resultant Objective Function Values
Finally, it is interesting to estimate the effectiveness of the initialization procedures proposed in this paper, as well as the differences between the starting objective function values for the test cases and the resultant optimization criteria estimations. It seems to be useful to visualize these results as violin plots. Figure 11 presents the distribution of the resultant objective function values ( ̅ ) calculated using the differential evolution method for the test system in the following two cases.

•
User-defined coefficients, set arbitrarily for the differential evolution method, were applied.

•
The coefficients corresponding to the best values found in terms of the previous experiment were applied. Visualization of the experimental results presented in Figure 11 demonstrates the following features of the application of the differential evolution method. First, this technically simple method requires a certain approach to its initialization, since, for the greater dimensionality of the scheduling problem, it becomes a really tough assignment to set appropriate method parameters. For instance, user-defined coefficients for 3TS do not allow even one solution that is better than that of the starting point to be obtained. Nevertheless, as the second conclusion, even the optimal method parameters for the third case do not guarantee results that are better than the initial values.
A similar experimental procedure was followed for the probability distributions for using the very fast simulated annealing method: both user-defined and best solution coefficients were considered. Figure 12 represents violin plots for this investigation. It is reasonably simple to see that the resultant objective function values ( ̅ ) are very similar for 2TS and 3TS, verifying the effectiveness of the initialization procedure for the very fast simulated annealing method. Nevertheless, it seems that a question remains not about initialization but about the method itself, as obviously, in the case of 1TS, there is no generator maintenance scheduling result better than the initial one and, as a result, this method appears to be unreliable and not universal for maintenance planning procedures. Last but not least, Figure 13 illustrates dispersions of the optimization criteria estimations for both heuristic methods observed in the research. It demonstrates both the advantages and disadvantages of these methods. The differential evolution method always provides certain results that are different from starting point of the search, but an obvious shortcoming is represented by the significant dispersion of results and the lack of guarantee for obtaining results better than the initial values. At the same time, the very fast simulated annealing method provides results for which the objective function values are relatively closely localized, but it is likely to be strongly dependent upon the initialization and, as a result, is not reliable for the different dimensionalities of the problem under investigation. Figure 13. Boxplots for the resultant objective function value deviations from the initial points for the test cases under consideration calculated using the very fast simulated annealing (SA) and differential evolution (DF) methods.

Conclusions
In terms of this paper, deterministic and heuristic techniques for the generator maintenance scheduling procedure were compared. An initialization approach for the simulated annealing method was proposed, and possible perspectives for the development of a new methodology to the repair planning problem were considered. On the one hand, the directed search techniques provided advantages over the heuristic ones, as comparatively insignificant computational efforts were required to achieve appropriate maintenance scheduling results. However, it is hardly possible to adopt this approach for complicated mathematical formulations of the problem, or for considering different elements to be repaired, as well as power flow constraints. Consequently, application of the directed search method is very limited. On the other hand, the heuristic differential evolution and very fast simulated annealing methods provided relatively effective means to complete maintenance scheduling. However, both of them have certain disadvantages, demonstrated by this research, which again provides a significant barrier for their independent application in practice.
Nevertheless, despite the fact that a question remains about how to effectively consider and implement the relevant heuristic methods in industrial application, more efforts are needed to incorporate these approaches in terms of a single procedure. A possible way to implement this is to use the pattern of bee colony methodology. Here, the differential evolution method could play a role in the first initialization stage, while the very fast simulated annealing could be used for a local solution search. Such an approach may neutralize the existing issues connected to the individual unreliability of these methods and could provide a powerful tool for the maintenance scheduling of different types of equipment in power systems.