Optimal algorithm for the demand routing problem in multicommodity flow distribution networks with diversification constraints and concave costs

1 OPTIMAL ALGORITHM FOR THE DEMAND ROUTING PROBLEM IN MULTICOMMODITY FLOW DISTRIBUTION NETWORKS WITH DIVERSIFICATION CONSTRAINTS AND CONCAVE COSTS Abstract.- Distribution problems are of high relevance within the supply chain system. In real life situations various different commodities may flow in the distribution process. Furthermore, the connection between production and demand centres makes use of complex mesh networks that can include diversification constraints to avoid overcharged paths. In addition, the consideration in certain situations of economies of scale gives rise to non-linear cost functions that make it difficult to deal with an optimal routing scheme. This problem is well represented by the multicommodity flow distribution networks with diversification constraints and concave costs (MFDCC) problem. Here we present an optimal algorithm based on the Kuhn-Tucker optimality conditions of the problem and capable of supplying optimal distribution routes in such complex networks. The algorithm follows an iterative procedure. Each iteration constructive solutions are checked with respect to the Kuhn-Tucker optimality conditions. Solutions consider a set of paths transporting all the demand allowed by its diversification constraint (saturated paths), a set of empty paths, and an indicator path transporting the remaining demand to satisfy the demand equation. The algorithm reduces the total cost in the network in a monotonic sequence to the optimum. The algorithm was tested in a trial library and the optimum was reached for all the instances. The algorithm showed a major dependency with respect to the number of nodes and arcs of the graph, as well as the density of arcs in the graph. Keywords: demand routing; concave cost; distribution network; Kuhn Tucker conditions 1. Introduction Distribution problems are of high relevance within the supply chain system. One of the principal root problems in the family is the transportation problem that was introduced in 1941 by Hitchcock (Hitchcock, 1941). This has been researched extensively due to its importance and the wide range of applications it has. The transportation problem involves a network with a group of production centres willing to send their commodity to a group of demand centres through a set of arcs with linear costs. However, in real life situations various different commodities may flow instead of just one commodity (e.g. Kengpol et al., 2011). Furthermore, each production centre is not connected with each demand centre using a single arc but rather with a mesh network of different arcs (Ubeda et al., 2010). Certain cases cannot be correctly modelled using a linear cost function and require a more complex way of modelling (Tsao et al. 2010). Sometimes more accurate results can be achieved with concave cost functions. This situation arises when dealing with economies of scale which are very common in many fields such as freight distribution and demand routing systems, and even other technological fields such as telecommunication networks. 2 In this paper, we consider a mesh network with a set of nodes as the origin of demand and a set of nodes as the demand destination. Each pair of nodes constitutes an origindestination pair which is known as a specific commodity, converting the problem into a multicommodity flow problem. Additionally, the costs in arcs are modelled by concave costs. Lastly, we propose the diversification constraints needed to prevent the arcs from becoming extremely loaded which can produce malfunction effects on the network and may arise when concave costs are involved as a result of the benefits offered by economies of scale. This type of constraint increases the survivability and reliability of the networks. Distribution problems with concave costs are NP-Hard and their optimal solution cannot be found in polynomial time. Several authors have proposed approximate methods to tackle this problem, such as the linear approximation by Thach (1992), the lagrangian relaxation by Larsson et al. (1994), or the dynamic programming approach by Zangwill (1968), Florian and Klein (1971) and Burkard et al. (2001). Metaheuristic approaches have also been tested although mainly in telecommunication network problems (Kapsalis et al, 1993; Dengiz et al, 1997; Altiparmak et al, 2003; Zhou and Gen, 2003). In a more general case, Yan and Luo (1999) developed a heuristic based on simulated annealing and threshold acceptance, and Altiparmark and Karaoglan (2008) proposed a tabu approach to solve the transportation problem with concave costs. Here, we have developed an optimal approach based on the optimality Kuhn-Tucker conditions of the mathematical formulation of the problem. The paper then defines the problem and its mathematical formulation in Section 2. Section 3 presents the Kuhn-Tucker optimality conditions for the demand routing problem in multicommodity flow networks with diversification constraints and concave costs that has been described in the previous section. Next, the solution based on the optimality constraints demand routing algorithm is proposed is Section 4. The results for a trial library with 27 generated instances are presented in Section 5 and the detailed procedure for an example network is developed. Finally, the main conclusions are presented in Section 6. 2. Problem definition and formulation Given a graph G=(N,E) where N is the set of nodes and E the set of arcs, a set of commodities K is considered and represents every origin-destination pair of demand. We define each arc of the network as e∈E and the set of paths connecting every origindestination pair is given by P(k). A subset of this set is given by the arc-disjoint paths, P d (k). Let γ k be the demand volume for each origin-destination pair, k, and δ k the diversification parameter for each pair k, so that the demand of pair k is routed along 1 k δ     different paths. As paths and arcs are pre-computed, ϕ eh is a binary parameter that determines if arc e is on path h. The variables of the problem are the demand fraction related to the origin-destination pair k that is routed on path h, which is a continuous variable given by p kh , and the total 3 flow transported on arc e, including the total amount of flow from all the commodities on the arc, which is a continuous variable given by x e . Owing to the nature of the problem described, numerous concave cost functions can be considered. We used the objective function given by ( ) e e e e c x d x = for arc e, where x e has been previously explained and d e means the per unit variable cost corresponding to the arc. We selected this function because it provides a widely acknowledged basis of analysis of scale economies’ effect in transportation problems, and because it is taken as a representative cost function very similar to many transportation cost functions in practice (see detailed arguments in LeBlanc 1976; Larsson et al. 1994; Yan and Luo 1999; and Altiparmak and Karaoglan 2008). The cost concave function is, therefore, a differentiable function. Note that the objective function is used to provide comparable numerical examples in the subsequent results section. The selection of this function does not condition the methodology and developed algorithm, because the objective function only affects the procedure for the calculation of the derivative providing numerical data. As a result, the demand routing problem in multicommodity flow networks with diversification constraints and concave costs (MFDCC problem) can be formulated as: ( ) ( ) ( ) Minimize , 1 , 2 , e e e E e eh kh k K h ( k ) kh k h ( k ) kh k k c x s.t.: x p e E p k K p h ∈ ∈ ∈Ρ ∈Ρ = ϕ ∀ ∈ = γ ∀ ∈ ≤ δ γ ∀ ∈Ρ ∑ ∑ ∑ ∑ ( ) ( ) , 3 0 , 0 d e kh k k K x p ∀ ∈ ≥ ≥ Constraint (1) calculates the total amount of flow on the arc e; constraint (2) ensures that the demand for every origin-destination pair is met; and constraint (3) ensures that the demand is routed along 1 k δ     alternative paths. Set P d (k) represents all the feasible disjoint paths for each origin-destination pair. The “disjoint” concept is important and can be weakly imposed on the arcs or strongly imposed on the nodes. Here, we have chosen the arc diversification concept because in most cases it is normally only necessary to take this weak constraint into consideration. Using an example from telecommunications, disjoint arcs can be used to model reality, since the optical crossconnect (OXC) in an all-optical WDM mesh network is seldom broke. The OXC is a device which switches an optical signal from an incoming fibre to an outgoing fibre on the same wavelength. So, constraint (3) is equivalent to:  ( ) ( ) , d kh k p i, j h , h P k : k K ≤ δ ∀ ∈ ∈ ∈ , constraint for arc-disjoint paths  ( ) , d kh k p j h , h P k : k K ≤ δ ∀ ∈ ∈ ∈ , constraint for node-disjoint paths depending on what is being considered and as explained above, we are looking at 4 constraints for arc-disjoint paths. 3. Kuhn Tucker conditions for the demand routing problem in multicommodity flow networks with diversification constraints and concave costs Every solution verifying the necessary and sufficient Kuhn-Tucker conditions of any problem will be the optimal solution of problem. This fact is proven for linear or nonlinear problems since they are global optimality conditions. Additionally, duality theory can be regarded as a particular case of the Kuhn-Tucker conditions for linear problems. See classic theory such as Hillier and Lieberman (1977), Nocedal and Wright (1999), Avriel (2003) for detailed explanations on Kuhn-Tucker conditions. Appendix 1 depicts a general overview on the Kuhn-Tucker necessary optimality conditions and their detailed corresponding application to MFDCC problem, and appendix 2 depicts the equivalent general overview and detailed application of the Kuhn-Tucker sufficient optimality conditions to the MFDCC problem. Next, we re-formulate MFDCC problem to specify the objective function in terms of p kh variables. This can be done because p kh and x e are inter-related by constraint (2). Therefore, constraint (2) is included in the objective function to re-formulate the problem and subsequently formulate the Kuhn-Tucker necessary optimality conditions: ( ) ( ) ( ) ( ) = 0 , (free multipliers) 0 , , (non-negative multip ∈ ∈ ≡ ≡ γ − ∀ ∈ ← µ − δ γ ≤ ∀ ∈ ∀ ∈ ← ∑ ∑ kh h P(k) e e e E k k d kh k k kh p Minimize c x Minimize C x Minimize C p s.t.: k K p · h P k k K v ( ) liers) 0 , , (non-ne gative multipliers) − ≤ ∀ ∈ ∀ ∈ ← d kh kh p h P k k K u Kuhn-Tucker conditions state that all ( ) d h P k ∈ and k K ∈ can be written as: kh kh k kh k kh k kh kh kh C(p) = u v p ( · ) · = 0 pv · = 0 pu - - ∂µ ∂ γ δ - (4) By calculating the cost function variation related to the demand fraction routed along path h, we replace the expression of the total cost in the network C(p) with the sum of the costs of the flows in the arcs, and then we substitute the value of the flow x e accordingly with equation (1): 5 ϕϕ eh ee Ee ehkh P(k)hKk kh ee Ee kh e e ee Ee kh ee Ee kh ee Ee kh ) x ( c = p p ) x ( c = = p x x ) x ( c = p ) x ( c = p ) x ( c = p C(p) ' · '∑           ∑∑ ∂ ∂ ∑ ∂ ∂ ∂ ∂ ∑ ∂ ∂ ∑ ∂ ∑ ∂ ∂ ∂ ∈∈∈∈ ∈∈ ∈ (5) that corresponds to the sum of cost variations on the arcs belonging to path h. The optimality conditions are thus enforced on the increase of the cost and not on its absolute value. By analysing these conditions, we can differentiate the meaning of the conditions and three different types of paths are found:  Case 1: the origin-destination pair does not use path h to establish a communication. kh kh kh kh k kh e kh eeh k e E e eeh k e E = 0 0 and = 0 pu v C(p) - = 0 u p c' ( ) - = 0 x u c' ( ) x ∈ ∈ ⇒≥ ∂≥ µ ∂ ≥ ϕ µ ≥ ϕ µ ∑ ∑ (6) Therefore, the sum of the cost variations corresponding to the arcs on the empty or void paths that are not used by the origin-destination pair k is greater than k µ .  Case 2: the origin-destination pair uses path h to establish a communication. Here, two possibilities arise: o Case 2.a: the arc transports all the flow allowed by the diversification constraint. There are 1 1 k     −       δ     paths in this case. In this case, the Kuhn- Tucker conditions are reduced to: k kh kh kh kh k kh e kh eeh k e E e eeh k e E = = 0 and 0 pu v C(p) - = - 0 v p c' ( ) - = - 0 x v c' ( ) x ∈ ∈ ⇒≥ δ ∂≤ µ ∂ ≤ ϕ µ ≤ ϕ µ ∑ ∑ (7) The sum of the cost variations corresponding to the arcs on the saturated paths used by the origin-destination pair k, i.e. ( ) d h P k ∈ , is upper bounded by k µ . 6 o Case 2.b: the arc transports less flow than the limitation imposed by the diversification constraint. There is only one path that verifies this situation, which is: k kh kh kh k kh eeh k e E < = 0 and = 0 pu v C(p) - = 0 p c'( ) = x ∈ ⇒ δ ∂µ ∂ ϕ µ ∑ (8) The sum of the cost variations corresponding to the arcs on the path that transports a positive amount of flow but below the maximum value imposed by the diversification constraint related to the origin-destination pair k is equal to k µ . This path will be known as the indicator path. This is a concept that appears in similar terms in the traditional transportation problem (when links in the graph are grouped into three sets; a set of slacks equal to zero; a set of slacks equal to the capacity, that is empty arcs; and one link with an amount of flow equal to the difference between the demand and the sum of all the flow in arcs with slack equal to zero). We originally introduced this term in Cortes et al. (2006). The existence of only one indicator path is also intuitive. Note that if there was a cheaper path, that path would be full, and if there was a more expensive path it should be empty, and so the path in the middle separating cheap paths from expensive paths should transport the rest of demand to satisfy the demand equation. There could be only two particular cases: (i) the rest of demand is equal to the capacity, so it would be full, (ii) there is a draw between two paths, corresponding to a case of alternative solutions, but the solution cost would be the same, and these two cases do not condition the algorithm procedure that will be later detailed in section 4, because either of the drawing paths could act as indicator as a matter of fact. In summary:  For every origin-destination pair, k, cheap, expensive and indicator paths can be considered. The cheap paths are saturated and transport as much flow as possible, i.e.: kh k k p · = δ γ . The expensive paths are empty and do not transport flow, i.e.: 0 kh p = The indicator path transports the necessary amount of flow to satisfy the demand equation (eq. 2), which is less than k k · δ γ and more than 0.  The transported product that flows through the communication paths in order to satisfy the demand of all origin-destination pairs will be optimum if the “separator cost concept” is given by e k eh e E c'( ) x ∈ = µ ϕ ∑ in the indicator path. In this way, the Kuhn-Tucker multipliers include the following meanings: 7 ' h k i c ∼ µ It is a concept associated with the value of the derivative of the cost function along the indicator path. ( ) hshkh c -c v i ''∼ It is a concept associated with the deviation between the values of the derivative of the cost function of the indicator path and each saturated path. The Kuhn-Tucker optimality conditions determine that this deviation must be positive for every ( ) d h P k ∈ and k K ∈ . ( ) iv hhkh c -c u ''∼ It is a concept associated with the deviation between the values of the derivative of the cost function of each empty path and the indicator path. The Kuhn-Tucker optimality conditions determine that this deviation must be positive for every ( ) d h P k ∈ and k K ∈ . 4. Demand routing algorithm After having explained the Kuhn-Tucker optimality conditions above, this section describes an optimal iterative algorithm that can be used to solve the demand routing problem in multicommodity flow networks with diversification constraints and concave costs. The algorithm converges to the optimal process by reducing the non-feasibility after every iteration. The algorithm needs to pre-compute all the paths of the graph between a pair of nodes. This is done with an algorithm based on a depth first search routine that is shown in appendix 4. Although the problem of finding all the paths between every pair of nodes in a graph is a NP-Hard problem, we consider sparse networks, so that results can be obtained quickly as the result section shows next. Step 1. Initialisation. Greedy heuristic The procedure is based on the k shortest disjoint path algorithm provided by Kleinberg (1996). 1. Order the origin-destination pairs according to the demand volume, γ k . Subsequently, following that order: 2. For each origin-destination pair, k: a. For a number of times equal to 1 k δ     -1 i. Calculate the shortest path between the origin o(k) and the destination d(k).This is a saturated path h s ∈ P d (k). ii. Assign a volume of demand equal to δ k ·γ k . iii. Erase the path to search for a new disjoint path. b. Calculate the shortest path between the origin o(k) and the destination d(k), and assign the remaining demand to satisfy equation (2). This is the candidate indicator path h i ∈ P d (k). In the greedy heuristic the costs of the arcs are simplified to ( ) e e e c x d = in order to calculate the shortest paths. Dijkstra’s algorithm can, therefore, be applied. This initial flow load provides a determined level of flow along the paths, p kh0 , corresponding to a 8 specific level of flow load on the arcs of the network, x e0 . We used this greedy heuristic with the objective of providing an initial set of disjoint paths. It starts selecting a set of edge-disjoint paths that provides an initial set that is only considered as the seed of the algorithm. After that, the iterative algorithm (next steps of the iterative routine that forces to accomplish the Kuhn-Tucker conditions) guides the flow in the paths from that seed to optimality by testing the optimality conditions for all the paths of the graph connecting such pair k. The main reason to use this shortest path based initialization algorithm is to generate a seed quickly. Furthermore, as the shortest path tries to minimise the cost between a pair of nodes, this can be taken as a reasonably good initial solution following the arguments detailed in Kleinberg (1996). Step 2.Checking the Kuhn-Tucker optimality conditions The Kuhn-Tucker conditions are checked for every origin-destination pair, k. The condition is given by (9): ( ) ( ) ( ) ( ) ( ) ( ) { } ( ) ( ) ( ) 0 0 0 0 for every saturated path indicator path empty or void path where s i v h kh h kh h kh d s i v h kh e e eh e E c' p c' p c' p h h ,h ,h P k c' p c' x ∈ ≤ ≤ = ∈ = ⋅ϕ ∑ (9) After the saturated and indicator paths have been obtained (as a result of step 1 or the iterative process), we have to check the conditions for all the empty paths of the graph (it has to be taken into account that the condition must be verified for all the empty paths in the graph as well. These do not necessarily have to be disjoint paths). The iterative process comes to an end if the optimality conditions are checked; otherwise, the procedure follows step 3. Step 3. Iteration All the origin-destination pairs that do not check the optimality conditions shape the set K (we say that k K ∈ if pair k does not verify the Kuhn-Tucker conditions, i.e. eq. 9). For these pairs, the set H(k) is defined as the paths that do not check the condition of such pair, k K ∈ . Next procedure is followed: 1. For every k K ∈ , paths in H(k) are arranged in increasing order according to the derivative of the cost function assessed for the previous (step 2) demand routing pattern, 0 kh p . That is, ( ) ( ) 0 0 ' ' h kh e e eh e E c p c x ∈ = ⋅ϕ ∑ is calculated for every h∈ H(k), for every k K ∈ . 2. Origin-destination pair k * is selected as it is the pair that provides the highest absolute deviation with respect to the optimality conditions, as calculated in (10). 9 ( ) ( ) ( ) { } ( ) ( ) { } { } * ' 0 ' 0 ' 0 ' 0 : max max ,max s i i v s v h kh h kh h kh h kh h H k h h k k K c p c p c p c p ∈     = ∈ − −       (10) 3. For * k K ∈ , the maximum demand allowed by constraint (3), δ k ·γ k , is assigned to the “ ( ) 1 1 k δ −     ” disjoint paths with lower value of ( ) * ' 0 h k h c p . These are saturated paths for the new iteration. 4. The next disjoint path, following the order determined by the derivative, transports the rest of the required demand to satisfy constraint (2). This is the indicator path for the new iteration. Subsequently, the procedure follows step 2 to check the optimality of the new demand routing pattern. Although we have tackled for this algorithm only the Kuhn-Tucker optimality necessary conditions, the sufficient condition is always checked for every intermediate feasible (non-optimal) solution of the algorithm given the special characteristics of the problem. In fact, the sufficiency condition is always checked independently from its optimality condition due to every solution composed of an indicator path, a number of ( ) 1 1 k δ −     saturated paths and a set of empty paths checks it (see appendix 2 for a detailed argument and proof). Therefore, a feasible solution [ ] * kh p which also checks the necessary optimality conditions previously described in equation (9) will be the global optimum of the MFDCC problem. See Hanson and Mond (1987) for detailed arguments and explanations for necessary and sufficient conditions in optimisation problems. In addition, the procedure guarantees that the algorithm follows a monotonic sequence that reduces the cost of the objective function step-by-step. Proof and arguments are shown in appendix 3 arguing how the total cost is reduced after each iteration. 5. Computational results 5.1. Application example We present here the next case as an application example to demonstrate the performance of the algorithm. Consider the eight nodes, and eleven arcs network in Figure 1. Data over the arcs represent the value of d e in the cost function ( ) e e e e c x d x = . The remaining data appear in Table 1. 16 linear relation that shoots up the distribution cost in the network, being the parameter with a more significant repercussion (see Figure 2). Figure 2. Relation between cost and main design parameters for instances 1 to 26 The time required by the algorithm to get the optimal solution was mainly dependent on the number of arcs and nodes (Figure 3). The time consumption increased at a greater rate than the size of the problem as expected for a NP-Hard problem such as the analysed. Instances 8, 16 and 18 showed the higher times to reach the optimum of the problem because they represent the major dimension problems. So, worse results were obtained for graphs which combined a large number of nodes, arcs, origin-destination pairs, and demand. 0 500 1000 1500 2000 2500 0 5 10 15 20 25 30 35 Cost Nodes 0 500 1000 1500 2000 2500 0 20 40 60 80 Cost Arcs 0 500 1000 1500 2000 2500 0 50 100 150 200 250 300 350 400 Cost Number of origin-destination pairs 0 500 1000 1500 2000 2500 0 2000 4000 6000 8000 Cost Total demand volume 0 500 1000 1500 2000 2500 0 50 100 150 200 250 300 Cost Diversification 17 Figure 3. Relation between computational time and main design parameters for instances 1 to 26 The number of iterations of the algorithm is strongly linked to the computational time of the algorithm, and it is mainly affected by the density of the graph, as well as the number of arcs and nodes too. So, the problems that take a greater time in the optimum convergence were the larger problems and also those problems with a higher number of pre-computed paths, which required a greater number of iterations (Table 8). Following this later aspect, Figure 4 shows the relation between the number of paths that were pre-computed for each graph, and the computational time and the number of iterations of the algorithm. The figure shows that the increase in terms of computational time and number of iterations follows an increasing tendency with a rate lower than when it is analysed with respect to other design parameters (see how the increasing rate is higher when the computational time is analysed with respect to other design parameters in Figure 3). This fact contributes to validate and enhance the real potential applicability of the approach. 0 500 1000 1500 2000 2500 3000 0 5 10 15 20 25 30 35 Time Nodes 0 500 1000 1500 2000 2500 3000 0 10 20 30 40 50 60 70 80 Time Arcs 0 500 1000 1500 2000 2500 3000 0 100 200 300 400 Time Number of origin-destination pairs 0 500 1000 1500 2000 2500 3000 0 1000 2000 3000 4000 5000 6000 7000 8000 Time Total demand volume 0 500 1000 1500 2000 2500 3000 0 50 100 150 200 250 300 Time Diversification 18 Figure 4. Time and number of iterations evolution with respect to the number of paths in the graph for instances 1 to 26 Finally, Figure 5 provides an interesting analysis related to the algorithm convergence. The optimising curve depicts a monotonic decrease of the total cost value that the algorithm achieved, as shown in the figure. We have selected instances 1, 4, 8 and 18 as representative examples to appreciate the algorithm convergence. It is interesting to note how this phenomenon that is argued in appendix 3 is empirically checked when solving the instances of the trial library. Figure 5. Convergence to the optimum curves for graph instances 1, 4, 8 and 18 5.4. Analysis of the results for dense graphs Finally, we analyse the main parameter dependencies related to dense graphs. Starting with the objective function, the cost analysis in dense graph follows a very similar structure to the analysis in sparse or medium density graphs, depicting a clearly linear tendency. However, the linear dependency is stronger with respect to the nodes than the arcs (see left side of Figure 6). Furthermore, the dependency with respect to the diversification parameter ( k δ ∑ ) than with respect to the number of origin-destination 0 500 1000 1500 2000 2500 3000 0 5000 10000 15000 20000 Time Number of paths in the graph 0 10 20 30 40 50 60 0 5000 10000 15000 20000 Number of iterations Number of paths in the graph 2325,00 2330,00 2335,00 2340,00 2345,00 2350,00 2355,00 0 1 2 3 4 5 6 7 8 9 101112131415161718192021222324252627282930313233343536373839404142434445464748495051 CT(total cost) iterations Graph #1 CT(Nº de Iteraciones) 1985,74 1985,57 1985,28 1983,23 1982,99 1981,03 1980,98 1980,21 1980,19 1980,14 1979,48 1979,47 1979,33 1979,28 1979,18 1978,96 1978,89 1978,59 1977,95 1977,00 1976,99 1976,93 1976,90 1976,86 1976,00 1977,00 1978,00 1979,00 1980,00 1981,00 1982,00 1983,00 1984,00 1985,00 1986,00 1987,00 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 CT(Total cost) Iterations Graph #4 CT(Nº d e Iteraciones) 1844 1846 1848 1850 1852 1854 1856 1858 1860 1862 1864 1866 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 CT(Total cost) Iterations Graph #8 CT(Nº de Iteraciones) 483,07 482,12 479,65 479,05 478,74 478,50 479,00 479,50 480,00 480,50 481,00 481,50 482,00 482,50 483,00 483,50 01234 CT(Total cost) Iterations Graph #18 CT(Nº de Iteraciones) 19 pairs (see right side of Figure 6). Figure 6. Relation between cost and main design parameters for instances 27 to 30 The analysis of the computational time against the main parameters is also interesting and relevant. Here again major dependencies appear for nodes and diversification. Also, a higher than linear tendency is appreciated with respect to number of nodes (left side of Figure 7), and especially with respect to the diversification parameter (right side of Figure 7). Figure 7. Relation between time and main design parameters for instances 27 to 30 Finally, similar comments to those provided for sparse and medium density graphs attending to the computational time and number of iterations with respect to the number of paths is also valid for dense graphs. Here again, the increase both in computational time and number of iterations is not so high with respect to the number of paths in the graphs as with respect to other design parameters as can be viewed in Figure 8. Figure 8. Time and number of iterations evolution with respect to the number of paths in the graph for instances 27to 30 0 500 1000 1500 2000 2500 3000 3500 4000 0 10000 20000 30000 40000 50000 60000 70000 Time Number of paths in the graph 0 10 20 30 40 50 60 0 10000 20000 30000 40000 50000 60000 70000 Number of iterations Number of paths in the graph 20 6. Conclusions In this paper we have proposed a novel optimal algorithm based on the Kuhn-Tucker optimality conditions to determine the optimum of the demand routing problem in multicommodity flow networks with diversification constraints and concave costs. The aforementioned problem is considered NP-Hard because of its concave cost function. In addition, it is an ambitious problem that takes most of the elements from traditional concave cost transportation problems literature such as the objective function and the multicommodity constraints, but also includes new elements as the diversification constraints that have not been considered in concave cost paper literature, depicting a more complex problem that those tackled in most of the concave cost transportation papers. Such diversification constraints are required many times to model situations that require diversifying elements in alternative routes of a distribution network, such as freight in transportation networks or cells in telecommunication networks, modelling a fact that appear in a number of real practice situations. We have presented an iterative algorithm based on the Kuhn-Tucker optimality conditions of the problem. Feasible solutions composed of an indicator path, a number of ( ) 1 1 k δ −     saturated paths and a set of empty paths are analysed for checking the Kuhn-Tucker necessary conditions each iteration, since the sufficient condition is always checked for every intermediate feasible solution defined as the previous one and independently from its optimality. The algorithm enables the optimum to be reached after a number of iterations. The method is a constructive step-by-step method because after each iteration the total cost in the network is reduced in a monotonic sequence. Therefore, given a MFDCC problem, a sufficiently high number of iterations of the algorithm will allow reaching the optimum of the problem. It is clear that as the MFDCC problem is NP-Hard, the number of iterations of the algorithm increases in a non-polinomial way when the size of the problem increases. For very great size problems this exact algorithm could not be efficient. The algorithm showed a major dependency with respect to the number of nodes and arcs of the problem as well as the level of diversification, which define the size of each instance. Also the density of the graph depicted a clear influence for problems of the same size increasing significantly the computational time. The number of origindestination pairs, as well as the number of alternative paths for connecting each origindestination pair revealed a not so relevant dependency levelling off its tendency. Also, it is relevant to note that tendencies were appreciated in a similar way both for sparse, medium density and dense graphs. The optimum was obtained for all the instances studied, and the time consumption was feasible when analysed from a non-real-time perspective, which is normal in network design problems. Although in real time situations, the algorithm may not be suitable depending on the time response required, alternatives such as stop criteria after a number of iterations or time consumption could be considered because the method guarantees a monotonic approach to the optimum. 21 Acknowledgements The authors acknowledge an anonymous referee for his/her helpful comments on earlier versions of this manuscript that have significantly contributed to enrich the quality of the paper in its final form. Appendix 1 Appendix 1 states the Kuhn-Tucker optimality necessary conditions for every optimisation problem, and its application to the MFDCC problem. Consider the next generic optimization problem: The Kuhn-Tucker conditions are given by: The Kuhn-Tucker conditions for the MFDCC problem are given by: A solution verifying the necessary condition is only optimum if it verifies also the sufficiency condition in problems with non-convex objective functions. Appendix 2 Appendix 2 deals with the Kuhn-Tucker optimality sufficient conditions. In addition, here is proven that every solution for the communication of K commodities, each of them with a demand γ k and forced to a diversification given by δ k , that is composed by a set of ( ) 1 1 k δ −     saturated paths, an indicator path and being the rest void paths verifies the Kuhn-Tucker sufficient conditions. Thus, every solution associated to each iteration of the algorithm described in section 4 will verify the Kuhn-Tucker sufficient condition.      (13) s.t.:      ≤ 0 , ∀  (14) ℎ     = 0 , ∀  (15)  ∈ ℛ  ∇     + # $  ∇      +  # %  ∇ ℎ      = 0 (16)      ≤ 0 , ∀  ℎ     = 0 , ∀  (17) $  ·      = 0 , ∀  (18) $  ≥ 0 ∀  and %  *++ (19) ,-  .  , .  ℎ − %  = $  ℎ − /  ℎ (20)  .  ℎ − 0  · 1   · /  ℎ = 0 (21) .  ℎ · $  ℎ = 0 (22) $  ℎ , /  ℎ ≥ 0 and %  *++ (23) 22 2 ℎ 3 ′ > 2 ℎ  ′ Let consider the generic problem given by equations (13-15). The sufficiency condition for such generic problem is given by (24): Where: The sufficiency condition for the MFDCC problem is given by (25): And this expression must be satisfied in set V(x) that is defined as: Therefore, any flow pattern verifying equations (20 to 23) verifies the sufficient condition too, and therefore the necessary conditions becomes the global optimality conditions for the specific structure of the analysed problem. Appendix 3 Appendix 3 proves that the algorithm described in section 4, guarantees the reduction of the cost of the objective function after each step following a monotonic sequence to optimality. Suppose that after iterating, the Kuhn-Tucker conditions cannot be satisfied for every origin-destination pairs. Therefore, the origin-destination pair k * is given by (9). For such pair three options can appear: a) the derivative is higher in a saturated path than in the indicator path 6 7 ∇  8   , $ , %  6 ≥ 0 ∀ 6 ∈ 9    (24) ∇  8   , $ , %  = ∇     + # $  ∇      +  # %  ∇ ℎ      and: 9=6∈ℛ  :∇ℎ   7 ·6=0 ,∀ ; ∇   7 · 6=0 ,∀∈ where A(x) is the set of active constraints (that is, verifies constraint      ≤ 0 with zero slack). 6 7 < , 2 -  .  , .  ℎ 2 = 6 ≥ 0 , ∀ 6 ∈ 9    (25) 9    = > 6 ∈ ℛ  : # 6  ℎ = 0 , ∀  ∈  ℎ∈? and 6  ℎ = 0 for those paths ℎ ∈ ?    so that .  ℎ = 0   saturated paths  or so that .  ℎ = 0  empty paths  K which is verified by every feasible solution associated to the corresponding iteration of the algorithm described in section 4. In effect, equation (25 ) is verified in every saturated and empty path because 6 ℎ =0 for such paths , and the indicator path verifies also 6 ℎ =0 because ∑6 ℎ =0,∀∈ ℎ∈? and there is only one indicator path in each iteration, so being the other paths (saturated and empty) equal to zero, the indicator path must be equal to zero too. 23 2 ℎ  ′ > 2 ℎ / ′  2 ℎ 3 ′ − 2 ℎ  ′  and  2 ℎ  ′ − 2 ℎ / ′   2 ℎ 3 ′ − 2 ℎ  ′  >  2 ℎ  ′ − 2 ℎ / ′  .  2 ℎ 3 ′ − 2 ℎ  ′  <  2 ℎ  ′ − 2 ℎ / ′  . b) the derivative is higher in the indicator path than in a void path c) Both situations occur but one of the expressions: is higher than the other, so we select the corresponding case. So we need to analyse the recirculation of flow for the two cases given below as option A and option B. Option A: It corresponds to recirculate the flow from the saturated path to the indicator path for cases a) or c) when Total cost is given for current situation (26) and, after recirculating, new situation (27): So, there is an unchanging part given by the flown on the links for every link not belonging to the saturated or indicator path, and a changing part given by the links in the saturated and indicator paths. Calculating the variation of cost between new and current situations, equation (28) is reached. Option B It corresponds to recirculate the flow from the indicator path to the void path for cases b) or c) when Total cost is given for current situation (29) and, after recirculating, new situation (30): - 2  ?  = # 2 +   +  + ∉  ℎ 3 ; ℎ   + # 2 +   +  + ∈ ℎ 3 + # 2 +   +  + ∈ ℎ  (26) -   ?  = # 2 +   +  + ∉  ℎ 3 ; ℎ   + # 2 +   + − △   + ∈ ℎ 3 + # 2 +   + + △   + ∈ ℎ  (27) -   ?  − - 2  ?  = O # 2 +   +  + ∈ ℎ 3 + # 2 +   + − △   + ∈ ℎ 3 P + O # 2 +   +  + ∈ ℎ  + # 2 +   + + △   + ∈ ℎ  P (28) From the checking of the Kuhn-Tucker conditions we know that 2 ℎ 3 ′ > 2 ℎ  ′ (that was the reason for non-optimality). Therefore, the decrease in the cost experimented when reducing flow in the saturated path corresponding to the evaluation of the first term of the sum is higher than the increase in the cost experimented when increasing flow in the indicator path corresponding to the evaluation of the second term of the sum (this argument is based on 2 ℎ 3 ′ > 2 ℎ  ′ ). And thus, -   ?  − - 2  ?  < 0 producing a decrease of the total cost after recirculating the flow consequently . - 2$**+R  ?  = # 2 +   +  + ∉  ℎ  ; ℎ /  + # 2 +   +  + ∈ ℎ  + # 2 +   +  + ∈ ℎ / (29) 24 So, there is an unchanging part given by the flown on the links for every link not belonging to the indicator or void path, and a changing part given by the links in the indicator and void paths. Calculating the variation of cost between new and current situations, equation (31) is reached. These arguments prove that after each iteration the total cost in the network is reduced. Therefore, given a MFDCC problem, a sufficiently high number of iterations of the algorithm allow reaching the optimum of the problem. It is clear that as this problem is NP-Hard, the number of iterations of the algorithm increases in a non-polinomial way when the size of the problem increases. Appendix 4 This appendix includes the procedure followed to find all the paths in a graph. There are several alternatives to find all the paths between a pair of nodes in a graph. Most of them follow Depth-First Search (DFS) or Breadth-First-Search (BFS) structures. We followed a procedure based in the DFS algorithm. This is a well-known algorithm that can be described as next: Input: A graph G and a vertex v of G Output: A labelling of the edges in the connected component of v as discovery edges and back edges procedure DFS(G,v): label v as explored for all edges e in G.adjacentEdges(v) do if edge e is unexplored then w ← G.adjacentVertex(v,e) if vertex w is unexplored then label e as a discovery edge recursively call DFS(G,w) - +S  ?  = # 2 +   +  + ∉  ℎ  ; ℎ /  + # 2 +   + − △   + ∈ ℎ  + # 2 +   + + △   + ∈ ℎ / (30) -   ?  − - 2  ?  = O # 2 +   +  + ∈ ℎ  + # 2 +   + − △   + ∈ ℎ  P + O # 2 +   +  + ∈ ℎ / + # 2 +   + + △   + ∈ ℎ / P (31) From the checking of the Kuhn-Tucker conditions we know that 2 ℎ  ′ > 2 ℎ / ′ (that was the reason for non-optimality). Therefore, the decrease in the cost experimented when reducing flow in the indicator path corresponding to the evaluation of the first term of the sum is higher than the increase in the cost experimented when increasing flow in the void path corresponding to the evaluation of the second term of the sum (this argument is based on 2 ℎ  ′ > 2 ℎ / ′ ). And thus, -   ?  − - 2  ?  < 0 producing a decrease of the total cost after recirculating the flow consequently . 25 else label e as a back edge References 1. Altiparmak F, Dengiz B, Smith AE., 2003. Optimal design of reliable computer networks: A comparison of metaheuristics. Journal of Heuristics 9, 471-487. 2. Altiparmak F, Karaoglan I., 2008. An adaptive tabu-simulated annealing for concave cost transportation problems. Journal of the Operational Research Society 59, 331-341. 3. Avriel, M., 2003. Nonlinear Programming: Analysis and Methods. Dover Publishing. 4. Burkard RE, Dollani H, Thach PT., 2001. Linear approximations in a dynamic programming approach for the uncapacitated single source minimum concave cost network flow problem in acyclic networks. Journal of Global Optimization 19, 121-139. 5. Dengiz B, Altiparmak F, Smith AE., 1997. Local search genetic algorithm for optimal design of reliable networks. IEEE Transactions on Evolutionary Computation 1; 179–188. 6. Cortes, P., Muñuzuri, J., Onieva, L., Larrañeta, J., Vozmediano, J.M., Alarcon, J.C., 2006. Andalucía assesses the investment needed to deploy a fiber-optic network. Interfaces 36 (2), 105-117 7. Florian M, Klein M., 1971. Deterministic production planning with concave cost and capacity constraints. Journal of Management Science 18, 12–20. 8. Hanson MA, Mond B., 1987. Necessary and sufficient conditions in constrained optimization. Mathematical programming 37; 51-58. 9. Hillier, F.S., Lieberman, G.J., 1974. Operations Research. Holden-Day Inc. (San Francisco, California, US) 10. Hitchcock, F.L., 1941. The distribution of a product from several sources to numerous locations. Journal of Mathematical Physics 20, 224–230. 11. Kapsalis A, Rayward-Smith VJ, Smith GD., 1993. Solving the graphical Steiner tree problem using genetic algorithms. Journal of the Operational Research Society 44, 397– 406. 12. Kengpol, A., Meethom, W., Tuominen, M., 2012. The development of a decision support system in multimodal transportation routing within Greater Mekong sub-region countries. International Journal of Production Economics 140(2), 691–701. 13. Kleinberg, J. 1996. Approximation Algorithms for Disjoint Paths Problems. Ph.D Thesis, Department of EECS, MIT.