## Serviços Personalizados

## Journal

## Artigo

## Indicadores

- Citado por SciELO
- Acessos

## Links relacionados

- Citado por Google
- Similares em SciELO
- Similares em Google

## Compartilhar

## Pesquisa Operacional

##
*versão impressa* ISSN 0101-7438

### Pesqui. Oper. vol.33 no.2 Rio de Janeiro maio/ago. 2013 Epub 16-Jul-2013

#### http://dx.doi.org/10.1590/S0101-74382013005000004

**Global optimization of capacity expansion and flow assignment in multicommodity networks **

**Ricardo Poley Martins Ferreira ^{I,}^{*}; Henrique Pacca Loureiro Luna^{II}; Philippe Mahey^{III}; Mauricio Cardoso de Souza^{IV} **

^{I}Departamento de Engenharia Mecânica, Universidade Federal de Minas Gerais - UFMG, Brazil. Phone: +55(31) 3409-3503. E-mail: poley@demec.ufmg.br

^{II}Instituto de Computação, Universidade Federal de Alagoas - UFAL, Maceió, Brazil. E-mail: henrique.luna@pq.cnpq.br

^{III}LIMOS UMR 6158 CNRS, Université Blaise Pascal, Clermont-Ferrand, France. E-mail: philippe.mahey@isima.fr

^{IV}Departamento de Engenharia de Produção, Universidade Federal de Minas Gerais - UFMG, Brazil. E-mail: mauricio@dep.ufmg.br

**ABSTRAC**

This paper describes an exact algorithm to solve a nonlinear mixed-integer programming model due to capacity expansion and flow assignment in multicommodity networks. The model combines continuous multicommodity flow variables associated with nonlinear congestion costs and discrete decision variables associated with the arc expansion costs. After establishing precise correspondences between a mixed-integer model and a continuous but nonconvex model, an implicit enumeration approach is derived based on the convexification of the continuous objective function. Numerical experiments on medium size instances considering one level of expansion are presented. The results reported on the performance of the proposed algorithm show that the approach is efficient, as commercial solvers were not able to tackle the instances considered.

**Keywords:** capacity expansion, flow assignment, global optimization, implicit enumeration, multicommodity flow problems.

**1 INTRODUCTION **

We consider a model for the joint problem of capacity expansion and flow assignment in multicommodity flow networks which takes into account congestion effects. The resulting optimization problem relies on the combination of two conflicting criteria: the expansion and the congestion cost functions. Capacity expansion cost functions are discrete as only a finite number of capacity sizes are available. On the other hand, congestion cost functions are nonlinear convex increasing functions as they try to capture queueing effects on the network [3, 4, 18, 23]. For instance, Ishfaq & Sox [18], in the context of an intermodal hub network, deal with shipment delays due to limited resources at logistics hubs for measuring service performance. The authors conducted a study of strategies to integrate hub operation queuing model and the hub locationallocation model on a 25-city road-rail intermodal logistics network. Our aim is to study the general multicommodity flow problem under such assumptions, with no particular applicationin mind.

Typically, the two decision levels of capacity expansion and flow assignment are decoupled to treat the corresponding difficulties separately, and few works, at least those considering queueing effects, have been done on the joint optimization problem. Gerla & Kleinrock [14] stated the general Capacity and Flow Assignment problem (CFA), decomposed it into simpler subproblems, and then suggested a heuristic procedure that alternates a capacity assignment phase with a flow assignment phase. Similar approaches have been proposed in successive papers by Gavish and co-authors [11, 12] and Gerla *et al. *[15]. Most of that literature concerns heuristic procedures, but the paper by Mahey *et al. *[20] is an exception for this rule, showing that a generalized Benders decomposition method can find exact solutions for the CFA problem.

The capacity expansion problem is a special case of CFA where initial capacities are already installed on each arc of the network. Given a traffic requirement matrix between origin-destination pairs, the problem consists in jointly deciding which arc capacities, if any, should be expanded and the flow assignment leading to a feasible routing that minimizes expansion and congestion costs. Thus, the problem results in finding a trade-off between investment and routing costs.

