A Bio-Inspired Method for the Constrained Shortest Path Problem

The constrained shortest path (CSP) problem has been widely used in transportation optimization, crew scheduling, network routing and so on. It is an open issue since it is a NP-hard problem. In this paper, we propose an innovative method which is based on the internal mechanism of the adaptive amoeba algorithm. The proposed method is divided into two parts. In the first part, we employ the original amoeba algorithm to solve the shortest path problem in directed networks. In the second part, we combine the Physarum algorithm with a bio-inspired rule to deal with the CSP. Finally, by comparing the results with other method using an examples in DCLC problem, we demonstrate the accuracy of the proposed method.


Introduction
The shortest path problem (SPP) is one of the most classic problems, which is widely used in many fields, like network optimization [1,2], road navigation [3,4], and so on [5][6][7]. The constrained shortest path (CSP) problem aims at finding the shortest path from a predetermined source node to a predetermined destination node in a directed network under the constraint of being less than or equal to a given upper limit. Such constraints on the classical shortest path problem generally result in the problem becoming NP-hard. CSP is often considered as a variant of SPP and is frequently put into practice, like tail assignment problem in aircraft scheduling [8][9][10], the quality of service routing in communications networks [11,12], and so forth [13].
Strategies of solving the CSP can be classified into three main categories: Lagrangian relaxation methods, dynamic programming (DP), and path ranking methods. The first method of CSP relies on the solution to a Lagrangian dual problem and uses different approaches to close the duality gap. Handler and Zang [14] proposed a method to solve the Lagrangian dual problem. They closed the duality gap by using the shortest paths algorithm, terminating with the first path satisfying the constraint. They also proposed a method which searched in a direction. The method was based on combination of the original objective and the constraint. Carlyle et al. [15] presented a new algorithm for CSP. The improvement of this method is that the generated Lagrangian function was optimized and the optimality gap was closed by enumerating near-shortest paths. Recently, Zhang et al. [16] used Lagrangian relaxation combined with an adaptive amoeba algorithm to solve the CSP. In fact, if the constraint was relaxed, it made the CSP as the classical SPP, which could be solved repeatedly by plentiful computational efforts. What is more, the efficiency of this approach relied on the effectiveness of the underlying unconstrained shortest path algorithms.
A number of previous works were based on DP. All of these methods were based on label-setting or labelcorrecting algorithm. Starting from the the work of Joksch [17], the approaches of node labeling were deeply investigated. Dumitrescu and Boland [18] extended this idea and solved this problem by using preprocessing techniques. They [18] showed how the performance of the label-setting method could be further improved by making use of all Lagrange multiplier information collected in a Lagrangian relaxation's first step. Recently, Likhyani and Bedathur [19] developed SkIt index structure, which supported a wide range of label 2 The Scientific World Journal constraints on paths, and returned an estimation of the shortest path that satisfied the constraints. However, all of these methods mentioned above had a disadvantage: they might fail to well-estimate very large networks due to the "curse of dimensionality. " As a result, it could not apply to real world with a large size of nodes.
A path ranking method was presented by Handler and Zang [14] by shortest paths algorithm. Santos et al. [20] extended this idea by improving the way of sorting the shortest paths. Its direction was based on the relative tightness of the constraint. Jia and Varaiya [21] presented two new methods based on the -shortest-path approach. These heuristic methods, where one is centralized and the other is distributed, were both polynomial. In numerical experiments the proposed algorithms almost found the optimal paths. The main drawback of this method is that large value for is required before a feasible solution of CSP is found.
Inspired by the intelligence of the amoeba, we combine the Physarum algorithm with a bioinspired rule to deal with the CSP. The main idea behind the approach is to employ Physarum algorithm to find the shortest path. Then a check procedure is processed to ensure it satisfies the constraint. Once it does not meet the constraint, a penalty function is conducted. This procedure will go on executing until the path satisfying the constraint is found.
The paper is organized as follows. In Section 2, the basic theories, including CSP and the Physarum model, are introduced. In Section 3, the bioinspired method which is based on amoeba algorithm for CSP is proposed. In Section 4, an example is used to illustrate the proposed method. In comparison with other methods, the validity of the proposed method is demonstrated. In Section 5, the conclusion is given.

