A hybrid algorithm based on state-adaptive slime mold model and fractional-order ant system for the travelling salesman problem

The ant colony optimization (ACO) is one efficient approach for solving the travelling salesman problem (TSP). Here, we propose a hybrid algorithm based on state-adaptive slime mold model and fractional-order ant system (SSMFAS) to address the TSP. The state-adaptive slime mold (SM) model with two targeted auxiliary strategies emphasizes some critical connections and balances the exploration and exploitation ability of SSMFAS. The consideration of fractional-order calculus in the ant system (AS) takes full advantage of the neighboring information. The pheromone update rule of AS is modified to dynamically integrate the flux information of SM. To understand the search behavior of the proposed algorithm, some mathematical proofs of convergence analysis are given. The experimental results validate the efficiency of the hybridization and demonstrate that the proposed algorithm has the competitive ability of finding the better solutions on TSP instances compared with some state-of-the-art algorithms.

Particularly, the pioneering ACO algorithm, ant system (AS) [4], was originally developed to solve the TSP [5]. On this account, the TSP is used mostly as a benchmark to test ACO variants. ACO is a swarm intelligence algorithm which simulates the foraging behavior of ant colony [11]. With the characteristics of positive feedback and information sharing among ants, it has become one effective method for solving combinatorial optimization problems. A considerable number of its variants have been designed to boost the performance, including the popular Max-Min ant system (MMAS) [12] and the rank-based ant system (ASrank) [13]. Meanwhile, the convergence proofs of AS and MMAS were presented to explain the rationale behind the ACO theoreti-cally [14]. Recently, two modified ACO algorithms based on fractional calculus have been proposed [15,16], which integrate the long-time memory property of fractional calculus into updating pheromones from a different perspective.
Generally, higher convergence speed and lower risk of trapping into local optima are the main targets of these extensions. In this context, how to make a dynamic balance between the exploration and exploitation of a swarm, has always been an important issue worthy of study and discussion. Considering the various advantages of different algorithms and search techniques, researchers have proposed hybrid intelligence algorithms by combining multiple methods. In [17], a parallel cooperative hybrid algorithm, PACO-3Opt, is developed to address the TSP, where the main idea is generating solutions by ACO in parallel, and then integrating the 3-Opt local search heuristic to modify these solutions. Another hybrid method called PSO-ACO-3Opt [18] utilizes PSO algorithm to adjust the parameters of ACO, and applies the 3-Opt technique to avoid premature convergence of ACO algorithm. These works show the effectiveness of the hybridization experimentally. However, there remains much room for improvement in accelerating convergence and jumping out of local optimal solutions. Moreover, they have made no attempt to provide the theoretical analysis of the convergence.
On the other hand, recently, there have been increased interests in slime mold (SM) model, due to its unique biological mechanism and the intelligent behavior of network design and path finding. An SM is an amoeba-like singlecelled organism with a tubular network structure and flowing protoplasm, which can adaptively adjust its body shape to construct a shortest path [19]. During foraging, the tube network, which contains nutrients and chemical signals that circulate throughout the organism, forms by connecting the growing point to food sources. These tubes are disassembled and reorganized according to external conditions, such as the quantity and position of nutrient sources, the terrains, and the light intensity distribution. In this way, the network structure of the organism is optimized to facilitate the efficient access to available nutrients [20]. For instance, when an SM is placed in a maze with two nutrient sources positioned at the entrance and exit, it modifies its body shape dynamically to absorb nutrients efficiently, and eventually connect the entrance and exit (a pair of food sources) with the shortest path [21]. Biological experiments have indicated that the simulation network constructed by SM model is very similar to the real Tokyo, Mexico, and Greece railway networks in terms of transportation performance, robustness, and cost [22][23][24]. Therefore, a biologically inspired mathematical model of SM has been developed with the positive feedback characteristic, which means the tube diameter increases with the increase of the protoplasm flux, and decreases otherwise [25]. Thereafter, the SM model has been successfully applied to solve combinatorial optimization, including the knapsack problem [26], traffic network optimization [27,28], load-shedding problem [29], supply chain network design [30,31], etc.
This model attracts our attention because apparently, the path optimization ability of SM model makes it a good candidate approach to solving the TSP. We noticed that a recent research has presented a slime mold-ant colony fusion algorithm (SMACFA) to deal with TSP, which uses SM model to determine the connections of some edges before ACO is implemented [32]. This approach is more like the initialization process before ACO, which takes the advantage of SM in shortest path construction. However, in this way, the variety of path selection is limited, and the advantages of SM model are not brought into full play.
The purpose of our work is to develop an effective approach to hybridize the AS algorithm with the SM model to bring the ability of SM model of constructing shortest paths into full play. To this end, a novel hybrid algorithm based on state-adaptive SM model and fractional-order AS (SSMFAS) is developed in this paper for addressing the TSP. Moreover, some convergence proofs are given to ensure the feasibility of the proposed hybrid algorithm theoretically. The main contributions of this paper are summed up as follows.
1. To make the SM model more suited for solving the TSP cooperating with the ant colony, we have developed a state-adaptive multi-entrance-exit SM model where each pair of nodes is regarded as the entrance/exit of an SM subsystem, equipped with two auxiliary strategies. The adaptive conductivity strategy changes the contraction rate smoothly over iterations, so that SM model shifts its state from exploration to exploitation to achieve a dynamic balance between local and global optimization. The maximum-minimum flux strategy is introduced to limit the flux through the tubes within an appropriate range matching the pheromone in the AS. 2. A hybrid algorithm based on state-adaptive SM model and fractional-order AS (SSMFAS) is proposed, where a fractional-order neighborhood probability is used in the node transition of ants to take the information of the neighborhood of candidate cities into account by fractional-order calculus, and a modified pheromone update rule is proposed which adds the protoplasm flux through each pipe as the pheromone on the corresponding edge dynamically to make the ant colony obtain information from the SM model. 3. The convergence properties of SSMFAS are analyzed in this paper. In simple terms, first, it has been proved that the probability of SSMFAS in finding an optimal path arbitrarily approaches to 1 given enough iterations. Second, we have proved that the pheromone trails of the optimal edges increase over iterations after the optimal solution is found, while the pheromone trails of other edges decrease. In addition, the lower bound of the probability of finding an optimal solution is given. 4. Abundant experiments have been conducted with promising results. First, ablation studies are implemented to validate the rationality of every module of the hybridization as well as the proposed two auxiliary strategies. Subsequently, comparisons with several state-of-the-art algorithms, including ACO-based hybrid algorithms and other heuristic approaches, are made to illustrate the effectiveness of SSMFAS.
The remainder of this paper is organized as follows. In "Background" section, we briefly introduce some algorithms and theories related to our work. The details of SSMFAS are presented in Methodology section. The theoretical convergence properties and the proofs are given in Theoretical convergence proofs section. Experiments section demonstrates the simulation results and the relevant analysis. Finally, Conclusions section summarizes the conclusion.