Luna & Mahey [19] modeled the capacity expansion and flow assignment as a piecewise convex multicommodity flow problem. Congestion is modeled by a convex increasing function for a given capacity and, at given breakpoints, which represents the maximum tolerable congestion for the users (thus strictly lower than the available capacity), expansion to a higher capacity is decided, decreasing the marginal congestion cost in a discontinuous way. The combinatorial nature of the problem, related to arc expansion decisions, is therefore embedded in a continuous objective function that encompasses congestion and investment costs. The resulting objective function is continuous, but it is nonconvex and nonsmooth. Mahey and Souza [22] derived local optimality conditions for the model proposed in [19]. By exploiting complete optimality conditions for local minima, Souza *et al. *[25] give the convergence analysis of the negative-cost cycle canceling method. Remark that these former works only consider the case of simple expansions from one installed capacity to a new one. The case of the general expansion problem where several capacity expansion values are available for each arc is analyzed by Ferreira & Luna [7], where a method to find solutions with performance guarantee is introduced.

The contribution of this paper is to provide an algorithm to assure the global optimality of the problem, also including computational experiments to show the efficiency of the method to cope with the case of simple expansions. After showing the correspondences between the continuous model of Luna & Mahey and a mixed-integer model, we analyze the numerical behavior of an implicit enumeration algorithm. The discrete model is used to assign capacities to a subset of arcs and hence to define what is called a partial solution. The continuous model is used to compute, given a partial solution, a lower bound on the value of the best solution that can be obtained with the assigned capacities. The resulting procedure is tested on different types of networks and shown to be quite efficient to solve the mixed-integer nonlinear model.

**2 RELATIONS BETWEEN THE CONTINUOUS AND MIXED-INTEGER MODELS **

This section presents the network expansion model. The basic component of the model is a digraph G = (*V, E*) with *n *nodes and *m *arcs. Any kind of traffic between a given pair of nodes is treated as a separate commodity *k*. Let *T *be a (*n ×n*) traffic requirement matrix such that *t _{ij} *is the traffic between origin

*i*and destination

*j*. We will consider the problem of deciding which arcs should be expanded from a given installed capacity

*c*0 to a greater capacity

*c*while minimizing the total congestion and expansion costs.

_{l}Given a commodity *k*, we consider the set of directed paths *P _{k} *joining the corresponding origin and destination. Let

*x*be the amount of flow of commodity

_{kp}*k*through the path

*p ∈ P*and

_{k}*a*its arc-path incidence vector defined by

_{kp}The vector *x *is composed by the component *xe *which denotes the total flow on arc *e*, and also by the component *x _{kp} *which denotes the flow of commodity

*k*routed through path

*p*. These two components are related by

The set of multicommodity flow vectors, denoted by M(*T *) can be described by the arc-path formulation, *i.e.*, for each commodity *k *flowing between nodes *i *and *j*, the active paths must satisfy

That implicit formulation (as the paths are not known in advance) is generally preferred to thenode-arc formulation where *x _{e}* = and each

*x*is a flow vector on G satisfying flow constraints for commodity

^{k }*k*.

We assume now that for each arc *e *in the topology is assigned a positive capacity *c _{0e}* that is expandable to a larger capacity chosen among a given set of capacities

*c*< ... <

_{le }*c*at given fixed costs π

_{Ne}*= 1,...,*

_{le}, l*N*. Let δ

*= 1,...,*

_{le}=c_{le}-c_{0e}, l*N*, be the increment of capacity to the

*l*-th capacity value. The capacity expansion model will minimize the total congestion cost plus the expansion fixed costs. Let Φ(

*c*) be the arc congestion function for a given capacity

_{e}, x_{e}*c*. It is assumed that Φ is expressed in terms of monetary values and that it is convex smooth and increasing up to infinity on the interval [0,

_{e}*c*]. A common choice is the Kleinrock's average delay function valid for M/M/1 queues which is proportional to, see [3, 14].

_{e}A mixed integer model makes use of a binary variable *y _{le}, l *= 1,...,

*N, e ∈E*, that assumes 1 if capacity of arc

*e*is to be expanded from

*c*to

_{0e}*c*and 0 otherwise. We can now define a mixedinteger nonlinear model for the capacity expansion problem (DCE):

_{le}We will now study the relationship between (DCE) and a continuous model which does not makeuse of any boolean decision variables *y *(CCE):