Basic Theories
In this section, we will introduce some basic theories, including the constrained shortest path problem and Physarum polycephalum model. Each edge has two nonnegative weights and . and represent the generalized cost and the constrained variable, respectively. The CSP can be expressed by the following mathematical formulation: subject to where is binary variable, which is defined as follows: The parameter indicates the maximum value allowed for the sum of . Equations (1), (2), and (4) define the shortest problem. The constraint (3) causes the CSP to be an NP-hard problem.

Physarum polycephalum Model. Physarum polycephalum
is a single-celled amoeboid organism that can form a dynamic tubular network to link the discovered food sources during foraging. The physiological mechanism behind the phenomenon is that tubes thicken in a given direction when the flux through it persists in that direction for a certain time. By this mechanism, a mathematical model for path finding has been constructed to describe the adaptive dynamics of tubular network. In what follows, we will give a brief introduction to this path finding model [56].
Suppose the shape of the network formed by the Physarum is represented by a graph; the edge and the junction in the graph can be treated as a plasmodial tube and the node, respectively. Two special nodes labeled as 1 and 2 act as the source node and sink node, respectively. The other nodes are labeled as 3 , 4 , 5 , 6 , and so forth. The edge between nodes and is expressed as . The parameter denotes the flux through tube from node to . Assuming the flow along the tube is an approximate Poiseuille flow, the flux can be expressed as where is the pressure at the node , is the conductivity of the tube , and is its length. By considering that the inflow and outflow must be balanced, we have The Scientific World Journal 3 As for the source node 1 and the sink node 2 the following two equations hold: where 0 is the flux flowing from the source node, which is set as a constant in the model.
Then the network Poisson equation for the pressure can be derived from (6)-(8) as follows: By setting 2 = 0 as a basic pressure level, all can be determined by solving (9) and can also be obtained according to (6).
With the mechanism between the flux and tube thickness (described as the conductivity of the tube), the adaptation equation for conductivity is constructed as follows: where is a decay rate of the tube. It can be obtained from the equation that the conductivity will finally vanish if there is no flux along the edge, while it is enhanced by the flux. The is monotonically increasing continuous function satisfying (0) = 0.

Proposed Method
The original amoeba algorithm can only find the shortest path in undirected networks. In order to solve the CSP in directed networks, we have to solve two problems. The first one is how to find the shortest path in directed networks. The second one is how to efficiently solve CSP by modified amoeba algorithm.

Modified Amoeba Model for the Shortest Path Problem in
Directed Networks. Let = ( , , ) be a directed network, where represents the set of nodes, represents the set of directed edges, and represents the weight set of . We suppose that there is only one direction between node and node in the network . Assume node and node are source node and sink node, respectively. The shortest path problem in the network can be defined as how to find a path from node to node , which only consists of directed edges of , with the minimum sum of weights on the edges.
In the original amoeba model, every arc is bidirectional. It means that the flux can flow from node to node . They can also flow in the opposite direction. It is no different for original amoeba model, while the direction should be strictly distinguished in the directed networks. In particular, there is only one direction between two nodes in network . It means that if the flux can flow from node to node , they cannot flow in the opposite direction. In order to take the factor of direction into account, each edge is regarded as two tubes with opposite direction and equal weight. Here, we assume the weight represents the length between two nodes represented by . If there is a link between nodes and , then there are two tubes with equal length and , respectively, denoting two opposite directions from node to node and from node to node . The initial value of conductivity along with each tube is used to imply whether the tube is accessible or not. If arbitrary tube satisfies ∈ , then the corresponding initial value of conductivity is assigned an arbitrary value among [0. 5,1]. Besides, is assigned as 0. The matrix implies not only the conductivity but also the direction of each edge. By this way, the factor of direction is considered. In order to make original amoeba model that can adapt the directed networks, (9) is improved as follows: With the purpose of keeping the validity of conductivity, the flux along with each tube needs to be adjusted. Equation (10) is modified as follows: Amoeba model has a positive feedback mechanism: the higher the pressure difference is, the more the flux will be. The more the flux is, the higher the pressure difference will be. When the iteration continues, the shortest path gradually appears through this mechanism. In directed networks, the positive feedback mechanism also exists. The method for finding the shortest path in directed network is described as follows.
(1) Initialize values of each parameter. According to the weight of the edge in , the length of is initialized. According to the direction of each edge, the conductivity is initialized. Here, the initial value is assigned equally as 0.5. If there is an edge from node to node , = 0.5. Otherwise, = 0. The pressure of the node is represented by . The basic information of this directed network is initialized.
(2) According to current conductivity, the pressure of each node can be calculated using (11). (3) According to (6), we can obtain the flux of each edge. (4) According to (12), we can calculate the conductivity of the next moment. (5) If the conductivity of each edge is no longer changed, the termination criterion is met. The directed shortest path consists of the edges with conductivity approximately equaling 1. Otherwise, go back to Step 2 and repeat until convergence.