Background
In this section, we first give the definition of TSP and the notations used in the paper. Then we briefly introduce the framework of the AS algorithm, the mathematical model of SM and the background knowledge related to fractional calculus.

Problem definition and notations
In a TSP, assume the number of cities is N ∈ N + and the cities are labeled as i = 1, 2, · · · , N . The distance between city i and city j is given as L i j . The objective is to find the shortest route visiting each city once and returning to the starting point. Therefore, each solution to the problem is a sequence with length of N, i.e., C = (C 1 , C 2 , · · · , C N ). The cost function of the solution is For the sake of the reader's convenience, the mathematical notations used in the paper are listed in Table 1.

Ant system algorithm
Route construction and pheromone update are two critical steps of the AS algorithm. The process of the AS is summarized as follows: 1) Initialization: The population of ants is set to M ∈ N + . Every ant is placed at a random city as a start position. The initial value of pheromone on each edge is set to τ 0 > 0. 2) Route construction: Each ant m (m = 1, 2, . . . , M) in the current city i selects the next visiting city j depending on a random proportion rule of the transition probability, which is defined as where τ i j represents the pheromone concentration on edge (i, j), η i j is called the heuristic information which is calculated as η i j = 1 L i j , where L i j is the length of edge (i, j), α and β are two preset parameters used to balance the effects of heuristic information and the pheromone concentration, J m i represents the set of cities which have not been visited by ant m. Every ant keeps moving from one city to another guided by this probability and records the visited cities until it has visited all the cities.
3) Pheromone update: At time (t + 1), after each ant completing its tour, the pheromone on edge (i, j) is updated as where ρ (0 < ρ < 1) denotes the evaporation rate of the pheromone, τ m i j represents the volume of pheromone released by ant m on edge (i, j), which is calculated as where S m is the length of solution constructed by ant m. 4) Repeat the above steps until the termination condition is met, and then output the optimal solution.