where we assume π* _{0e} *= 0, ∀

*e ∈ E*, and

*f*(

_{e}*x*) = min{Φ(

_{e}*c*) + π

_{le}, x_{e}*= 0,...,*

_{le}, l*N*}.

**Remarks. **

1. Thanks to the feasibility assumption above and the fact that Φ(

c) → +∞ whenever_{Ne}, x_{e}x↑_{e}c, we do not need any capacity constraint in the continuous model._{Ne}2. As shown on Figure 1, where the nonconvex resulting arc cost function of (CCE) is represented by a bold line, we denote by γ

_{(l-1)e},l= 1,...,Nthe breakpoint at which expansion occurs fromc_{(l-1)e }toc. The breakpoint can thus be interpreted as the capacity where congestion is such that the network manager is willing to pay for a new expansion. Thus, π_{le}- π(_{le }l-1)= Φ(_{e }c_{(l-1)}e, γ_{(l-1)e}) - Φ(c, γ_{le}_{(l-1)e)}is the new expansion cost converted in congestion cost units.3. The arc cost function in (CCE) is continuous but nonconvex and nonsmooth at the break-points γ

. It is shown in [19] how one can easily compute a lower bound on the optimal value of (CCE) by taking the convex envelope of each arc cost function._{le}

**Proposition 1. ***If *(*x, y*) *is feasible for (DCE), x is feasible for (CCE); If x is feasible for (CCE), then there exists y such that *(*x, y*) *is feasible for (DCE); If one of both problems is infeasible, so is the other one. *

The proof is straightforward and it is omitted here for sake of simplicity, as the correspondence between (DCE) and (CCE) works with feasibility the same way it works with optimality, as developed below.

The following lemma is a direct consequence of the cost structure of (DCE).

**Lemma 1.** *Let (x*, y*) be an optimal solution of (DCE); then, we have the correspondences:*

*Moreover, if there exists an arc e with *= γ_{(l-1)e}, *then either y**_{(l-1)e} = 1 *and* = 0 *or y**_{(l-1)e} = 0 *and* = 1, *so the optimal solution is not unique. *

The two cases where is not a breakpoint are straightforward. If = γ(_{l-1)e}, we have:

Φ(*c*_{(l-1)e}, γ_{(l-1)e)} + π_{(l-1)e} = Φ(*c _{le}*, γ

_{(l-1)e)}+ π

_{le }

which shows that the value of the arc cost function does not change whenever *y**_{(l-1)e} = 1 and = 0 or, conversely, *y**_{(l-1)e} = 0 and = 1. The correspondence between optimal solutions of (DCE) and (CCE) follows immediately.

**Proposition 2. **i) *If *(*x** , *y** ) *is an optimal solution of (DCE), then x* is optimal for (CCE) and the cost values are equal. *

ii) *If x* is an optimal solution of (CCE), then *(*x** , *y** ) *is optimal for (DCE) with: *

*where l *= 1,..., *N, and the cost values are equal.*

Finally, we would like to point out that the tight relationship between the optimal solutions of both models does not mean that they are equivalent. In general, the continuous model is not able to take into account additional constraints on the topology which, unlike, can be generally done by the *y*-variables. Nevertheless, we will mention a few common situations where it is possible to convert such constraints from (DCE) to (CCE):

a. Many models of network design require the same capacity on arcs (i, j) and (j, i) between two adjacent nodesiandjin G. The orientation of arcs (i, j) and (j, i) being mainly to model the flow that pass fromitojand fromjtoiin a common physical link. As these two flows actually share a common link, it is required in such models symmetry between capacities of arcs (i, j) and (j, i). This is modelled in (DCE) by the constraintyfor some arc_{ij}= y_{ji}e= (i, j). To obtain the same effect, we must add the following constraints for eachlin (CCE):

(*x _{ij} *- γ

*)(*

_{lij }*x*- γ

_{ji}*)*

_{lji }__>__0

b.