Bioinspired Method.
In this section, we will introduce some nomenclatures and phenomena about amoeba. Then we will introduce the proposed method.

Nomenclatures and Phenomena
Segment: the tubular pseudopodia of amoeba linking node to node .
Growing segment: the segments with continuous thickening.
Potential segment: the growing segments with growing times exceeding given threshold.
Threshold: the critical value of judging whether growing segment is potential segment.
Important segment: the potential segments which can make up an unbroken path. Important path: the path which is made up by potential segments.
Conductivity: indicating the thickness of the segment.
As the iteration continues, some segments become thicker while others degenerate. Simultaneously, the corresponding conductivity tends to increase or decrease. When the conductivity value in the current iteration is greater than the previous one, it is called thickening segment. When a segment is thickening, it means that the segment will continue to thicken. The segment becomes a growing segment. If the growing times of the growing segment exceed the threshold, it becomes a potential segment. If the potential segment belongs to the important path, it becomes an important segment. Through many trials, we find that the segment which belongs to the shortest paths continuously thickens from the beginning to the end. At the beginning, growing segments may be dispersed in the network. Then, potential segments and important segments occur. At the end, there are only important segments left, which make up a shortest path. In other words, we think some of those important segments are the final part of the shortest path and the important path is the shortest path on the great probability. Figure 1, an example employed in [16] is used to illustrate the general flow of the proposed method. Each arc has two attributes: the length and the cost . The specific values can be obtained according  50 40 to Table 1. For instance, the numbers (120 and 60) along the arc between node 1 and node 2 mean the length along the arc is 120, and the cost along the arc is 60. We want to find a length shortest path from source node 1 to sink node 20 with the cost no more than 200. Inspired by the important segment, we propose a bioinspired method based on the internal mechanism of the amoeba model. Steps are introduced in detail as follows.