Slime mold model
The SM model has a tube-liked network, which contains N nodes and N (N − 1)/2 edges connecting the nodes, indicating the cytoplasmic tubules between the nodes. Each edge has a weight, denoted by Q, which means the flux through the edge. Figure 1 illustrates the path construction process of an SM network. V 1 and V 2 are the entrance and exist nodes, respectively. V 3 , V 4 and V 5 are transition nodes. As shown in Fig. 1a, the classic SM model has only one growth point (entrance node), V 1 , and one food source (exist node), V 2 .   [25,33]. Therefore, the flux through tube (i, j), Q i j , is formulated as where L i j is the length of edge (i, j), p i and p j denote the pressures at nodes i and j, and D i j is used to describe the conductivity of edge (i, j), which is determined by where ω represents the viscosity of the fluid, and a i j denotes the radius of pipe (i, j). Following the Kirchhoff's law [34], the total flux should be conserved. Therefore, the outflow and inflow of each internal node must be equal, except for the start and end nodes. Let I 0 be the fixed total flux from the entrance to the exit in the entire network. This conservation law is modeled as Given an initial conductivity D 0 for every edge, and p exit = 0 as a fundamental pressure at the exit node, the linear equation system with sparse symmetric matrix (6) can be solved for the pressure value at each node, as well as the flux Q i j through each edge.
The slime mold adapts itself when foraging in a way that high-flow pipes are thickened, while low-flow pipes shrink and disappear. Since the length L i j is a given constant, the adaptive characteristic of the model can be expressed by the variation of conductivity D i j , modeled as where f |Q i j | is a monotonically increasing function with f (0) = 0, which describes the positive effect of increased flow on conductivity, γ > 0 is the contraction rate of tubes, which denotes the tubes contracts over time.
The discretized expression of (7) is expressed as where t represents a time interval which is typically considered as 1. Therefore, (8) is rewritten as From (4) and (9), it is noteworthy that there is a positive feedback cycle between the conductivity and the flux. That is, the conductivity of the edge, D i j , increases when a larger amount of flux, Q i j , passes through it, and this is conducive to the increase of flux further.
Under this mechanism, the network structure of the SM model in Fig. 1a changes to the one in Fig. 1b after a few iterations, where some tubes are disappeared (marked as the dotted lines), while some other tubes are thickened and stabilized (highlighted with bold black lines).

Fractional calculus
Fractional-order calculus is an extension of the integer calculus, where the order extends from an integer to a real number. As a mathematical tool, it has been investigated intensively and applied successfully in various areas [35]. The most popular used definitions of fractional calculus are Grünwald-Letnikov, Riemann-Liouville, Caputo, and Riesz [36][37][38][39]. In this work, we employ Grünwald-Letnikov fractional form because it can be discretized. For a differentiable function g(x) in the duration of [a, x], its Grünwald-Letnikov fractional-order differential definition of order ν (ν > 0) is defined as when ν = 1, it is obvious that (11) can be revised as is the first-order difference. By comparing (11) with (12), we observe that the fractional-order difference includes more history information, which is generally regarded as the property of long-term memory.

Methodology
In this section, we elaborate on the proposed SSMFAS, including the improved SM model and details of key procedures. Figure 2 shows the overall procedure of the proposed SSMFAS. After the initialization, a state adaptive multientrance-exit SM model is developed, where the positive feedback mechanism is used to update the flux and the conductivity. It consists of two strategies. The adaptive contraction strategy shifts the state of SM from exploration to exploitation. The maximum-minimum flux strategy limits the flux by the upper and lower bounds of AS algorithm, making it accessible for linking with AS algorithm. For the ant colony part, when the ants constructing solutions, the fractional-order neighborhood transition probabilistic rule which uses the neighboring information is introduced. After that, the flux of SM model is added to pheromone update rule as a part of adaptive pheromone. This loop continues until the stopping criterion is satisfied.

State-adaptive multi-entrance-exit SM model
We have introduced the classic single-entrance-exit SM model in Background section. However, it may not be suitable for solving the TSP directly since there is no fixed pair of entrance/exit in the TSP [40]. Therefore, we propose a stateadaptive multi-entrance-exit SM model where each pair of nodes is considered as a single-entrance-exit SM subsystem to address the TSP.
Let the number of the nodes in the improved SM model be N , which corresponds to the number of cities in AS. Accordingly, there are N (N − 1)/2 single-entrance-exit SM Fig. 2 The flowchart of the proposed SSMFAS subsystems since each pair of nodes is considered as the entrance/exit. Given the pressure at exit node as p exit = 0, the pressure at any other node in the nth (n = 1, 2, · · · , N (N − 1)/2) sub-network is derived from where I 0 denotes the total flux in the model, which is fixed during the optimization process. Then the flux q n i j (t) through the tube (i, j) in the nth subsystem at time t is computed as The total flux Q i j through the tube (i, j) is calculated by the sum of the flux in each subsystem, expressed as According to (9), the conductivity of tube (i, j) at time (t + 1) in the nth sub-network is updated by: where γ is the contraction rate of tubes and f (·) is a monotonically increasing function.
To make the improved model more applicable to the TSP, instead of using the two classical types of f (·) [25], we introduce a function as where μ is a parameter which controls the influence of q i j . Its value is optimized by orthogonal experiments in Experiments section.