Cutset constraints: LetAbe a subset of nodes ofVandCthe corresponding cutset. If subset_{A}Acontains somehow crucial nodes for the network, the arcs in the cutset,i.e., those with one extremity inAand the other inV\A, may be considered bottleneck arcs as they are the only to carry flow between nodes inAand nodes inV\A. Thus, it may be of interest forcing the subsetAto be connected to the other nodes by at least one expanded arc. This is modelled in (DCE) by the constraint Σe ∈ C_{A}y_{le}>1, which is equivalent in(CCE) to:

Observe that both constraints derived in a. and b. define polyhedral nonconvex regions of * ^{m }*.Such a situation could not be treated in the following approach, that requires convexity to assureglobal optimality.

**3 GLOBAL OPTIMIZATION STRATEGY **

In this section, we follow the exposition provided by Geoffrion [13] to describe the application of an implicit enumeration algorithm to solve the capacity expansion and flow assignment problem. Geoffrion [13] presents, including his own developments, an unified approach of the works of Balas [1] and Glover [16] on the implicit enumeration scheme. See also Balas [2]. We consider that each arc *e *is expandable from an installed capacity *c*0*e *to a capacity *c*_{1e}, although the procedure can be generalized to deal with more than one possibility of expansion.

The implicit enumeration algorithm combines information from both discrete and continuous models. A partial solution *S *defines the capacities of a subset of arcs. Here, the discrete model is used to assign capacities to the arcs of . According to the notational convention introducedin [13], for each discrete variable associated to an arc in , the symbol *e *(resp. -*e*) denotes *y _{e} *= 1 (resp.

*y*= 0). The discrete variables associated to arcs not in are called free. As an example, suppose a small network with five arcs and a partial solution

_{e}*S*={1, 3, -5}. In this example,

^{t}*y*

_{1}= 1,

*y*

_{3}= 1,

*y*

_{5}= 0, and

*y*

_{2}and

*y*

_{4}are free. A completion of a partial solution

*S*is defined as a solution that is determined by

*S*together with a binary specification of the values ofthe free variables. It is said to be a feasible completion if the assignment of values to the binary variables leads to a feasible solution. The four possible completions for

*S*in the above example are {1, 3, -5, 2, 4} where the free variables assume

*y*

_{2}= 1 and

*y*

_{4}= 1; {1, 3, -5, -2, 4} where

*y*

_{2}= 0 and

*y*

_{4}= 1; {1, 3, -5, 2, -4} where

*y*

_{2}= 1 and

*y*

_{4}= 0; and {1, 3, -5, -2, -4} where

*y*

_{2}= 0 and

*y*

_{4}= 0.

A key feature of implicit enumeration is the ability to generate information that can be used to exclude all the completions of a partial solution *S *from further consideration. Here, the continuous model is used either to provide a lower bound on the value of the best feasible completion of *S, i.e.*, a feasible completion that minimizes the objective function among all feasible completions of *S*, or to show that *S *has no feasible completion. To do this, we solve a convex multicommodity flow problem *PS*:

where * _{e}, e ∈* , is fixed at partial solution

*S*and

*con*v(

*f*(

_{e}*x*)) is the convex envelope of function

_{e}*f*(

_{e}*x*),

_{e}*c.f.*, Figure 1. Any efficient algorithm designed for convex multicommodity flow problems (see [24] for instance) can be employed to solve

*P*. The value

_{S}*of the optimum solution of*

*P*is a lower bound on the value of the best completion of

_{S}*S*.

An upper bound is always possible to be derived if *P _{S} *has an optimum solution. Let be an optimum solution of

*P*, then

_{S}is an upper bound. Remark that for feasible problems, an upper bound is derived even if *S *= Ø. We store the best upper bound found during the search as the incumbent. If a partial solution *S *leads to an upper bound that improves upon the incumbent, then it replaces the latter as the new incumbent. If the lower bound is greater than or equal to or *P _{S} *is infeasible, then

*S*is fathomed. And in this case all completions of

*S*have been implicited enumerated as they can be excluded from further consideration.

A partial solution is said to be nonredundant if it cannot generate a completion equal to one generated with a previous solution that was fathomed. Geoffrion [13] gives a procedure for generating nonredundant partial solutions that terminates fathoming all feasible solutions. By starting with = Ø, *i.e., S*^{0 }has no capacity previously assigned to any arc, *P _{S}*0 gives a lower bound

