Metaheuristic or evolutionary algorithms are a family of optimization algorithms inspired by nature, art, or social developments. Despite of their lack in convergence guarantee, these kind of algorithms have been popular in various engineering applications and scientific researches because of their simplicity and efficiency in finding global optimum solution for relatively low-dimension problems. Metaheuristics can be divided to single-solution approaches and multi-solution methods. The well-known examples of the single-solution approach are simulated annealing and tabu search 
algorithms in which one solution vector is evolved through strategic steps. On the other hand, in multi-solution(agent) methods a number of solutions are appointed and their relations are defined in an evolutionary or optimistic way such that they approach hopefully to a globally optimal solution. Hence, the main difference among the algorithms in the multi-solution approaches is in the kind of relations (motions) defined among the solutions.
Multi-solution approaches are recognized by two prior works on genetic algorithm (GA)
and particle swarm optimization (PSO). Other algorithms, for example, differential evolution , imperialist competition algorithm (ICA) , and teaching-learning-based optimization (TLBO)  are evolved versions of the mentioned two algorithms. They evolve to handle an smart balance between exploration and exploitation. The main differences are in the update rule and selection criteria. Specifically, in the swarm-based approaches the main difference is in the exploration strategy, since a motion toward a globally best solution(s) is a fixed part of these algorithms. The exploration item in PSO is accomplished by assigning a memory to the particles to save the best result of their search. In the ICA, the exploration was attended by appointing several imperialists selected from best solutions. Finally, in the TLBO, the exploration is a part of an educational system which is accomplished by interaction between pairs of the students chosen at random. Our proposed algorithm is another swarm-based optimization method in which the relation and movements among dull air particles of the tornado are mimicked. At the proposed algorithm, the exploration rule is designed inspired by spiral motion of the tornado. Our motivation behind modeling of the tornado motions was its powerful capability in the reduction of ambient temperature (up to C) by its mysterious air currents.
In fact, in formulation of a metaheuristic algorithm, there is a wide spectrum between random search and deterministic search where the intensities of exploration and exploitation of an new algorithm is determined. It gives an opportunity to design different algorithms to handle various tradeoffs in facing real-word problems. Two main drawbacks among the metaheuristic algorithms are run-time and parameter adjustability. At this work, both of them, with an emphasis on the former, is addressed. Two search directions is defined to boost the exploration and exploitation capability of the proposed algorithm. Each one is assigned to just one particle. That is despite of conventional methods, e.g. TLBO and PSO algorithms, which apply both main motions to each particle. Consequently, in the proposed algorithm, computational complexity is considerably reduced. Regardless of the step size parameters, the only parameter of the algorithm is the number of particles that should be assigned to each motion. The simulation results indicate that, in most of the problem, this parameter can be floated without significant degradation in the performance. This simple self-adjustment works, since the assignment is done at random.
Ii Tornado Morphology
Tornado is a kind of extreme weather. The strongest wind that could be formed in the earth is related to this phenomenon. When a layer of warm and moist air is located under the cool and dry air, due to the lower density of the warmer one it attempts to climb up to the top of the cold air. Conversely cold air descent to replace the rising air. If this process is occurred quickly, spiral airflow (like a funnel) is formed that is called Tornado . At the beginning of the Tornado the axis of the rotating airflow is horizontal. At the present of wind, the rising warm and moist air tilts the rotating air to vertical. A schematic model of Tornado is shown in Fig. 1 . During the Tornado, the warm air blows towards the tornado and rises through the spiral paths. This continues until combination of warm and cold air together and establishment of thermal equilibrium. Two main and specific air currents of the tornado, i.e. updraft and spiral motions, form the basis of our proposed algorithm. Moreover, atmospheric turbulence is an important part of these motions that has significant influence in the proposed simulated tornado optimization (STO) algorithm.
Iii Proposed Algorithm
Let the columns of matrix be consist of solution vectors ). Each solution includes decision variables. The solutions are interpreted as the position of air particles in a tornado. Their corresponding cost function determines the temperature of particle at position . The particles are divided into three types of coldest (), spiral (), and updraft particles (), where . Fig. 2 demonstrates this terminology and illustrates the idea behind the direction of motion for each particle. The coldest particle (), in top of all other particles, is best solution at each iteration. For a minimization problem that is a position vector with minimum cost function among the particles. Its replacements, after each iteration, simulates the storm motion (see Fig. 1, top of the left-hand picture). Other positions () are randomly sorted at the columns of matrix to participate as a spiral or updraft particle.
Updraft particles (the particles inside the funnel of a real tornado) move directly toward the coldest particle. The following relation formulates this kind of motion:
where indicates an entry-wise multiplication and is the turbulence vector, consist of step sizes for each decision variable. That is responsible for slight divergences from direct motion toward the coldest particle , in order to avoid from
a pure exploitation. The mentioned formulation is a common rule in the attraction-based algorithms . On the other side, as an special motion - designed to promote the exploration property of the proposed algorithm - a spiral particle at position moves toward the position of another spiral particle or the coldest one, i.e. . This index is determined according to the following criteria:
where denotes Euclidian norm. Hence, each particle moves toward another nearest-better spiral particle or the coldest one (see Fig. 2). Corresponding update rule for a spiral motion can be formulated as:
After each iteration, the coldest position is refreshed and spiral-updraft assignments are renewed, randomly.
Iii-a On the Turbulence Model
The Turbulence vector
imposes an uncertainty to both spiral and updraft motion. It is modeled by i.i.d. normal distribution with zero mean and unit variance. This distribution for both kind of updates leads to more generality and acceptable performance respect to the uniform distribution. In addition, that is more general in dealing with different problems respect to the realistic atmospheric turbulence models like log-normal distribution utilized in the communication channel modeling. Although, that is important to mention, according to our observations, log-normal distribution had shown some promisingly better performance than the normal distribution, at the case of unimodal problems. Nevertheless, for the sake of generality, our simulations are done using normal distribution.
Iii-B Tornado Diameter
As it is illustrated at the simulations (Fig. 8), the number of spiral particles can considerably influence the performance of the proposed algorithm in solving different problems. Its ratio respect to the total number of particles is called tornado diameter . Each problem has its own optimal value for tornado diameter. As a parameter-free version of the proposed algorithm, we suggest to float this parameter at each iteration to be selected uniformly among the feasible range of . Simulation results confirm the efficiency of the randomized tornado diameter. The finalized procedure of the proposed algorithm is summarized as follows:
Iii-C On the Convergence and Complexity
The movement toward best solution (coldest particle) in the updraft current is inherently a converging motion . However, that has a risk of falling into a local minimum. The spiral current boosts the exploration capability of the algorithm. Akin to the rings of a chain, the spiral particles are partially connected together and conduce an implicit motion toward the best solution. Hence, the spiral motion itself is a kind of converging attraction. Intuitively, the algorithm still converges after adding the spiral current. However, an analysis of convergence would have insight in choosing the step sizes. It is worth to mention that the randomized assignments of the motions at each iteration is not an essential part of the algorithm. Assignment according to the temperature can also be effective in some problems.
On the other hand, STO has low computational complexity, mainly because the motion of each particle is just in one direction at each iteration. That is despite of the algorithms like TLBO and PSO in which the final direction of movement for each particle at each iteration is determined by the computation of superposition of two or more direction vectors.
Iv Simulation Results
At this section, optimization of six different benchmark cost functions are investigated. These functions are shown in Fig. 3 and their formulations, along with domains, optimal solutions, and minimum costs are postponed to appendix. These functions were chosen from . Simulations are conducted through five experiments. Firstly, the trajectory of the particles in the proposed algorithm is demonstrated. Secondly, the effect of tornado diameter on the optimization performance is investigated. Thirdly, influence of the randomization of tornado diameter was considered in a quantitative comparison with a fixed diameter version. Next, fourth experiment includes some qualitative comparisons with with GA, PSO, and TLBO. Finally, the performance of the STO was evaluated on two higher dimension problems.
In all of the experiments the number of particles or population size was fixed on 40 for all of the algorithms. Distortion function is defined as , where is the real global minimum and is the result of optimization. One run of the optimization algorithms was regarded as a perfect optimization, if the distortion was less than . For Modified Rosenbrock and Rastrigin functions this definition was considered as a cost value less than 36 and , respectively. Parameters of the PSO were set according to the constriction coefficients with 
. In the GA, 0.8 of the population were selected for crossover and the rest for mutation. Crossover coefficient and mutation probability were tuned on 0.05 and 0.08, respectively. The TLBO is a parameter-free algorithm.
At the first experiment, trajectory of the air particles on artificial landscape of EggHolder function was demonstrated through Fig. 4 to Fig. 7. The number of spiral particles was fixed to be equal with that of updraft particles (including coldest one as a spiral particle). At the second experiment, influence of the parameter of tornado diameter was investigated. For this aim, the algorithm was implemented in whole range of the possible number of the spiral particles . Fig. 8 illustrates the effect of this number in the probability of perfect optimization of three test functions. The success rate was computed over 1000 trials at each values of the . As shown, EggHolder and Ripple functions are two examples of the functions that behave at extreme; the EggHolder function was efficiently optimized by high diameter tornados while the optimization of Ripple function was more efficient when diameter was small. In between there is Beale function which had low sensitivity to the tornado diameter. The proposed algorithm, had best performance in the optimization of Beale function when the ratio of the two spiral and updraft particles was in middle points.
As mentioned in the previous section, possibility of random assignment of the kind of motion for each particle position, allows to float the number of particles for each air current at each iteration. This randomization is not degrading and in most functions leads to a comparable performance with a tuned version of the STO. Table I compares the probability of perfect optimization at the case of randomized tornado diameter (indicated by STO) with the best result obtained by tuning of the number of particles assigned to spiral motion (indicated by STO). As shown, the performances are near (regardless of Modified Rosenbrock function). Roughly speaking, possibility of the randomization of parameter at each iteration indicates a promising self-adjustment capability of the proposed algorithm. In addition, the performance of other three algorithms are included at this table. The results are obtained over 1000 trials and all of the algorithms are stopped after 100 iterations. As inferred, the proposed algorithm has comparable quantitative performance in most of the problems and a better one at the case of Modified Rosenbrock function. In the rest of experiments, STO in the randomized diameter version was evaluated, except of Fig. 12.
|STO111best performance by tuning||0.95||0.94||0.99||0.54 ()||-|
parameter-free version (estimated)
Fig. 9 to Fig. 12 show a qualitative comparison of the algorithms on optimization of EggHolder, Ripple, Beale, and Modified Rosenbrock functions, respectively. The convergence curves were obtained by averaging over 50 independent runs. As shown in the first three figures, the proposed algorithm has comparable performance with the TLBO as it was expected from the quantitative comparisons. In addition, STO in both tuned and randomized diameter versions, has best performance in the optimization of Modified Rosenbrock (Fig. 12).
Fig. 13 shows the percentage of perfect optimizations (distortion less than ) for the Styblinski-Tang function respect to the dimension of problem. In two dimensions, this function consists of 3 local minimums and 1 global optimum (see Fig. 3). As shown in Fig. 13, all of the algorithms are completely successful in finding the global minimum at 2 dimensions. However, as the dimension increases, the PSO, GA, and TLBO were trapped in the local minimums such that their probability of perfect optimization dramatically decreases, while, the STO survives to explore new solutions in significantly higher dimensions. Maximum number of the allowed iterations for all of the algorithms was 5000. It was a limiting factor for the proposed algorithm, since it rarely fall into the local minimums of this problem and a better performance was possible by more iterations.
In the other experiment, convergence curve of the STO in optimization of Rastrigin function with 5 dimensions was evaluated in Fig. 14. This function has a global minimum in the origin, with the cost function vlaue of zero. The results indicate a comparable performance of the proposed algorithm with the TLBO, both in quality and quantity (see Table I). In order to have a sense of the low computational complexity of the proposed algorithm, a normalized run-time of the algorithms for last experiment was provided at Fig. 15.
Synchronized air currents in a tornado leads to a significant reduction in the ambient temperature. It is a desired paradigm in optimization. Simulated tornado optimization (STO) mimicked two main air currents of a tornado in movement toward colder positions (better solutions). Participation of each air particle in only one current, or equivalently movement in just one search direction, led to low computational complexity in the proposed algorithm. Hence, it would be useful in the applications with time constraint. A simple self-adjustment is possible by randomization of the number of particles appointed to each current (randomization of the tornado diameter). Numerical results indicated the efficiency of the proposed algorithm in facing with different problems. As a future direction of research, more improvements can be achieved by adaptive reduction in the range of variations of the diameter. Also, a more realistic step size mechanism or turbulence model would be of great interest.
Problems: (a) EggHolder Function:
(b) Ripple25 Function:
(c) Beale Function:
(d) Modified Rosenbrock Function:
(e) Styblinski-Tang Function with decision variables:
(f) Rastrigin Function with decision variables:
-  S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by simulated annealing,” Science, vol. 220, no. 4598, pp. 671–680, 1983.
-  F. Glover, “Tabu search — part I,” ORSA Journal on Computing, vol. 1, no. 3, pp. 190–206, 1989.
-  J. H. Holland, Adaptation in Natural and Artificial Systems. Cambridge, MA, USA: MIT Press, 1992.
-  J. Kennedy and R. Eberhart, “Particle swarm optimization,” in Neural Networks, 1995. Proceedings., IEEE International Conference on, vol. 4, pp. 1942–1948 vol.4, Nov 1995.
R. Storn and K. Price, “Differential evolution – a simple and efficient heuristic for global optimization over continuous spaces,”Journal of Global Optimization, vol. 11, no. 4, pp. 341–359, 1997.
E. Atashpaz-Gargari and C. Lucas, “Imperialist competitive algorithm: An
algorithm for optimization inspired by imperialistic competition,” in
IEEE Congress on Evolutionary Computation, pp. 4661–4667, Sept 2007.
-  R. Rao, V. Savsani, and D. Vakharia, “Teaching learning-based optimization: A novel method for constrained mechanical design optimization problems,” Computer-Aided Design, vol. 43, no. 3, pp. 303 – 315, 2011.
-  J. Snow, “Tornado,” in Encyclopaedia Britannica (T. E. of Encyclopaedia Britannica, ed.), 2015 ed., 2015.
-  E. P. Krider, “Thunderstorm,” in Encyclopaedia Britannica (T. E. of Encyclopaedia Britannica, ed.), winter 2016 ed., 2016.
-  X.-S. Yang, S. Deb, T. Hanne, and X. He, “Attraction and diffusion in nature-inspired optimization algorithms,” Neural Computing and Applications, pp. 1–8, 2015.
-  X. Tang, Z. Ghassemlooy, S. Rajbhandari, W. Popoola, and C. Lee, “Performance of the coherent optical binary polarization-shift-keying heterodyne system in free space optical communications using a lognormal atmospheric turbulence model,” in Mosharaka international conference on Communications, Propagation and Electronics, pp. 1–1, Mar 2010.
-  M. Jamil and X.-S. Yang, “A literature survey of benchmark functions for global optimisation problems,” Int. J. of Mathematical Modelling and Numerical Optimisation, vol. 4, no. 2, pp. 150–194, 2013.
-  M. Clerc and J. Kennedy, “The particle swarm - explosion, stability, and convergence in a multidimensional complex space,” IEEE Transactions on Evolutionary Computation, vol. 6, pp. 58–73, Feb 2002.