Adaptive contraction strategy
Specifically, as the intention of integrating SM model with AS is to improve the global searching and convergence ability, we introduce an adaptive strategy, which changes the contraction rate of tubes, r , dynamically over iterations, to make the SM model maintain a dynamic balance between convergence speed and global searching ability. More specifically, at the beginning of the iteration when the global optimal has not been found, the algorithm should Fig. 3 The variation of contraction rate γ over iterations be more sensitive to various solutions, which is represented by a larger value of the contraction rate of tubes. Subsequently, the contraction rate of tubes should decrease with the number of iterations increasing, for the sake of accelerating the convergence of the algorithm.
To achieve this goal, we propose a smooth decreasing function within a range of (a, b), which is defined as where T is an inflection point when the algorithm shifts from the period of exploration to exploitation, and c is a coefficient to adjust the speed of the transition. The variation of the contraction rate γ (t) over iterations is plotted in Fig. 3, with a = 0.7, b = 0.2, and T = T max 3 . It shows that the value of γ is higher in the early search, which makes the conductivities of the tubes in SM sensitive to the flow, leading to the piping connections changing rapidly. In this way, the SM is in a more active state to look for various solutions at this stage. As the number of iterations increases, the value of the contraction rate decreases, which means that the conductivity changes slowly to form the optimal path gradually. In other words, the SM changes from active state to stable state over successive iterations.

Maximum-minimum flux strategy
Furthermore, there is another important point need to be consider when we try to combine the SM model with AS using the flux of SM model as a kind of pheromone. That is whether flux and pheromone, the two variables belonging to different models, keep similar orders of magnitude.
We notice that one simple but quite efficient mechanism of MMAS [12] is limiting the values of pheromone to the maximum and minimum ranges. Inspired by this, we consider an upper bound and a lower bound to limit the flux through each tube to a suitable range. More specifically, to better integrate with AS algorithm, the bounds should be related to the pheromone in AS, which is set to where τ min and τ max are the upper and lower bounds of the pheromone, respectively, δ > 0 is the weight coefficient, which is set based on some trial experiments in the algorithm.
The pseudo-code of the state adaptive multi-entrance-exit SM model is outlined in Algorithm 1.

Algorithm 1 State Adaptive Multi-entrance-exit SM Model
Input: The length matrix L; the maximum number of iterations T max ; the number of cites N ; Output: The total flux Q i j (t) at time t; 1: Initialize the conductivity matrix D n 0 ; the flow of the entire network I 0 ; 2: for t = 0 to T max do 3: Calculate p n (t) at each node based on (13); 5: Obtain q n i j (t) according to (14); 6: Update the conductivity D n i j using (16); 7: end for Calculate Q i j (t) using (15); 8: end for

Fractional-order neighborhood transition probability
In general ACO, an ant selects a movement using only the information between two nodes, which brings about a higher risk of trapping into local optima, as the best node for the next step is not always the best choice for the total path. Recently, a fractional-order ant colony algorithm (FACA) [16] has been proposed which calculates the probabilities of next few steps, and combines them with fractional-order calculus to get a more predictable probability. Similarly, another way of combined probability has been proposed in our previous work [15]. However, these algorithms need to calculate the transition probabilities of a few more nodes ahead of the current node. Therefore, the computational complexity is multiplied by the number of steps that the algorithm counts in.
To overcome this limitation, our goal is to take full account of the neighbor information without increasing too much computational cost. Here, we propose a fractionalorder neighborhood transition probability, which combines the moving probability of the candidate node with those of its neighboring nodes to form a new transition probability. Specifically, when an ant m (m = 1, 2, . . . , M) is at the ith node, the transition probability to the jth node at the tth iteration is defined as where ν is the fractional order, p m i j (t) is calculated by (1), (N − 1) is the number of nodes considered in the neighborhood of the node j, n k means the unvisited node which ranks k in descending order of Euclidean distance from node j, the denominator part is used to normalize this combination of probabilities.
For the candidate node j, the proposed fractional-order neighborhood transition probability considers the information in the neighborhood of j. This modification only needs the values of probabilities of the candidate cities for city i, which have been calculated already. Furthermore, it is logically reasonable that | | is a monotonically decreasing function of k, meaning the effect of its information decreases as the distance from the neighbor node to j increases.