_{0}to the optimum value and a first upper bound

*0 taken as the incumbent . Then, the procedure augment the partial solution by assigning a capacity to an arc at a time. Now suppose a partial solution*

_{S}*S*is fathomed. Nonredundancy is achieved by having at least one element of subsequent partial solutions complementary to

^{t}*S*. The last element that was added to

^{t}*S*

^{t-1 }to generate

*S*is then underlined and changed to its complement. In the above example, the sequence would be

^{t}*S*

^{1 }={1},

*S*

^{2 }={1, 3},

*S*

^{3 }={1, 3, -5}. If

*S*

^{3 }could be fathomed, then the next partial solution in the sequence would be

*S*

^{4 }={1, 3,

__5__}, which means change capacity assigned to arc 5 from

*c*

_{0,5}in

*S*

^{3}(

*y*

_{5}= 0) to

*c*

_{1,5}in

*S*

^{4}(

*y*

_{5}= 1). If

*S*

^{4 }could also be fathomed, then

*S*

^{5 }would be

*S*

^{5 }={1,

__-3__}. Otherwise,

*S*

^{5 }would be generated by assigning a value to a free variable and by adding it to

*S*

^{4}, for instance

*S*

^{5 }={1, 3,

__5__, -2}. After fathom a solution

*S*, the procedure locates the rightmost element of

^{t}*S*that is not underlined. If there is none, then end with the stored incumbent as being the optimal solution. Otherwise, replace this element by its complement underlined

^{t}*e*→

__-__and delete all elements that are on the right. In the example, if

*e**S*

^{5 }={1, 3,

__5__, -2} and also

*S*

^{6 }={1, 3,

__5__,

__2__} could be both fathomed, the next solution generated would be

*S*

^{7 }={1,

__-3__}.

We now discuss some strategies to choose an arc *e ∈ *\ *E *with which augment a partial solution *S ^{t}* that cannot be fathomed. Note we chose an arc among those such that

*con*v(

*f*(

_{e}*)) <*

_{e}*f*(

_{e}*).Such strategies rely on the flow distribution of an optimum solution of*

_{e}*P*. Various strategies have been tested:

_{St}• assign