Proposed Method. As shown in
(1) On the basis of Step 1 in Section 3.1, we need to do some extra operations as follows.
According to the information of the given CSP, we initialize the constraint . is a parameter which is used to judge whether the growing segment is the potential segment or not. is the divisor of the punishment. When the important path does not satisfy the constraint conditions, is used to decrease the conductivity of those important segments so that the amoeba can converge to the path satisfying the constraint. Besides, we set a variable flag representing the end mark and assign 0, which means the end mark does not satisfy the termination criterion.
(2) Step 2 to Step 5 are the same as in Section 3.1.
The Scientific World Journal 5 (3) According to the current and the previous conductivity, the growing segments can be found and marked. The growing times of each growing segment are also saved.
(4) According to the saved times and the threshold , we can judge whether the growing segment exceeds the threshold. If it exceeds it, it is the potential segment. Save these potential segments.
(5) Judge whether these potential segments can make up a complete path from the source node to the sink node. If they can, we mark these important segments and then go to Step 6. Otherwise, go to Step 7.
(6) According to these important segments, the important path can be obtained. Judge whether the important path satisfies the constraint . If it can, the original data such as , is kept and makes the flag equal 1, which means the flag satisfies the termination criterion. In other words, the constrained shortest path is found. Then go to Step 7. Otherwise, do the following. The conductivities of each important segment (the segment who belongs to this important path) are decreased. The updated conductivity is equal to the original value divided by the divisor of the punishment . The marks including times and important segments are initialized.
(7) Check whether the termination criterion is met or not. The termination criterion is that the value of flag is 1. If it is satisfied, tubes with maximal conductivity make up the constrained shortest path. Otherwise, go to Step 2 and repeat until termination is satisfied.
More detailed procedures are showed in Algorithm 1. The basic idea behind this method is that we take full advantage of Physarum algorithm to find the path. After finding the path successfully, the check procedure is conducted to insure the path satisfies the constraint. If it does not meet the constraint, the penalty mechanism will be processed. By this way, the successive path will fade out. The program will continue to execute until the path satisfying the constraint is found. As can be seen in Algorithm 1, it is necessary for us to solve the linear equations shown in (9). The time complexity of solving the linear equations is ( 3 ), where the | | is the number of equations. In the original amoeba model, the number of equations and nodes in the network are equal. Hence, the time complexity of the bioinspired method is ( 3 ), where the | | is the number of nodes in the network. According to Figure 1, we let be 1, 20, and 200. Parameters and are experimental values. Here, we let them be 3 and 10, respectively. As shown in Figure 2, the solution 1-5-9-16-20 to this problem is obtained easily. The length of the path is 340 and the cost of the path is 200.
It can be noted that the solution obtained by the proposed method is the same as the result presented by Zhang et al. [16]. Both of them can find the optimal solution. Therefore, the accuracy of the presented method is demonstrated. Different from the method proposed by Zhang et al. [16], the bioinspired method combines the Physarum algorithm and the penalty mechanism by modifying the internal mechanism

Examples
In this section, an example is introduced in order to prove the validity of the proposed method. The results are compared with the algorithm presented by Handler and Zang [14]. It is well known that the delay constrained least cost problem (DCLC) is a typical problem of CSP, which is widely used in network routing [57,58]. The DCLC problem is to find the least cost path in a graph while keeping the path delay below a specified value. Let = ( , ) be a directed network, where = 1, . . . , is the set of nodes and = ( , ) : , ∈ , ̸ = , is the set of edges. The edge ( , ) represents the direction from node to node . Each edge has two nonnegative attributes: and represent cost and transmission delay from node to node , respectively.
Here, ( , ) meets the following constraint: where ( , ) represents the set of all paths from node to node . As shown in Figure 3, it is a real network topological structure with source node 1 and sink node 23. It was frequently an example in fuzzy shortest path problem [54,59,60]. The delay and cost of each link in the network are randomly generated. The costs of links are randomly generated by normal distribution with mean value equaling 10 and standard deviation equaling 3. The delays of links are randomly generated using normal distribution with mean value equaling 10 and standard deviation equaling 5. The cost and delay of each link randomly generated by normal distribution are shown in Table 2.
The Scientific World Journal 7   In order to get significative constraint, the delay constraint is got as follows [20]: where is the least delay and is the delay of the path which has the least cost. They can be obtained by amoeba algorithm. The parameter is a problem specific ratio that measures the "tightness" of constraint (3) such that 0 < < 1. The more constraining (or "tight") the value of in constraint (3) is, the lower the value of will be. For simplicity, here, we let = 0.1.
Here, the least delay path is 1-3-8-13-19-22-23 and its delay is 44.0553. The least cost path is 1-4-11-17-20-23 and its delay is 54.1798 and its cost is 48.7660. The constraint 45.0680 can be obtained by using (15) The solution is exactly the same as the method [14] as well. Table 3 lists these shortest paths [14] and important paths.
As can be seen in Table 3, our searching order is different from Hanlder and Zangs. The reason is that the order of important path depends on the specific values of the parameters and . Generally speaking, smaller and larger lead to a bigger step when it searches for the following path in our method. Otherwise, the step will be small. It will search the alternative paths one by one. The parameters setting is in association with the scale of network and the parameter . In the case of a large scale of network, different parameters result in a big difference in executing time. The problem of how to set specific parameters will be further studied in the future work.