Algorithm 2 SSMFAS Algorithm
Input: The length Matrix L, maximum number of iterations T max ; Output: The best solution; 1: Set iteration t = 0; 2: Initialize the parameters; 3: Place ants in starting cities randomly; 4: while t ≤ T max do 5: Each ant completes the tour using fraction-order neighborhood transition probability (20); 6: Local optimization and solutions sorting; 7: Use Algorithm 1 to get Q i j (t); 8: Update the pheromone on each edge according to (22) For the component of pheromone trails, similar to a rankbased version of the AS (AS rank ) [13], only the pheromone trails of some elite ants with ranks count. First, we sort the solutions constructed by M ants at the tth iteration in ascending order, written as Then, only the M (1 ≤ M ≤ M) best ants are used to update the pheromone trials and their effects are based on the ranks of their solutions. We use a nonlinear decreasing function with 1 as the largest weight instead of the smallest used in AS rank , since the flux in SM model is also added as the pheromone in this procedure.
The pheromone update policy of SSMFAS is given by and where M (1 ≤ M < M) is the number of chosen ants, λ denotes a control coefficient which is set to 0.1 in the experiments, σ (t) is a control parameter used to adjust the influence of protoplasm flux, which is defined as Similar to (18), σ (t) is a decreasing function within the range of (0, 1), which means the effect of SM model decreases with the number of iterations increasing. The reason for this is that the addition of the flux of SM model changes the pheromone distribution of the ant colony, and thus influences the edge selections of ants. At the beginning of the iteration, the algorithm needs to use the flux of SM model as part of information (pheromone) to guide ants to find some good edges. However, when the path search in SM model converges gradually, the algorithm reduces the effect of the flux to prevent excessive pheromone at some edges which would make it difficult for ants to find other solutions and thus fall into the local minimum. The adaptive hybridization of the two intelligent algorithms is conducive to better solutions. The general framework of SSMFAS is shown in Algorithm 2.

Theoretical convergence proofs
We analyze the convergence properties of SSMFAS in this section from several aspects, including the pheromone values and the probability of finding the optimal solution. Proposition 1 For SSMFAS, the pheromone value τ i j (t) on an arbitrary edge (i, j) at any time t satisfies that where max is the maximum value of pheromone (including flux) added to an arbitrary edge, determined as where S min is the theoretical shortest solution, Q max is the upper bound of flux.
Proof Based on (22), the maximum pheromone value on edge (i, j) at time 1 is derived as Accordingly, the maximum pheromone value on edge (i, j) at time 2 is obtained by In the same manner, the maximum pheromone value on edge (i, j) at iteration t is derived as As the pheromone evaporation rate, ρ, satisfies 0 < ρ < 1, based on the formula of summation for geometric sequence, we obtain that Thus, we have Thus, the proof is completed.

Proposition 2
On the premise that an optimal solution S min is found at the t * th iteration, for the pheromone value, τ i j (t), it holds that where edge (i, j) ∈ S min , τ max = (S min ) −1 is the maximum pheromone concentration added by ants.
Proof The pheromone updating expression (22) can be rewritten as wherẽ Thus, we have Since Q i j (t − 1) has an upper bound according to (19). It holds that Next, we only need to prove lim t→∞τi j (t − 1) = τ max ρ . On the premise that an optimal solution S min is found at the t * th iteration, the value of pheromone trail on each edge of the optimal path at the (t * + 1)th iteration is given as It can be concluded from (38) that Thus, by taking the limits of both sides of (39), we have This then completes the proof.
Theorem 1 Let P * (t) be the probability of SSMFAS finding an optimal route S min at least once in the first t steps. Supposing that t is sufficiently large, for an arbitrarily small positive number 0 < ε < 1, it holds that and lim t→∞ P * (t) = 1.
Proof Since the pheromone concentration on each edge is bounded by τ i j (t) ∈ [τ min , τ max ], the minimum value of the probability of choosing an optimal edge, ν p min has a lower bound given by where p min is reached theoretically when the pheromone values are τ min on the optimal edges, and τ max on all the other edges, which is written as Accordingly, the probability of finding an optimal solution S min is expressed as Thus, the maximum probability of not finding the optimal solution S min at the tth iteration is expressed as Therefore, the probability of SSMFAS searching out the optimal route S min at least once in the first t steps is given as On the premise that t is sufficiently large, we conclude that there exists an arbitrarily small ε such that Also, we obtain that lim t→∞ P * (t) = 1.
Therefore, the proof is completed.
Proof Under the worst circumstance that the pheromone concentration is τ i j (t * ) = τ min on optimal edge (i, j) ∈ S min , and τ kl (t * ) = τ max on non-optimal edge (k, l) / ∈ S min at time t * , τ i j (t * + t ) is expressed as and τ kl (t * + t ) is given as Let τ kl (t * + t ) be (1 − ρ) t τ max for our purpose here. Therefore, τ i j (t * + t ) > τ kl (t * + t ) holds when which is equal to Therefore, the proof is completed.