cto the arc_{1e}ewith the highest valuef(^{'+}_{e}) of the right partial derivative of arc_{e}ecost functionfwith respect to_{e}._{e}• assign

cto the arc_{0e}ewith the lowest valuef(^{'-}_{e}_{e}) of the left partial derivative of arc costefunctionfewith respect to._{e}• the arc inducing the highest value of when

Pis defined for ∪{_{S}e} assigningcto_{0e}eif_{e }<γand_{e }cotherwise;_{1e}• the arc with the highest flow value;

The last rule was found to be the most efficient, and was adopted.

We employ algorithms proposed in [7], which have guaranteed performance, to find a good initial upper bound _{0}. Those algorithms find a feasible solution and then gradually reduce the objective value of the obtained solution until no better solution can be found. They are based in two phases: a common first phase where a multicommodity flow problem taking the convex envelope of arc cost functions is solved and a lower bound is found; and a second phase where the obtained routing is used as starting point for switching methods between capacity assignment and the application of a local search algorithm until no more improvement occurs.

The pseudo-code of the implicit enumeration algorithm is as follows:

**Implicit enumeration procedure**

Step 1 - | Initialize with S^{0 }=Ø. Set t = 0 and = _{0}. |

Step 2 - | Solve z*to the best completion of _{S}t S. If ^{t}Pis feasible, then obtain the respective upper bound _{S}t and check to update . _{S}t |

Step 3 - | If it is possible to fathom Sis infeasible or ^{t}, i.e., P_{S}t z*_{S}t > , then go to Step 5. Else go to Step 4. |

Step 4 - | Augment Sby adding a free variable ^{t }ye with e or -e to obtain S. Set ^{t+1}t =t + 1 and return to Step 2. |

Step 5 - | Locate the rightmost element of e →-and delete all elements to the right. Set e t =t + 1 and return to Step 2. |

**4 NUMERICAL TESTS **

Numerical tests were performed to analyze the numerical behavior of the proposed algorithm and the influence of the parameters on its performance. Three different network topologies were used for the computational tests: the C-NET introduced in [24]; the RING introduced in [11]; and the NTS100 generated using a special program driver [6, 7]. For the sake of illustration, Figure 2 and 3 present the topologies of networks C-NET and NTS100, respectively. Table 1 gives the main characteristics of the three networks used in our numerical experiments. Given a root node* r*, the hop-depth of a node *i ∈ V *\{*r*} is the number of arcs in the path between *r *and *i *that has minimum length. The hop-based diameter of the graph is the largest hop-depth among the nodes of the network. C-NET and RING are full duplex networks, *i.e.*, flow can be sent between two nodes in both directions and simultaneously.

Two sets of tests were made with these topologies:

- The first set concerns the C-NET network [24]. The aim is to verify the influence of thetraffic throughput increase and the capacity expansion factor (

c_{1e}/c_{0e}) on the number of iterations performed by the implicit enumeration algorithm. Results are compared with the ones obtained solving the mixed integer nonlinear formulation of the problem (DCE) using two commercial solvers: BARON and LINDOGlobal [5, 17, 27] and modeling the problem with GAMS [10]. The problems were solved with help of the NEOS Server [5, 17];- The second test set was performed on the other two topologies RING and NTS100. The aim is to assess the effectiveness of the proposed implicit enumeration algorithm to solve the capacity expansion problem with different scenarios of larger networks and heterogeneous traffic requirement demands.

The computation of the lower bound was performed by a specialized Flow Deviation algorithm for the convex multicommodity flow problem. That algorithm has been shown to be efficient since the early work of Fratta *et al. *[9], but it is generally known to become very slow when the algorithm approaches the optimal solution (see [3] for instance). To accelerate its convergence, a parallel tangent procedure (PARTAN, see [8]) was introduced in the direction finding step. The algorithm was coded in C.

The following notations were used:

•

Φ* is the global optimal value;• optimal value of the convexified problem;

•

Φis the solution obtained by the performance guaranteed algorithms [7];_{pg }•

Φis the solution obtained with BARON [27];_{Baron }•

Φis the solution obtained with LINDOGlobal [5, 17];_{LG }•

Nis the number of partial solutions enumerated to guarantee the global optimality;_{ggo}•

Nis the number of partial solutions enumerated until the optimal solution is found be not necessarily proven to be optimal;_{o}•

Nis the number of visited nodes made by the Baron solver;_{Baron}•

Nis the number of visited nodes made by the LINDOGlobal solver;_{LG}• α = ;

• α

= ._{pg }

**4.1 First test set **

In these experiments, the congestion cost is derived from the Kleinrock's average delay function for M/M/1 queueing networks and defined by ρ . We fix ρ = 1, see [11]. All arcs have the same initial capacity. And the expansion cost π_{1e}, for each arc *e ∈ E*, is given, for each value of the parameter γ = equal to 0.7 and 0.9, by

Tables 2 to 5 present the results for the C-NET network varying traffic requirement demands, parameter γ and the available capacity for expansion.

The results show that the algorithm is rather sensitive to the capacity expansion factor. When *c _{1e}* = 4

*c*the number of iterations increase significantly, as observed for

_{0e}*N*and

_{ggo}*N*in Tables 4 and 5 in contrast with

_{o}*N*and

_{ggo}*N*in Tables 2 and 3. Increasing γ affects the efficiency of the algorithm. Both effects can be illustrated by the values of α and α

_{o}*.*

_{pg}As the efficiency of the algorithm depends on the quality of the lower bound to exclude completions of a partial solution from further consideration, the weakest is the convexification bound, the biggest is the number of solutions enumerated. Observe that the global optimal solution was reached in many instances, some of them well before *N _{ggo} *»

*N*.

_{o}Table 6 shows the results obtained solving the mixed integer nonlinear formulation of the problem (DCE) using two commercial solvers: BARON and LINDOGlobal [5, 17] and the resultsof the proposed algorithm. The Baron and the LINDOGlobal solvers were not able to guarantee the optimality of the results within the adopted time limits and or the number of iterations limits. They were not able to obtain the guaranteed global optimum of the C-NET problem. The LINDOGlobal solver uses linear approximations of the original problem which explains the large number of iterations within the limits of runtime.

**4.2 Second set of tests **

The second set of tests concern the RING and the NTS100 topologies with heterogeneous trafficrequirement demands. The congestion costs, as in the precedent set of experiments, is given by ρ . The following scheme was adopted for all these experiments:

- An initial feasible solution is computed using algorithms proposed in [7].

- The initial demand between each OD-pair was fixed to 10 for the RING network and to 3 for the NTS100 network. The installed capacities are given in Table 7.

• In a first set of instances, the demand is uniformly multiplied by a throughput factor of 20%, 40%, 60%, 80% and 100% until capacity expansion becomes economically interesting.

- In the second set, the demand increase is no more uniform: half of the demands receive the same throughput factor, but 25% receive 50% more, and 25% receive 25% less. Theresults are exposed in Tables 8 to 11.

- An initial distribution of the increase of demand between the OD-pairs was defined randomly and fixed for each instance.

- The parameter

ρused to compute congestion costs is set to 500, 1000, 5000, and 10000.

A small number of iterations was needed to solve the proposed instances. Considering the complexity of the studied problem (CNET = 2^{34 }= 17,179,869,184 solutions, RING = 2^{60 }solutions, NTS100 = 2^{187 }solutions), one can evaluate the effectiveness of the proposed method. The second set of test problems was easier to be solved. This can be explained by noting that for these instances, the quality of the lower bound obtained by solving the problem with the convex envelope of the integrated function of congestion and expansion costs, *c.f.*, Figure 1, is very good.

At each iteration of the implicit enumeration algorithm a convex multicommodity flow problem is solved. A convex multicommodity flow problem can be difficult, especially if feasibility problems appear when an arc flow is close to the arc capacity. The decision to apply the Partan was primarily due to the gain in speed in solving problems is adequate for a precision of about 1% or better. Each multicommodity flow subproblem spent from a second fraction to several minutes.

All test problems have many local optima with objective function value close or equal. The implicit enumeration depth search approach has modest memory requirements. It needs to store only the path from the root node to the leaf node in the search tree. As the problem has many solutions, the depth search has a good chance of finding an optimal solution after exploring a small portion of the entire search space. The proposed algorithm can be parallelized and we believe that a good scalability can be get.

The proposed algorithm can be further specialized if more information about the problems are considered. For example, in real scenarios budget and operational constraints may restrict the set of candidate arcs for expansion. The proposed procedure can also be used for capacity reduction problems.

**5 CONCLUSION **

We have applied successfully an adapted implicit enumeration scheme to a mixed-integer nonlinear model for the capacity expansion and flow assignment in multicommodity networks. The key fact to explain the good performance of the method is the use of a reliable bounding procedure at each node of the enumeration tree, which is based on the convex hull of the equivalent continuous objective function.

The literature has scarce approaches to find exact solution of the discrete capacity expansion and flow assignment problem. Basically two lines of approaches have been adopted: Benders decomposition [20] and methods based on difference of convex functions - DC programming [21]. Recently a conic quadratic formulation was proposed [26]. The results show that the implicit enumeration algorithm proposed here can be considered as an alternative available to solve to global optimality such large scale problems.

Future research is to be done on the extension of the proposed approach to deal with multiple choices of available capacities for expansion on each arc. It remains to show that lower bounds obtained with the convex envelope of the integrated function of expansion and congestion costs remain sharp in the presence of multiple expansions.

**ACKNOWLEDGMENTS **

This work has been supported by CNPQ (305164/2010-4) and (473763/2011-7).

**REFERENCES **

[1] BALAS E. 1965. An additive algorithm for solving linear programs with zero-one variables. *Operations Research*, **13:**517-546. [ Links ]

[2] BALAS E. 1966. Discrete programming by the filter method. *Operations Research*,**13:**915-955. [ Links ]

[3] BERTSEKAS DP & GALLAGER RG. 1987. Data Networks, Prentice-Hall. [ Links ]

[4] DE CAMARGO RS, MIRANDA JR G, FERREIRA RPM & LUNA HP. 2009. Multiple allocation huband-spoke network design under hub congestion. *Computers and Operations Research*, **36:**3097-3106. [ Links ]

[5] CZYZYK J, MESNIER M & MOR J. 1998. The NEOS Server. *IEEE Journal on Computational Science and Engineering*, vol. 5. [ Links ]

[6] DOAR MB & NEXION A. 1996. A better model for generating test networks. *Proceedings of Globecom '96*, pp. 86-93. [ Links ]

[7] FERREIRA RPM & LUNA HPL. 2003. Discrete capacity and flow assignment algorithms with performance guarantee. *Computer Communications*, **26:**1056-1069. [ Links ]

[8] FLORIAN M, GULAT J & SPIESS H. 1987. An efficient implementation of the PARTAN variant ofthe linear approximation method for the network equilibrium problem. *Networks*, **17:**319-339. [ Links ]

[9] FRATTA L, GERLA M & KLEINROCK L. 1973. The flow deviation method: an approach to storeand-forward communication network design. *Networks*, **3:**97-133. [ Links ]

[10] www.gams.com, consulted on september 19, 2012. [ Links ]

[11] GAVISH B & NEUMANN I. 1989. System for routing and capacity assignment in computer communication networks. *IEEE Transactions on Communications*, **37:**360-366. [ Links ]

[12] GAVISH B & ALTINKEMER K. 1990. Backbone network design tools with economic tradeoffs. *ORSA Journal on Computing*, **2:**236-252. [ Links ]

[13] GEOFFRION AM. 1967. Integer programming by implicit enumeration and Balas' method. *SIAM Review*, **9:**178-190. [ Links ]

[14] GERLA M & KLEINROCK L. 1977. On the topological design of distributed computer networks.*IEEE Transactions on Communications*, **25:**48-60. [ Links ]

[15] GERLA M, MONTEIRO JAS & PAZOS R. 1989. Topology design and bandwith allocation in ATM nets. *IEEE Journal on Selected Areas in Communications*, **7:**1253-1261. [ Links ]

[16] GLOVER F. 1965. A multiphase-dual algorithm for the zero-one integer programming problem. *Operations Research*, **13:**879-919. [ Links ]

[17] GROPP W & MOR E´J. 1997. Optimization Environments and the NEOS Server. *Approximation Theory and Optimization*, M.D. Buhmann and A. Iserles, eds., pp. 167-182, Cambridge University Press. [ Links ]

[18] ISHFAQ R & SOX C. 2012. Design of intermodal logistics networks with hub delays. *European Journal of Operational Research*, **220:**629-641. [ Links ]

[19] LUNA HPL & MAHEY P. 2000. Bounds for global optimization of capacity expansion and flow assignment problems. *Operations Research Letters*, **26:**211-216. [ Links ]

[20] MAHEY P, BENCHAKROUN A & BOYER F. 2001. Capacity and flow assignment of data networks by generalized Benders decomposition. *J. of Global Optimization*, **20:**173-193. [ Links ]

[21] MAHEY P, PHONG TQ & LUNA HPL. 2001. Separable Convexification and DC Techniques forCapacity and Flow Assignment Problems. *RAIRO Operations Research*, **35:**269-281. [ Links ]

[22] MAHEY P & SOUZA MC. 2007. Local optimality conditions for multicommodity flow problemswith separable piecewise convex costs. *Operations Research Letters*, **35:**221-226. [ Links ]

[23] MORABITO R & DE SOUZA MC. 2010. Roteamento de multi-fluxos em redes de filas genericas. *Pesquisa Operacional*, **30:**583-600. [ Links ]

[24] OUOROU A, MAHEY P & VIAL JP. 2000. A survey of algorithms for convex multicommodity flow problems. *Management Science*, **46:**126-147. [ Links ]

[25] SOUZA MC, MAHEY P & GENDRON P. 2008. Cycle-based algorithms for multicommodity networkflow problems with separable piecewise convex costs. *Networks*, **51:**133-141. [ Links ]

[26] SINAN G. 2011. A conic quadratic formulation for a class of convex congestion functions in networkflow problems. *European Journal of Operational Research*, **211:**252-262. [ Links ]

[27] TAWARMALANI M & SAHINIDIS NV. 2005. A polyhedral branch-and-cut approach to global optimization. *Mathematical Programming*, **103:**225-249. [ Links ]

Received December 9, 2011

Accepted January 22, 2013

* Corresponding author