Swarm optimization approach to non-stationary physical field survey problem using a group of autonomous underwater vehicles Anton Tolstikhin, Igor Bychkov Matrosov Institute for System Dynamics and Control Theory of Siberian Branch of Russian Academy of Sciences (IDSTU SB RAS), Irkutsk, Russia E-mail: madstayler93@gmail.com Abstract. The paper considers the problem of searching for the source of a non-stationary physical field. We assume that the use of swarm algorithms may be applicable in this case. A hybrid of the Whale Optimization Algorithm and Grey Wolf Optimizer is proposed in this paper. The algorithm has several advantages over its origins: a more precise solution of the optimization problem for low-dimensional functions and a higher convergence rate of the first iterations. Two modifications were made to adapt the algorithm to the requirements of the problem. The proposed algorithm is used as a basis for a control strategy for a group of autonomous underwater vehicles. As a result, in the vast number of cases, the group can find the source within the given number of search iterations. 1. Introduction Although humanity has studied the Earth quite well, our knowledge of the underwater world is limited. Several negative factors cause this lack: the inability of a man to be in the underwater environment without special equipment, the high labor and resource intensity of research, and other physical limitations, particularly on the distance of communication and vision. Only relatively recently, with the invention of submersibles, the active phase of ocean research began. There are many different challenges in this area. One of them is the physical field survey problem (in other related works – adaptive sampling [1]). This problem includes many particular cases: mapping, search for given objects, search for areas with extreme values of the physical field, and many others. The solution of this problem has both fundamental importance - a general study of the underwater environment and improvement of our understanding of the processes occurring in it, and practical importance, for example, the search for oil fields or localization of the consequences of human-made disasters. Besides, the results of solving this task can be used as auxiliary information in other areas of marine research. Undoubtedly, this task can be solved with human resources. However, a more reasonable approach is to use autonomous underwater vehicles (AUV) or a coordinated group. Researchers often use a tack search algorithm or algorithms based on gradient descent to solve the trajectory physical field survey problem. The undoubted advantages of these approaches are simplicity Copyright c 2020 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). of implementation, computational efficiency, and acceptable work result in the vast majority of cases. Nevertheless, each of them has several disadvantages. Tack search implies a movement of AUV or their group along a predetermined zigzag trajectory and making measurements with a given frequency. As a result, the entire survey space is covered with a network of measurements. Processing of the obtained results allows AUVs to find an area with an extreme value of the physical field and to build field level lines. As can be seen from the description, this approach is not intelligent. Robots do not make any decisions during the mission and do not respond to environmental changes. Other shortcomings of the method include: • accuracy dependence on both the frequency of measurement and the width of the tack; • high execution time, depending on the size of the search space and the number of robots in the search group; • the need to store the values and the coordinates of all taken samples in memory; • the impossibility to apply the algorithm for non-stationary physical fields, because, by the moment of results analysis, the measurements may lose their relevance. These shortcomings are partly solved, using gradient descent. However, this algorithm implies only one agent, and the work of a group of robots can be considered as a multistart. In this case, interaction and information exchange within the group is minimal or absent. Finally, if the physical field has several local extrema, there is a high probability of hitting them. Population (in particular swarm) optimization algorithms may become a promising approach taking the above into account. Wildlife and inanimate nature processes underlying these algorithms often imply close interaction within a group to achieve a common goal. It allows one to transfer them to the AUV subject area with minor efforts. Despite the difficulty of proving the optimality of the solutions found, the swarm algorithms demonstrate high efficiency in solving practical problems. At the same time, it is often possible to dynamically influence the process of solving a problem, increasing or decreasing the number of iterations (as a consequence, measurements) necessary to obtain an acceptable result. This article suggests a swarm intelligence approach to a cooperative survey of non-stationary physical fields. Chapter 2 provides a brief overview of existing solutions. Chapter 3 describes the formulation of the problem with several accepted limitations and assumptions. Chapters 4 and 5 are devoted to a direct description of the developed algorithm as well as to some modifications vital for the non-stationary physical fields survey and the problem-solution process improvement. 2. Related works The physical field survey problem is a rather popular area of scientific research. In general, the works can be divided into two broad groups: the direct development of AUV for the physical field survey problem and the development of problem-solving methods. The first direction [2] includes the creation and improvement of sensors, devising robots for functioning in extreme conditions, and developing approaches to ensure communication within the group. Nevertheless, the second direction is of interest to this research. The most significant number of studies relate to two-dimensional physical fields in the water body and on the surface [3]. Significantly fewer works touch upon three-dimensional search space [4]. First of all, it is connected with the difficulty of robot vertical movement organization. Secondly, the three-dimensional space can be divided into two-dimensional layers, each of which may be examined independently. The survey of non-stationary physical fields is of particular interest; despite this, many studies are devoted to physical fields that do not change in time [5]. Depending on the specific task, the problem can be divided into searching for the source of a physical field [6; 7], searching for the area with the extremal value of the physical field (extremum) [2; 8], and mapping. Nevertheless, this division is conditional. For example, the source can simultaneously be an area with the extremum of the field in some particular cases. A parallel solution of more than one of the described tasks is the most preferable one. The most common approaches to the survey are the tack search [9] and algorithms based on gradient descent [10] as well as their combinations. It should be noted that swarm and other bio- inspired approaches are also used [6; 11] but in a much smaller number of cases. In particular, when it is known that the physical field is stationary or has a single source or extremum. Establishing sustainable communication between AUVs within a research group is quite a challenge, so most of the studies are aimed at using only one robot [3; 12]. Such an approach generates limitations on the size of the search space and may lead to a complete failure to solve the task. Therefore, despite the increasing requirements to the algorithms and control strategies, some researchers offer various centralized or decentralized approaches to solve the physical field survey problem using a homogeneous group of AUVs [13; 14]. Depending on nature, The following types of the physical fields [1] are naturally distinguished by their spatial characteristics, distribution features, and sensors required. • Chemical sources, for example, gas hydrates, dissolved salts, or oil spills. The task of finding the source or extremum value is most often considered [15]. • Physical sources presented by thermoclines or upwelling. The most significant interest, in this case, is mapping the phenomenon and building its two or three-dimensional structure [16]. • Biological sources, describing the distribution of certain biological species. It is challenging to choose the prevailing formulation of the problem because all of them are of equal interest to researchers [17]. 3. Problem formulation In this article, a two-dimensional case of the physical field survey problem is considered. In this case, the physical field is a kind of scalar mathematical field described by the function f (x, y, t) (Figure 1). The nature of the function is not known in advance for search agents (AUV). The function itself can be continuous, but to simplify calculations, we consider it to be defined on a grid. The value of each cell (x0 , y 0 ) ∈ F 2 (F ⊂ Z) is the maximum value of the objective function in the area bounded by cell borders and normalized as: f (∗) − fmin (∗) f 0 (x, y, t) = × 255 (1) fmax (∗) − fmin (∗) In this study, we discuss physical fields that are generated by several sources. Each source moves around the search space at speed, not exceeding the maximum speed of AUVs. The source trajectories are not known in advance. At the same time, the sources represent local extrema with the values fixed over time. The function, in this case, describes the distribution of physical field values at each moment. We do not identify any particular type of physical field at this stage of the study, considering it abstract. This assumption allows us to move away from considering specific features of studying a particular type of physical field, such as the sensors and sampling techniques used. Instead, we focus on solving the problem as a whole. Another assumption often found in works in this field is the absence of obstacles. Path- planning and collision avoidance are the objectives of the AUV navigation system, which lies at a lower control system level than is discussed in this study. For the same reason, we will consider that robots can accurately determine their position and are not affected by the environment, in particular currents. Figure 1. Example of the examined physical field representation The value of the physical field at the point (x, y) can be obtained only if any robot reaches these coordinates and takes samples (measurement). We assume that there are no restrictions on communication between AUVs, and each vehicle has information about all taken measurements (measurement array). In the current study, we consider the task of determining the source that produces the maximum value on the physical field. We assume that there is only one such source, which is simultaneously the extremum of the physical field. The problem is to determine the control strategy that leads by the end of the mission the group to such (x∗ , y ∗ ) ∈ F 2 at which the objective function takes the maximum value: f (x∗ , y ∗ , T ) = max f (x, y, T ) (2) 4. Algorithm development The article [18] discusses the solution of a similar problem - trajectory survey of a stationary physical field. The chosen Grey Wolf Optimizer demonstrated acceptable results. The framework of this research was proposed to create a hybrid algorithm to apply it to the current problem. The existing swarm optimization algorithms were reviewed and compared, to select two candidates for hybridization: the Whale Optimization Algorithm (WOA) [19] and the Gray Wolf Optimizer (GWO) [20]. Both algorithms have a similar structure and also showed good results on preliminary tests. It should be noted that there were earlier attempts to combine these algorithms [21; 22]. However, the proposed approach has its unique features and is not their full analog. 4.1. Algorithm description The proposed hybrid algorithm (WOA-GWO) has three behavioral patterns of agents inherited from the Whale Optimization Algorithm: exploration, exploitation, and the ”bubble-net attack” method. Switching between the patterns is independent for each agent at each search iteration and depends on two random variables A ∈ [−a; a] and P ∈ [0; 1]. In this case, a is the time function linearly changing from 2 to 0 as time a = 2 × (1 − ), (3) N where time is the number of the current iteration, N is the predefined maximum number of iterations. The exploitation pattern implies the specification of the found solution; in other words, attacking the prey. Activation occurs at P < 0, 5 and |A| < 1. The method for calculating the coordinates of the agent at the next iteration (target coordinates) corresponds to the one suggested in the Grey Wolf Optimizer. Three ”leaders” with the maximum value of the objective function are selected among all known measurements. F = f (x, y, t) (4) Next, the algorithm chooses temporary points located on lines passing through the current coordinates of the agent and the leaders (Figure 2). Thus, each temporary point is linked to one of the leaders. The location of these points is calculated separately for each dimension and depend on random variables Ah ∈ [−A; A] and C ∈ [0; 2]. Target coordinates in all patterns are defined as the arithmetic mean of temporary points. Xit = Xleadert − Ah × D, (5) D = |C × Xleadert − Xt |, (6) P3 i=1 Xit Xt+1 = , (7) 3 where Xt is the current coordinate of the agent, Xleadert is the leader’s coordinate, Xit is the coordinate of a temporary point. Figure 2. Exploitation example (blue area is the location of all possible target points) The exploration pattern describes the agents’ search for a potentially more promising solution, moving away from the current ones. Activation occurs at P < 0, 5 and |A| > 1. In contrast to exploitation, the target coordinates depend on three randomly selected measurements. The location of temporary points (Figure 3) depends on random variables As ∈ [−A; −1] ∪ [1; A] and C ∈ [0; 2] and is calculated separately for each dimension as follows Xit = Xleadert − As × D, (8) D = |C × Xrandomt − Xt |, (9) where Xt is the current coordinate of the agent, Xrandomt is the coordinate of a randomly chosen measurement among all previously made ones, Xit is the coordinate of a temporary point. Figure 3. Exploration example (blue area is the location of all possible target points) Finally, at P > 0, 5, the ”bubble-net attack” pattern is activated. It generates the behavior of agents similar to the way the whales are fed. The temporary points are selected on a spiral trajectory passing through the current position of the agent and the coordinates of each leader (Figure 4). These are selected as in the exploration pattern. Each coordinate of a temporary point depends on a random variable l ∈ [s; 1] and is calculated as follows Xit = Xleadert + D × eb∗l × cos 2πl, (10) D = |Xrandomt − Xt |, (11) l = (s − 1) × r + 1, (12) time s = −1 − , (13) N where Xt is the current coordinate of the agent, Xleadert is leader’s coordinate, r is a random value in the range of [0; 1], b is an adjustable parameter, adopted equal to 1, time and N are the number of the current iteration and the total number of iterations, respectively. 4.2. Testing Several tests aimed at constrained optimization of the given functions were carried out to analyze the proposed algorithm’s quality of work. Testing was performed on five functions (Table 1) that differ in the number of dimensions, the presence of local extrema, and belonging to the classes of convex and non-convex functions. Figure 4. Bubble-net attack method example Table 1. Test functions Function Dimensions Extremum Restrictions Beale 2 0.0 |xi | 6 4.5 Colville 4 0.0 |xi | 6 10.0 Goldstein 2 3.0 |xi | 6 2.0 Griewank 30 0.0 |xi | 6 600.0 Rastrigin 20 0.0 |xi | 6 5.12 In addition to hybrid WOA-GWO, the following swarm optimization algorithms participated in the study: • GWO — Grey wolf optimizer • GSA — Gravitation search algorithm • PSO — Particle swarm optimization • Bat — Bat algorithm • WOA — Whale optimization algorithm • ABC — Artificial Bee Colony • Firefly — Firefly algorithm • Gbest — Global best algorithm • MVO — Multiverse optimization • SOS — Symbiotic organism search There were three test series: the maximum number of iterations equal to 10, 100, and 1000. All series consisted of 50 independent runs of each algorithm with the same population size – 20 agents. Tests were conducted on Intel Pentium 4415U 2.3 GHz processor. Averaged results of each series are presented in tables 2–4. Table 2. Test results (1000 iterations) Beale Colville Goldstein Griewank Rastrigin ABC 0,000207724 0,570546723 5,70003257 0,005874714 0,004165207 Bat 0,158232589 88,54602738 11,10000009 230,6783114 79,05037751 Firefly 0,41311407 4230,039614 71,255498 57,94803945 44,97483681 GBest 0,09716027 8,604754176 12,46691855 3,833091684 5,143775522 GSA 0,001994162 1,546038046 3,000000612 639,2439958 18,80666382 GWO 0,038104855 1,785953011 3,000024301 0,004956 3,607712592 MVO 1,03163E-07 0,375888244 3,000003664 0,638474251 46,17326725 PSO 0 0,93648611 3 0,0317319 20,84438185 SOS 0,076289711 0,00664622 3 0 2,48452841 WOA 0,038103498 0,27585088 3,000408127 0,012729906 0 WOA-GWO 6,97403E-11 0,336746966 3,000000354 0,00651162 19,53432142 Table 3. Test results (100 iterations) Beale Colville Goldstein Griewank Rastrigin ABC 0,009523968 8,335194637 4,360152047 24,03029792 41,12014741 Bat 0,53042767 1168,542825 11,10000447 264,1918865 158,2655057 Firefly 1,276268951 6535,897502 47,74242201 443,6661571 146,766054 GBest 0,319842665 6,960217762 3,315806735 5,504292769 8,022998354 GSA 0,235132751 12,67997259 3,716502763 642,1539761 27,8563406 GWO 2,38421E-06 2,892877898 3,002142073 0,092457822 32,60208121 MVO 0,038130836 6,83868057 3,000348677 1,535644389 65,00815161 PSO 6,16152E-11 2,274432695 3 7,410641293 55,61478289 SOS 0,080755805 0,214794468 4,35000906 0 3,616272916 WOA 0,038110367 0,230467149 11,44800172 1,17611E-10 9,80638E-13 WOA-GWO 0,038107482 1,864989201 3,000098334 1,839233067 53,77797659 Table 4. Test results (10 iterations) Beale Colville Goldstein Griewank Rastrigin ABC 0,27544608 89,43936975 17,07573762 523,6649817 233,9436255 Bat 0,729765246 4096,957512 19,66647571 259,6968516 242,2274788 Firefly 2,464412748 7095,270527 89,35236493 594,9731594 230,884993 GBest 0,112107625 4,195715054 9,750612016 5,023009216 12,81697171 GSA 1,329928403 2415,862627 20,66805846 646,3561849 203,2148343 GWO 0,000775818 20,37462891 7,395616633 52,34030842 145,0470816 MVO 0,23339376 54,53605516 5,786769322 228,4235186 173,9305158 PSO 0,008820197 69,37418228 3,076063515 238,1131025 198,3255501 SOS 0,003617677 2,59503179 3,002812242 1,169603777 53,46858212 WOA 0,000280919 4,209769416 22,68463839 0,989905173 3,492248855 WOA-GWO 0,068984931 17,66576739 3,04665358 84,92464484 153,4045372 The following conclusions can be drawn from the results of the tests: (i) Regardless of the specified number of iterations, the proposed algorithm cannot cope with the optimization of Griewank and Rastrigin functions. It can be assumed that such a result will be achieved on all multidimensional functions. Since the physical field inspection task is considered only in two-dimensional and three-dimensional versions, this disadvantage can be ignored. (ii) On Beale and Goldstein functions, WOA-GWO shows the most accurate results among most of the considered algorithms, slightly inferior only to PSO and SOS, depending on the given number of iterations. On these tasks in the vast majority of runs, the proposed algorithm is superior to both of its ancestors. (iii) The algorithm shows average results on the Colville function. In this case, there is a high dispersion, so the results are inconsistent. On average, the proposed algorithm is superior to at least one ancestor and the PSO algorithm. (iv) In average, there is an increase in algorithm runtime, depending on the task (for 1000 iterations) by 2-4 times relative to WOA and 3-6 times relative to PSO. (v) WOA-GWO has a higher convergence rate in first iterations than all considered algorithms except PSO. However, then the speed is significantly reduced or stopped. (Figure 5). a) b) Figure 5. Best found values over time (a - for the Beale function, b - for the Colville function) In summary, we can conclude that a positive result of hybridization has been achieved. The existing shortcomings of the algorithm are either not significant for the problem or can be eliminated by further modifications. The high speed of obtaining the first acceptable approximation and the high stability of the work with the small size functions allows us to conclude that WOA-GWO is applicable to the physical field survey problem. 5. Modifications As it was mentioned earlier, a survey of a non-stationary physical field has one specific feature. Measurements may lose their relevance over time. That is why the algorithm needs a technique distinguishing acceptable and not acceptable measurements in necessary calculations. The proposed modifications are based on a three-dimensional Voronoi diagram, where the third dimension is time. Voronoi diagram is a division of space into subareas (cells), in which the distance from any point of the cell to a given point (site) is less than to any other site. In this case, the sites are the measurements made by agents during the problem solution. There is no need to build and store a three-dimensional diagram for the proposed approach; it is enough to have its two-dimensional cut at particular moments. If the area of any Voronoi cell on the given cut is equal to zero, the measurement that gave birth to it is thrown out of the list. Thus the ”forgetting” of old measurements, which in conditions of non-stationary physical field lose their relevance, is simulated. It was suggested to introduce a time factor allowing to adjust the speed of ”forgetting.” If its value is equal to 0, this process is disabled. Thus, all performed measurements are considered. If this parameter is increased, a smaller number of measurements is considered at the same time, up to the measurements obtained only at the current iteration. This process also indirectly allows one to regulate the labor intensity of the diagram cut calculation by reducing the number of measurements known to agents. The first detected barrier in solving the physical field survey problem, which is typical for both stationary and non-stationary formulations, is agent initialization. Usually, robots move to the examination area as a single group, so all agents should be situated on a small (in comparison with the whole search space) area at the first iteration. Such restriction leads to hitting the local extremum in the vicinity of the starting point. It was proposed to change the way of selecting ”leaders” in order to avoid the mentioned barrier. In the basic WOA-GWO algorithm, as in its parent algorithms, the leaders are three measurements with the maximum value of the function (4). A new function was proposed to replace it. F = k1 × f (x, y, t) × k2 × s(x, y, t), (14) where f (x, y, t) is the value of the objective function, s(x, y, t) is the value of the Voronoi cell area with the site at (x, y), k1 and k2 are weight coefficients (temporarily taken equal to 1). Function s(x, y, t) returns a value in range [0, 1] and represents the normalized area value of the corresponding cell. Thus, an increased interest in the inspection of cells with a large area is achieved. It provides a more uniform inspection of a given area without skipping zones that have received a bad initial approximation but may potentially contain an extremum. However, this approach negatively affects the final iterations of the algorithm. Therefore, the choice of the leader determination function has been synchronized with the basic patterns of WOA-GWO behavior: exploitation (survey) and operation (”attack”/refinement of the solution). In the ”exploration” behavior pattern, the leaders are selected according to (14), otherwise – according to (4). Such an approach allows solving both described problems. Another direction of work was to change the time function (3), according to which the choice of agent‘s behavior pattern is made. It has been noticed that during the final iteration, the agents examine a small area around a presumptive extremum (which size also depends on the time function) for a sufficiently large number of iterations. In the original algorithms, this allows specifying the value of the extremum with high precision. However, within the physical field survey problem (with the limitations and assumptions adopted), it makes no sense. It was proposed to change the time function time k a = 2 × (1 − ) , (15) N in such a way as to increase the exploration iterations percentage without changing the range of accepted values to eliminate this disadvantage. There is another indirect application of the Voronoi diagram for physical fields known to be stationary. With the time coefficient equal to 0, the result is a reasonably good approximation of the diagram to a given prototype (Figure 6). Thus, the mapping problem is solved in parallel. This additional information may be useful for further researches in this area. Figure 6. Comparison of the initial physical field represented by a height map and a Voronoi diagram 6. Conclusion The proposed control strategy is to build such a trajectory of each AUV on a dynamically generated set of points to gather the entire group in a small neighborhood of the current physical field source location. If there are several sources, the trajectory is constructed so that the source generating the maximum value of the physical field is found. A hybrid meta- heuristic optimization algorithm called WOA-GWO is responsible for the generation of key points in the trajectory. Based on a series of tests, it was found that on functions with a small number of unknowns WOA-GWO shows higher accuracy and speed than most currently popular algorithms. Since the considered problem is non-stationary, several modifications have been proposed, in particular, the mechanism of forgetting by robots of old measurements, based on the construction of the Voronoi diagram. Initial tests have shown that the underlying algorithms are unable to find a sufficiently good approximation to the extremum for non-stationary physical fields with any predicted probability. At the same time, proposed modifications allow finding a value close to the extremum with an average error equal to 5–7% and not exceeding 15% (at worst) in the overwhelming number of runs. The next step of this study is to check how the algorithm works on more complex physical fields, including both changing the coordinates of the source and its value. The reduction of the total survey time and assumptions made should also be considered. Finally, there may be a problem of losing a candidate solution on the final iterations of the search. Work is underway to make WOA-GWO more adaptive to avoid it. Acknowledgments The research was partially supported by the Russian Foundation of Basic Research, project no. 20-07-00397 A. The development of modifications for surveying non-stationary physical fields was also supported by the Russian Science Foundation, project no. 16-11-00053. References [1] Hwang J, Bose N and Fan S 2019 Applied Sciences 9 3145 [2] Zhang Y, McEwen R S, Ryan J P and Bellingham J G 2010 IEEE Journal of Oceanic Engineering 35 785–796 [3] Fahad M, Saul N, Guo Y and Bingham B 2015 2015 IEEE International Conference on Robotics and Automation (ICRA) pp 2654–2659 [4] Petillo S, Schmidt H, Lermusiaux P F J, Yoerger D R and Balasuriya A P 2015 OCEANS 2015 - Genova 1–10 [5] Medvedev A, Kiselev L and Tolstonogov A 2017 2017 IEEE Underwater Technology (UT) pp 1–6 [6] Pang S and Farrell J A 2006 IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics) 36 1068–1080 [7] Tian Y, Kang X, Li Y, Li W, Zhang A, Yu J and Li Y 2013 Sensors (Basel, Switzerland) 13 3776–3798 [8] Zhang Y, McEwen R, Ryan J, Bellingham J, Thomas H, Thompson C and Rienecker E 2011 J. Field Robotics 28 484–496 [9] Kiselev L V 2009 Autom. Remote Control 70 692–698 [10] Naeem W, Sutton R and Chudley J 2007 Journal of Navigation 60 173 – 190 [11] Farrell J A, Shuo Pang and Wei Li 2005 IEEE Journal of Oceanic Engineering 30 428–442 [12] Ai X, You K and Song S 2016 2016 14th International Conference on Control, Automation, Robotics and Vision (ICARCV) pp 1–6 [13] Khoshrou A, Aguiar A P and Pereira F 2015 [14] Chen B, Pandey P and Pompili D 2012 IFAC Proceedings Volumes 45 352 – 356 9th IFAC Conference on Manoeuvring and Control of Marine Craft [15] Baker E T, German C R and Elderfield H 2013 Hydrothermal Plumes Over Spreading- Center Axes: Global Distributions and Geological Inferences (American Geophysical Union (AGU)) pp 47–71 ISBN 9781118663998 [16] Zhang Y, Godin M A, Bellingham J G and Ryan J P 2012 IEEE Journal of Oceanic Engineering 37 338–347 [17] Das J, Rajany K, Frolovy S, Pyy F, Ryany J, Caronz D A and Sukhatme G S 2010 2010 IEEE International Conference on Robotics and Automation pp 4784–4790 [18] Tolstikhin A, Bakhvalov S, Dorofeev A and Bazhenov R 2019/05 7th Scientific Conference on Information Technologies for Intelligent Decision Making Support (ITIDS 2019) (Atlantis Press) pp 184–190 [19] Mirjalili S and Lewis A 2016 Advances in Engineering Software 95 51 – 67 [20] Mirjalili S, Mirjalili S M and Lewis A 2014 Advances in Engineering Software 69 46 – 61 [21] Mohammed H and Rashid T 2020 Neural Computing and Applications [22] Korashy A, Kamel S, Jurado F and Youssef A R 2019 Electric Power Components and Systems 47 644–658