Proposition 3 After the optimal solution S min has been discovered, the value of pheromone τ kl (t) on any non-optimal edge (k, l) / ∈ S min decreases as the number of iterations increases. It holds that
Proof After the t * th iteration when the optimal path is found, the pheromone on the non-optimal edges evaporates without addition. The pheromone on edge (k, l) at the (t * + t )th iteration is derived as Thus, when t → ∞, τ kl (t) → τ min .
Proof Now suppose that t is the first time such that It is derived that This then implements the proof.

Corollary 1
Suppose that the first optimal solution S min is found at the t * th iteration. For the probability of ant m constructing the optimal route S min at the tth iteration, P * m (t) , it holds that where C is a constant and f (·) is a function of τ max and τ min .
where ν p * i j (t) is the transition probability for edge (i, j) ∈ S min after the first optimal solution S min is found, determined by where F is a constant defined as and p * i j (t) satisfies that Accordingly, the probability of ant m finding the optimal solution S min at the tth iteration satisfies Based on the binomial expansion theorem, we derive that Clearly, this then results in the consequence, (66), by setting The proofs of the above four propositions, two theorems and one corollary show the convergence properties of SSM-FAS. Proposition 1 reveals that the pheromone value on any edge has an upper bound. Proposition 2 proves that after the optimal path has been discovered, the pheromone concentrations on edges of the optimal route increase to the maximum with the increase of iterations. Theorem 1 demonstrates that the probability of finding an optimal path tends to 1 with the increase of iterations. Theorem 2 illustrates that finite times later after the optimal solution is discovered, the value of pheromone on any edge of the optimal route keeps larger than that on any other edge. Proposition 3 indicates that after the optimal path has been discovered, the value of pheromone on the edge that does not belong to the optimal route approaches the minimum value with the increasing number of iterations. Proposition 4 specifically declares that finite times later after the optimal solution is discovered, the value of pheromone on the edge that does not belong to the optimal route maintains the minimum. Corollary 1 exhibits that the probability of finding an optimal solution has a lower bound.

Experiments
This section presents extensive experimental results of the proposed SSMFAS, as well as some advanced peer methods. The experimental setup is given first. Then, the convergence curves of SSMFAS on different TSP instances are

Experimental setup
All experiments have been implemented in MATLAB 2019a environment in Window 10 on the PC equipped with an Intel Core i5 processor and 8GB of RAM.

Evaluation metrics
TSP instances obtained from TSPLIB of Heidelberg University 1 are used in the experiments. The number in each instance's name represents the number of cities. Each algorithm is tested 20 runs for one instance with random initialization, and each run includes 300 iterations. The minimal, maximal and average solutions acquired from 20 runs are recorded as the crucial indicators to evaluate the overall performance of algorithms.
The standard deviation (SD) is adopted in parameter tuning to evaluate the dispersion of the solutions, which is computed as where R (R = 20) is the number of runs, S i and S AVE are the solution in the ith run and the average solution of 20 runs, respectively. SD is considered as a measure of the robustness of algorithms since it is sensitive to data outliers. The smaller SD indicates that the algorithm is more robust.
1 http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/. The relative error (ER) is used to reflect the reliability of the results, which is defined as where S AVE and S MIN are the average and minimal solutions of 20 runs, respectively.