Computational Comparisons of the Solution Algorithms
In order to demonstrate the accuracy of proposed method, a number of experiments have been conducted on different networks. In addition, the CPU run-time has been compared between our method and the method proposed by Handler and Zang [14]. The main idea of the algorithm proposed by Handler and Zang is to use the shortest paths algorithm to find the next shortest path terminating with the first path satisfying the constraint. This algorithm requires solving the shortest path problem. The method finding the shortest loopless paths in a network proposed by Yen [61] and the Dijkstra algorithm for finding the shortest path are used to solve the shortest path problem.
Sixty-five networks were randomly generated. There are 5 different networks sizes based on the number of nodes in the networks. For each network in the same size, four networks are randomly generated. For each of the four 8 The Scientific World Journal random networks, three different attributes, but the same topologies networks, are generated. These attributes are randomly generated by normal distribution.
As for the random networks, we use the BA model [62] to generate a randomly directed network. In BA model, there are two assumptions: firstly, the network is formed by the continuous addition of new vertices to the system. Therefore, the number of vertices increases throughout the lifetime of the network. Secondly, when there is a new node added in the network, it tends to be connected with the vertex who has more connections.
Under this assumption, we use the following steps to construct the network.
(i) Beginning with a small network 0 , it has 0 nodes and 0 edges. At each time step, we add only one new node.
(iii) indicates the degree of node (1 ⩽ ⩽ ). The linking probability between nodes and +1 is defined as follows: Because the network is generated by this model, its degree distribution follows the power-law distribution with model = 2.9 ± 0.1. In this paper, we let 0 , 0 , and be equal to 3.
All parameter values are listed in Table 4.
The tests are executed on Intel(R) Core(TM) i3-2120 processor (3.30 GHz) with 2 GB of RAM under Windows Seven. Each instance is run for 3 times, and we compute the average executing time and average accuracy. Results of 65 test problems are summarized in Table 5.
After experiment, the accuracy is 95.38% in the proposed method. In all the 65 networks, there are only 3 networks we cannot get optimal solution when = 2, = 30. However, we have verified that once these two parameters are adjusted, the optimal solution can be obtained. As can be seen in Table 5, the method proposed by Handler and Zang [14] is faster than the approach proposed by this paper. However, the method combining amoeba model and penalty mechanism has three potential advantages. First of all, the CSP problem is an NP-hard problem and it is difficult for us to solve this problem. The method proposed by this paper is to contribute an innovative idea and to enrich the solution strategies for solving this NP-hard problem. Secondly, the reason why the proposed method is relatively slow is that the amoeba model in this paper is implemented in a tradition way. In fact, there are already many investigations on parallel amoeba algorithm which have dramatically decreased the computational time. For example, worthwhile works [63,64] have been achieved. Amoeba algorithm is constantly improving and it has a great potential for the development of the bioinspired artificial intelligence. It is certain that the computational time of the proposed method will be greatly reduced if there was a sophisticatedly parallel amoeba algorithm. At the current stage, we just aim at its methodology. Finally, the proposed method is easy to combine with other algorithms. By combining with other approaches, we can take the advantages of both of them. To sum up, although the executing time of the proposed method is relatively slow, it can be improved with the implementation of the parallel amoeba algorithm. What is more, the proposed method has many advantages over the existing algorithms.
The Scientific World Journal 9

Conclusion
In this paper, we propose an innovative method to solve the CSP, which employs the internal mechanism of the amoeba model. The solution to the CSP consists of two parts. Firstly, we extend Physarum algorithm into the shortest path problem in directed networks. Secondly, we combine the Physarum algorithm with a bioinspired rule to deal with the CSP. In addition, we have shown the validity of the proposed method by comparing with other existing approaches using an example in the DCLC problem.
In the future, on the one hand, we will further investigate how to determine the associated parameters, such as and . Currently, these parameters are given according to the many experimental results. In the near future, we will try to provide theoretical analysis to these parameters and to provide an interval to these parameters. On the other hand, we will explore how to implement the proposed method in parallel environment, which will greatly reduce the executing time of the proposed method.