Parameter settings
For SSMFAS, α, β and ρ are set to 1, 5 and 0.1, respectively, referring to [41]. N and M are set to 8. The upper and lower (−ν) (m) | 1 S * and τ min = 2τ max /N , respectively, where S * is the approximate estimated shortest route length.
The fractional order, ν, the coefficient in conductivity update rule, μ, and the weight coefficient in maximumminimum flux strategy, δ, are searched to get appropriate values, since the performance of SSMFAS is sensitive to them. To determine an appropriate combination of values of these parameters, four alternative values are tested for each parameter, i.e., ν ∈ {0.7, 0.8, 0.9, 1.0}, μ ∈ {0.3, 0.55, 0.8, 1.05}, and δ ∈ {0.25, 0.50, 0.75, 1.00}. Accordingly, a four-level and three-factor orthogonal experimental design (OED) is utilized to choose a good combination of values of the parameters.
Let L a (b c ) represent the orthogonal array, where a is the number of tests performed, b and c are the level and number of parameters, respectively. Therefore, an orthogonal array L 16 (4 3 ) is procured, where 16 typical combinations are shown in Table 2.
One small-sized TSP instance, eil51, is used for this experiment. The average solutions of 16 combinations are listed in Table 3, where the best solutions are marked in bold font. We note that z 1 − z 4 represent the average of the results from 20 runs employing some certain parameters with specific values. For example, as δ = 1.00 is obtained by calculating (428.2 + 428.6 + 426.85 + 429.15)/4 = 428.2, the result 428.2 is placed in the corresponding position (z 1 , δ). From Table 3, it is clear that the best result is obtained when the parameter combination is δ = 0.25, μ = 0.8, and ν = 0.9.
Then, the optimal combination obtained by orthogonality,  Table 4 where the best solutions are emphasized in bold font. It is obvious that the combination of {0.25, 0.80, 0.9} wins on 4 of 5 instances. Therefore, we set δ = 0.25, μ = 0.8, and ν = 0.9 in later experiments.

Convergence property
The convergence curves of SSMFAS on 10 TSP instances are displayed in Fig. 4. It can be seen that the minimum solutions are reached within 150 iterations for all instances, which verifies the convergence of the algorithm experimentally.

Effectiveness of hybridization
An ablation experiment is designed to prove the effectiveness of the hybridization of AS, SM model and fractional-order difference components. Different versions of SSMFAS is developed for comparison. By setting the fractional order to ν = 1, SSMFAS degrades into an integer-order version, marked as SSMAS. By setting the parameter to σ (t) = 0, SSMFAS degrades into a version without SM model involved, marked as FAS. By setting both ν = 1 and σ (t) = 0, SSMFAS degrades into a version without the fractional calculus and SM model, marked as AS.
The average solutions as well as the best known solutions (BKS) are reported in Table 5, where the best solutions are marked in bold font. SSMFAS achieves the best solutions on all instances compared to other versions, which verifies the success of hybridization in SSMFAS. Furthermore, the results also show that both FAS and ASMAS perform better than AS, while there is no significant difference between them.
On the other hand, Table 6 presents the computation time of AS, FAS, SSMAS, and SSMFAS on different TSP instances. As can be seen from the table, insertions of the fractional calculation and SM model make the calculation time longer. This may be due to the long-term memory character of fractional-order difference is more computationally intensive, along with the fact that parallel computation has not been adopted in SM module. This limitation will be a major issue for future improvement.

Strategy evaluation
Two strategies proposed in the paper, adaptive conductivity strategy and maximum-minimum flux strategy, are evaluated separately to verify the effectiveness. The results are summarized in Table 7, where SSMFAS 0 represents the version without two strategies, SSMFAS 1 represents the version with the adaptive conductivity strategy but without the maximumminimum flux strategy, and SSMFAS 2 is the opposite version to SSMFAS 1 . The best results in each case are emphasized in bold font.
From the data in Table 7, it is apparent that SSMFAS achieves the best performances in terms of all evaluation metrics on the four TSP instances and has significant advan- tages than the others. It validates that the two specifically designed strategies improve the optimization capability and robustness of SSMFAS significantly. Second, the respective positive effects of the two strategies are evaluated by comparing SSMFAS 1 and SSMFAS 2 with SSMFAS 0 , respectively. SSMFAS 1 is superior to SSMFAS 0 in all TSP instances, which illustrates that the adaptive conductivity strategy has a positive impact on SSMFAS for solving TSPs. The comparison between SSMFAS 2 and SSMFAS 0 draws a similar conclusion.
Therefore, the employment of the two auxiliary strategies is beneficial to improve the performance of SSMFAS.

Comparison to other ACO-based hybrid approaches
To further investigate the efficiency of SSMFAS, several advanced ACO-based hybrid algorithms are used for compar- ison, including fractional-order ant colony algorithm (FACA) [16], heterogeneous adaptive ACO with 3-Opt local search (HAACO) [41], parallel cooperative hybrid ACO with 3-Opt local search [17], new hybrid PSO and ACO with 3-Opt [18], and two variants of MMAS algorithm [41] which are the standard MMAS incorporated 3-Opt and greedy, and only 3-Opt, respectively. We note that like the proposed SSMFAS, the six compared algorithms also include the 3-Opt local search procedure. The results of competitors are taken from the corresponding references. For a fair comparison, we set the values of the common parameters in our algorithm to be the same as those in [41]. Table 8 presents the common parameter settings of the algorithms. Other specific parameters are set following the original setting.

Performance evaluation on small-scale problems
First, SSMFAS are evaluated on 10 Small-scale TSP instances. Tables 9 and 10 present the minimal and average solutions of SSMFAS and other algorithms, respectively. For the purpose of analyzing the performance of the algorithms, the ranks of the solutions are computed and listed in the parentheses in those tables. We note that if any values are tied, we compute their average rank. Furthermore, the average ranks of each algorithm on the 10 instances are also computed and reported in the last lines of tables. The best solutions in each case are written in bold and highlighted in gray, and the second-best results are also written in bold font. As shown in Table 9, the minimal solutions of SSMFAS on all TSP instances are the shortest over seven algorithms. As a result, SSMFAS obtains the smallest average rank of 2.65, and FACA with an average rank of 3.25 comes next. And, perhaps more tellingly, in terms of the average solutions, Table 10 indicates that SSMFAS outperforms other algorithms with an average rank of 2.20, followed by PACO-3Opt with 3.15.
Now we check whether the solutions constructed by these algorithms are significantly different. A common statistical method, Analysis of Variance (ANOVA) [42], is utilized to test if there is significant differences in these algorithms.
The null-hypothesis is tested and the generated p-value is 9.95e − 8, which is less than the significance level of 0.1. Thus, the null-hypothesis is rejected. The Wilcoxon signed-ranks test [43] is then used to verify the significant improvement of the proposed SSMFAS in pairs. The results are reported in Table 11. As can be seen, the generated pvalues are all less than the significance level of 0.1 except for FACA. This may due to the fact that the two algorithms are tied on two instances. More specifically, Fig. 5 shows the shape of distribution of ranks in terms of average solution. The distribution of SSMFAS is more concentrated with the smallest median, which clarifies the effectiveness and robustness of the proposed SSMFAS.

Performance evaluation on larger-scale problems
Furthermore, to illustrate the ability of the proposed SSM-FAS to handle large-scale TSPs, 11 lager-scale TSP instances from TSPLIB whose numbers of cities are between 400 and 800 are used in this subsection.
Tables 12 and 13 present the minimal and average solutions of SSMFAS and its competitors, respectively. The best solutions in each case is highlighted in bold. It can be seen from the data in Table 12 that SSMFAS wins on 8 of 11 instances. From the average solutions in Table 13 we can see   It is worth mentioning that the superiority of the proposed SSMFAS becomes more obvious compared with the competitors when the problem size becomes larger. This may illustrate the proposed SSMFAS has good robustness and adaptability to large-scale problems.

Comparison to other heuristic approaches
In addition to ACO-based hybrid methods, some advanced heuristic methods are involved in the comparison, including a discrete stochastic population-based optimization algorithm (DJAYA) [44], an approach of improvement heuristics based on 2-opt operators by deep reinforcement learning [45], and an open source software for combinatorial optimization, Google OR-Tools [46], which includes 2-opt and LKH   (Lin-Kernighan-Helsgaun) [47] as improvement heuristics [45,48]. The results are reported on Table 14, where the best result in each instance is highlighted in bold. We note that the results of SSMFAS are the average solutions of 20 runs, and the results of competitors are taken directly from the corresponding references.
As can be seen from Table 14, SSMFAS outperforms other approaches on 9 of 10 instances. It further verifies the superiority of the proposed SSMFAS in addressing TSPs.

Conclusion
In this paper, we proposed a novel algorithm called SSMFAS combining a state adaptive multi-entrance-exit SM model and a fractional-order based AS algorithm for solving TSPs.
The adaptive conductivity strategy is developed to shift the state of SM in different periods. The maximum-minimum flux strategy limits the upper and lower bounds of the flux to match with the concentration of pheromone of AS algorithm. Fractional-order neighborhood transition probability which uses the neighbor information is introduced for path construction by ants to improve the performance. Varying degrees of flux in SM model are added as pheromone in the pheromone update process to provide more information.
The convergence properties of the SSMFAS have been verified. A multitude of experiments are carried out on TSP instances. First, the effect of the hybridization of AS, SM, and fractional-order calculus is verified. Then, the effectiveness of two auxiliary strategies is demonstrated. At last, comparisons with several state-of-the-art algorithms illustrate the competitive performance of SSMFAS.
In the future work, we plan to optimize the combination of SM model and AS algorithm, and improve the computation cost through some parallel strategies, to address some dynamic multi-objective TSPs. More broadly, further research could also be conducted to explore the potential of combining other path planning methods and metaheuristic algorithms.