NANZAN-TR-2002-02

On the hub-and-spoke model with arc capacity constraints

Mihiro Sasaki and Masao Fukushima

September 2002

Technical Report of the Nanzan Academic Society
Mathematical Sciences and Information EngineeringOn the hub-and-spoke model with arc capacity
constraints yMihiro Sasaki z
Department of Mathematical Sciences
Faculty of Mathematical Sciences and Information Engineering
Nanzan University
27 Seirei, Seto, Aichi, 489-0863, Japan
Tel: +81-561-89-2000 Fax: +81-561-89-2014
Email: mihiro@ms.nanzan-u.ac.jp
Masao Fukushima
Graduate School of Informatics
Kyoto University
Kyoto, 606-8501, Japan
Email: fuku@amp.i.kyoto-u.ac.jp

Abstract

In this paper, we present a new formulation of one-stop capacitated hub-and-spoke
model as a natural extension of the uncapacitated one-stop model. The model involves
arc capacity constraints as well as hub capacity constraints, which enables us to incorporate some practical factors into the model. We also present a branch-and-boundbased exact solution method with Lagrangian relaxation bounding strategy, and report
some results of numerical experiments using real aviation data. Computational results
show the practical signi cance of the proposed capacitated model.

yz

The work was supported by Nanzan University Pache Subsidy I-A.
Corresponding author11 Introduction
The hub location problem has attracted much attention since O'Kelly [17] formulated
the single allocation hub-and-spoke model as a quadratic integer program. Campbell [4]classi es the p-hub location problem into four classes, and presented basic formulations andformulations with ow thresholds for spokes for each of them. The four classes are the p-hub
median problem, the p-hub center problem, the uncapacitated hub location problem and the
hub covering problem. Besides the basic single allocation and multiple allocation models,
Campbell [4] presents models involving ow thresholds for spokes, and points out that the phub median problem and the uncapacitated hub location problem have mainly been studied
and other classes of problems have not been considered seriously. Skorin-Kapov, SkorinKapov and O'Kelly [23] considered the uncapacitated p-hub median problem and developed
linear programming formulations of both single and multiple allocation models that yieldtight relaxations to the original problem. For the same class of problems, O'Kelly et al. [18]introduced a more compact formulation and discussed sensitivity analysis of the model withrespect to various parameters. On the other hand, Aykin [1] formulated the capacitatedhub-and-spoke model as a 0-1 integer programming problem and proposed two heuristics tosolve the problem. Ebery et al. [11] also addressed the capacitated multiple allocation hub
location problem that arises in a postal delivery system.
The above mentioned models are of two-stop type in the sense that a trip between eachorigin-destination (OD) pair uses at most two hubs. Therefore, the number of variables andconstraints in the problem rapidly grows with the size of the model and hence we may expect
to solve problems of practical size only approximately. Sasaki, Suzuki and Drezner [21, 22]considered the uncapacitated one-stop model, in which a trip between any OD pair uses only
one hub. The one-stop model appears to be useful in certain situations and is simple enoughto be solved exactly. In this paper, we propose a capacitated one-stop model as a naturalextension of the uncapacitated one-stop model of [21, 22] and present a branch-and-boundbased exact solution method with Lagrangian relaxation bounding strategy. Speci cally,
we introduce capacity constraints on both arcs and hubs. Similar arc capacitated network
design problem have been considered in [13, 14]. However, studies on arc capacitated hub
location problems are scarce. Hub capacity may be regarded as the maximum number ofpassengers a node can deal with when it is selected as a hub. For example, the number of
runways and the number of slots assigned to the airline company may constitute the hub
capacity. On the other hand, arc capacity may represent the number of available aircrafts
for the airline company on that arc. By introducing arc capacity constraints as well as hub
capacity constraints, we can enrich the model so as to take into account more practicalsituations. Bryan [2] rst incorporated arc capacities into a hub location model. The model2

is an extension of the model proposed by O'Kelly et al. [19], which uses a piecewise linear
cost function. Assuming hub locations were given, they only focused on the allocation phaseand examined how the arc capacity levels a ect the total network cost.
This paper is organized as follows. In Section 2, we formulate the capacitated one-stop
hub-and-spoke model. In Section 3, we present a branch-and-bound method for solving the
proposed model and particularly describe Lagrangian relaxation strategy for nding goodlower bounds for partial problems. In Section 4, we report some numerical experience. In
Section 5, we give concluding remarks and brie y mention some future work.

2 Capacitated One-Stop Model
2.1 Model Description

Jaillet et al. [15] propose an airline network design problem. They examined the resultant
network structure and gured out that the one-stop policy could bring results as good as
those of two-stop policy, especially in the situation with relatively high level of demands.
In some countries like Japan, a majority of domestic routes are provided by one-stop services via a hub. Even in the United States, a large number of trips seem to use one-stoproutes. Motivated by these observations, an uncapacitated one-stop model has been proposed in [21, 22]. In this paper, we generalize this model to construct a capacitated one-stop
model. We formulate the problem as a 0-1 mixed integer programming problem and proposea branch-and-bound method with Lagrangian relaxation bounding strategy.
Aykin [1] formulated the capacitated two-stop hub-and-spoke model and proposed two
heuristics for the problem. In the model, capacity constraints on hubs are primarily considered. However, since airline companies may sometimes fail to provide ight service to all
passengers who want to use some arc, it is also practically signi cant to consider capacity
constraints on arcs as well as on hubs. Concerning route patterns, we suppose that nonstop
services between non-hub nodes are not available, while passengers whose origin or destination is a hub are permitted to travel via another hub. Possible seven route patterns according
to this rule are shown in Figure 1.
One of the major factors that characterize a hub location model is its allocation rule.
Basically, there are two possible choices of allocation rules; the single allocation rule andthe multiple allocation rule [3]. Allocating each demand node to more than one hubs is
prohibited in the former rule and allowed in the latter rule. The rst hub-and-spoke model
presented by O'Kelly [17] adopts the single allocation rule. Recent papers often adopt themultiple allocation rule because it allows more exible routing between OD pairs [5, 8, 18],
while some others adopt single allocation rule [7, 10]. Skorin-Kapov et al. [23] and Ernst et3

hub node

non-hub node

Figure 1: Route Patterns
al. [9] employ both rules in the uncapacitated p-hub median problem. These models are all
two-stop models. The choice of the allocation rule makes a sensible di erence in a two-stop
model. However, in a one-stop model like the presented model, it is necessary to adopt themultiple allocation rule. Because, if we adopt the single allocation rule in a one-stop model,
then any two demand nodes allocated to di erent hubs are not reachable.
The routing rule is another important factor in the model. There are also two possible
choices of routing rules; the single routing rule and the multiple routing rule. Only one routeservice provided for each OD pair in the former rule and more than one route services areallowed in the latter rule. In other words, all passengers whose origin and destination are
the same have to travel using the same route in the single routing rule. We note that themultiple routing rule is meaningful only in a capacitated problem, because in an uncapacitated problem, all passengers for each OD pair will use the least cost route even if multiple
routes are available.
In this paper, we adopt the multiple routing rule. The objective of the model is to nd
hub locations and passenger routing for each OD pair that minimize the total transportationand location costs subject to the constraints derived from the assumed rules.

2.2 Model Formulation
We now present a mathematical formulation of the hub location problem with the multiple
allocation and multiple routing rules. The following notation is employed:

N : the set of demand nodes, jN j = n.
A : the set of arcs.
4

H  N : the set of hub candidates, jH j = h  n:
 : the set of OD pairs  = hi; j i; i 2 N; j 2 N; i 6= j . jj = n(n ; 1).
Ai   : the set of OD pairs whose origin is i 2 N .
Bj   : the set of OD pairs whose destination is j 2 N .
p : the number of hubs to be selected.
d : the trip demand (the number of passengers) for OD pair  2 .
ck : the travel cost per passenger between OD pair  2  via hub k 2 H .
fk : the xed cost incurred by selecting hub k.
w : weight parameter for the xed cost.
ak : capacity of hub k 2 H .
bij : capacity of arc (i; j ) 2 A.
Decision variables in the model are as follows:xk : the number of passengers who travel between OD pair  2  via hub k 2 H .
yk : binary variable such that yk = 1 if node k is selected as a hub, and 0 otherwise.
We note that hi; j i denotes the OD pair whose origin is node i and destination is node j , while(i; j ) denotes the arc connecting node i and node j . For an OD pair  = hi; j i, if i 2 H , thenxi is the number of passengers traveling directly from hub i to node j . Similarly, if j 2 H ,then xj is the number of passengers traveling directly from node i to hub j . Moreover, if
i 2 H and j 2 H hold simultaneously, then both xi and xj are the number of passengers
from i to j . In this case, the total number of passengers on the route consisting of the single
arc (i; j ) is given by xi + xj .
Now we are ready to formulate the model as the following mixed 0{1 integer programming
problem.

P0:

XX

X
ck xk + fk yk
2 k2HX
k2H
subject toxk = d ;
 2 ;k2
H
X
X
xj + xi  bij ; i; j 2 H; i 6= j;
2Ai
X 2Bjxj  bij ;i 62 H; j 2 H; i 6= j;2AiXxi  bij ;i 2 H; j 62 H; i 6= j;2BjXyk = p;minimize

k 2H

5

(1)
(2)
(3)
(4)
(5)
(6)

X

xk  ak yk ;xk  0;yk 2 f0; 1g;

2

k 2 H;
 2 ; k 2 H;
k 2 H:(7)

The objective function (1) is the sum of the total travel costs and xed costs. Constraints(2) imply that each passenger between OD pair  travels via one of the hubs. Constraints(3){(5) represent arc capacity constraints; (3) applies to the case where both end points ofarc (i; j ) belong to the hub candidate set, while (4) and (5) apply to the case where oneof the end points belongs to the hub candidate set. Constraint (6) requires that exactly phubs have to be selected. Constraints (7) represent capacity constraints for hubs, implying
that, when hub k is selected, the number of passengers who use hub k cannot exceed the
hub capacity ak . Constraints (7) also assure that any node cannot be connected to non-hub
nodes. Namely, xk for all  2  are forced to vanish unless yk = 1.

3 Branch-and-Bound Method
3.1 Lagrangian Relaxation

We note that P0 can naturally be decomposed into two problems by xing an arbitraryyk at 0 or 1. We call such variable yk and the generated problems a branching variableand partial problems, respectively. A partial problem P may have one of the following two
properties:
1. An optimal solution of P is obtained.
2. It is detected that P cannot produce an optimal solution of P0 .
In either case, we can fathom P , because it is not necessary to decompose P further. To
check whether or not P has property 2, a lower bound test is applied to P . The branch-andbound method repeatedly applies decomposition and lower bound tests until there remainno partial problems that are active, i.e., generated but neither fathomed nor decomposed.
To describe a partial problem P , we partition the hub candidate set H as H = HF [ H1 [
H0 , where HF ; H1 and H0 are speci ed as follows.

HF : the set of indices k such that yk is a free variable.
H1 : the set of indices k such that yk is xed as yk = 1.
H0 : the set of indices k such that yk is xed as yk = 0.
For the partial problem P , either jH1j = p or jH1j < p holds. If jH1j = p, then P becomesa linear programming problem by xing all free variables to be zero. We can then obtain an
6upper bound of P0 by solving this problem, since its optimal solution is feasible to P0 . If thisupper bound is smaller than an incumbent value, we let it be a new incumbent value of P0.
If jH1j < p, then we attempt to obtain a good lower bound for the optimal value of P .
Let us denote H+ = HF [ H1 . Note that we eliminate variables xk ; yk ; k 2 H0 from Pbecause it immediately implies that xk = 0 for all  2  if k 2 H0 . Then partial problemP can be written as follows.

XXck xk +
fk +
fk yk
1
F
2 k2H+k2H
X k2H
subject toxk = d ;
 2 ;+
X k2H Xxj + xi  bij ;i; j 2 H+ ; i 6= j;
2Ai
X 2Bjxj  bij ;i 62 H+ ; j 2 H+;
2Ai
X
xi  bij ;i 2 H+ ; j 62 H+;X2Bj
yk = p ; jH1 j;
k2HFXxk ; ak yk  0;
k 2 HF ;2Xxk  ak ;k 2 H1 ;P : minimize

X X

 2

xk  0;yk 2 f0; 1g;

(8)
(9)
(10)
(11)

 2 ; k 2 H+ ;k 2 HF :We consider relaxing all the arc capacity constraints and the hub capacity constraints(11) by bringing them into the objective function. To achieve this, we introduce Lagrange
multipliers  = (ij ) and  = (k ) corresponding to the arc capacity constraints (8){(10) andthe hub capacity constraints (11), respectively. Then we obtain the following Lagrangian
relaxation problem RP [; ] for partial problem P .
RP [; ] :minimizeX X

Xck xk +fk yk
2 k2H+
k2HF
X
X
X X
ij ( xj + xi )
+
2Ai
2B
H+ j 2H+ ;j 6=i
X
X X j XX i2X
ij xi+ij xj +i2H+ j 62H+ 2Bj
i62H+ j 2H+ 2Ai
X X
k ( xk ; ak yk ) + L 0()
+
k2HF

2

7

X

subject to

k2H+

X

 2 ;xk = d ;yk = p ; jH1j;k2HFXxk  ak ;k 2 H1 ;2

xk  0;yk 2 f0; 1g;

whereL 0() =

X
k2H1

fk ;

 2 ; k 2 H+;k 2 HF ;X X
i62H0 j 62H0 ;j 6=i

ij bij :It is well-known that the optimal value of RP [; ] provides a lower bound of the optimalvalue of P for arbitrary   0 and   0. We can obtain an optimal solution of RP [; ]easily by decomposing it into some simple problems. First, note that the objective functionof RP [; ] can be rewritten as

X X

2 k2H+

+

X

k2HF

ck xk +X X

2 k2HFk xk +(fk ; k ak )yk + L 0():X X

i2H j 2H+ ;j 6=i

ij

X

2Ai

xj +

X X

i2H+ j 2H;j 6=i

ij

X

2Bj

xi
(12)

Let  and  denote the origin and the destination of OD pair , respectively. Then thethird term of (12) can be rewritten asXX2 j 2H+ ;j 6=   j xj :(13)

Similarly, the fourth term of (12) can be rewritten as

XX2 i2H+ ;i6= i  xi :(14)

By substituting (13) and (14) into (12), we may rewrite RP [; ] as follows.
RP [; ] :minimizesubject to

0
X@ X

X

X

X

ck xk +  k xk +k  xk +++2 k2H+
k2HF
k2H ;k6= 
k2H ;k6= X(fk ; k ak )yk + L 0()
+
k2HF
X

xk = d ;

 2 ;xk  ak ;k 2 H1 ;

k2H+

X

2

8

1
k xk A

X

yk = p ; jH1j;

k2HFxk  0;

 2 ; k 2 H+;k 2 HF :yk 2 f0; 1g;

Let us de ne8Ck to represent the aggregate coecient of xk :>>ck +   k + k  + k , if k 2 HF ; k 6=  ; k 6= >>>
if k 2 HF ; k = >< ck +   k + k ,Ck = > ck + k  + k ,if k 2 HF ; k = >>ck +   k ,if k 2 H1 ; k = >>>: ck + k  ,if k 2 H1 ; k =  .
Then, problem RP [; ] can be decomposed into the linear programming problem:
minimizesubject to

X X

2 k2H+

X

k2H+

Ck xkxk = d ;  2 ;P

2 xk  ak ;xk  0;k 2 H1 ;k 2 H+ ;(15)

and the 0{1 integer programming problem:

Xminimizesubject to

k2HF

X

k2H

F


(fk ; k ak )yk

yk = p ; jH1j;yk 2 f0; 1g;(16)

k 2 HF :Note that the linear programming problem (15) reduces to a transportation problem, if weregard the nodes k 2 H+ as source nodes and OD pairs  2  as sink nodes. Moreover, an
optimal solution yk ; k 2 HF , of the 0{1 integer programming problem (16) is obtained in thefollowing simple way: Order the indices k 2 HF as k1 ; k2 ;    ; kjHF j in the increasing orderof the coecients fk ; k ak , and let

8>< 1; if k 2 K;yk = >
: 0; if k 62 K;

where K = fk1 ; k2 ;    ; kp;jH1 jg.
9(17)

3.2 Subgradient Method
It is practically important to use appropriate values of Lagrange multipliers so as to obtaina good lower bound for the optimal value of the partial problem P . The optimal Lagrange
multipliers can be obtained by solving the following Lagrangian dual problem for P .

LD : maximize G (; )subject to ;   0;
where G (; ) denotes the optimal value of RP [; ]. This is a problem of maximizing a
nondi erentiable concave function. Since it is expensive to obtain an optimal solution ofLD exactly, we use the subgradient method, which is practically useful because it can nd
an approximate solution of LD conveniently.
The subgradient method for solving LD consists of the following iterative process: Choose
initial Lagrange multipliers (0)  0; (0)  0, and then successively update ((l) ; (l) ) by(l+1) = maxf0; (l) + (l) g(l) g;
(l+1) = maxf0; (l) + (l) g(l) g;(18)

where (l) > 0 is a stepsize and (g(l) ; g(l) ) is a subgradient of G at ((l) ; (l) ). Using an
optimal solution (x(l) ; y(l) ) of RP [(l) ; (l) ], where x(l) is an optimal solution of the linear
programming problem (15) with (; ) = ((l) ; (l) ) and y(l) is determined by (17), we may
compute each element of the subgradient (g(l) ; g(l) ) as

8 X (l) X (l)>xj + xi ; bij ; if i 2 H+ and j 2 H+;
>
>
2Bj
2A
>
>
< Xi (l)h (l) ixj ; bij ;if i 62 H+ and j 2 H+;g ij = >
2Ai
>
X (l)
>
>
if i 2 H+ and j 62 H+;>: 2Bj xi ; bij ;h (l) ig k =(19)X (l)xk ; ak yk(l) ; k 2 HF : 2

Stepsize (l) is determined according to the rule(l) (l)^(20)
(l) =  (G ; (Gl)  ((l) ; 2 )) ;k(g ; g )kwhere G^ is an estimate of the optimal value of LD , and  is a parameter satisfying    2 ;  for some  > 0 [20]. If we use the exact optimal value G of LD instead of G^ , the10sequence generated by (18) converges to an optimal solution of LD . However, the optimalvalue of LD is usually unknown. So we iterate (18) with (20) using an appropriate estimateG^ of the optimal value G . By using an appropriate adjusting strategy, we may expect thatG^ approaches G and ((l) ; (l) ) converges to an optimal solution of LD .
As a matter of fact, the convergence rate of subgradient method signi cantly depends on
the strategy for updating an estimate of the optimal value. For simplicity, we take  = 1 in
the following discussion. Let G^ > G and 0 <  < 1. Then it can be shown that there is apositive integer S such thathi
 ; G^ )   G ((l+1) ; (l+1) ) ; (2G ; G^ ) ;l  S =) Gmax
;
(2
G







(21)

where Gmax
 is the maximum value of LD obtained up to the l-th iteration with the initialpoint ((0) ; (0) ) [16]. Based on this fact, we propose to update the estimated value G^ atevery S iterations. With this strategy, G^ is expected to decrease steadily and converge toG . A procedure to update the estimated value G^ is described as follows:[Updating G^ ]Input:

 2 (0; 1); S 2 f1; 2;   g.

Step 0: LetG^ be an arbitrary upper bound of G . Select an initial point ((0) ; (0) )．Set

Step 1: Set

Gmax := ;1.

l := 0.

l

Step 2: If mod( + 1; S ) 6= 0 then go to Step 3. Otherwise, go to Step 4.

maxmax(l) (l)max :=G ((l) ; (l) ) > Gmax
 , then update G as G := G ( ;  ), and let (l) and max := (l) . Compute (l+1) ; (l+1) according to (18){(20). Set l := l + 1 andStep 3: If

go to Step 2.

^G^ as G^ := G + (1 ;  )2G; 2; G ( ;  ) . Set
(l+1) := max , (l+1) := max , l := l + 1 and go to Step 1.
max

Step 4: Update the estimated value

(0)

(0)

If Gmax
 ; G^ becomes smaller than a predetermined tolerance, then we stop the iterations.
In general, it is impossible to obtain the exact values of  and S . However, it is seen that S
gets smaller as  gets larger because of (21). From some preliminary computations, we havedecided to set  = 0:5 and S = 5 in our computational experiments.
Although we may use an arbitrary upper bound of G as an initial estimate of the optimalvalue of LD in Step 0, it is not always easy to obtain it. In the above procedure, we use
an arbitrary upper bound for P rather than G , because any upper bound for P is greater
11than or equal to G . We note that the uncapacitated one-stop hub location problem can be
regarded as a shortest path problem if any p hubs are xed. Therefore, for an uncapacitatedproblem, we can easily obtain an upper bound by solving the shortest path problem. For a
hub capacitated problem, we can also nd a feasible ow by solving a shortest path problem,
if we select p hubs so that the total hub capacity is greater than the total demand. Thus, we
can obtain an upper bound of a hub capacitated problem by allocating ows appropriately
to each arc. However, it is not necessarily easy to obtain an upper bound for P in the casewhere both hubs and arcs are capacitated. The reason is that an arc capacitated problem
often becomes infeasible even if we select p hubs so as to satisfy the above-mentioned total
capacity condition. In general, any feasible solution of P is also feasible to P0, though theconverse does not always hold. However, if a feasible solution of P0 is still feasible to P , its
objective value also provides an upper bound of P . Therefore we keep a list of any feasible
solutions to P0 obtained during the branch-and-bound computation, and try to nd a feasible
solution of P from the list so as to use it as an initial estimated value of LD . If there are
more than one feasible solutions to P , then we may choose the one with smallest objective
value. On the other hand, if the list contains no such feasible solution for P , then we try toobtain a new upper bound by solving a new shortest path problem that is constructed fromP by adding p ; jH1j hubs arbitrarily chosen form HF . If the new shortest path problem
is feasible, we use its objective value as an initial estimated value and keep the solution forthe subsequent computation. Otherwise, we let G^ be a tentative estimated value, which is
chosen to be greater than the current incumbent value z  of P0 . We note that there is no
guarantee that such a tentative estimated value is an upper bound of G . If G^ < G , thenG ((l) ; (l) ) will converge to a smaller value than G [16]. Nevertheless, we may fathom
the partial problem P if we can nd a ((l) ; (l) ) such that G ((l) ; (l) ) is greater than theincumbent value z . After fathoming a partial problem and backtracking to another partial
problem, some feasible solutions of P0 in the list may become invalid in the sense that they
can never be feasible to any partial problem generated in the subsequent computation. By
removing such invalid data, we can reduce the area of searching for feasible solutions of apartial problem.

3.3 Upper Bound
Upper bounds for the optimal value of problem P0 play an important role in developing abranch-and-bound method. An upper bound can be obtained by nding any feasible solutionof P0, which is not always easy as described in the previous subsection. However, since wecannot expect to fathom any partial problem without knowing an upper bound for P0, we
need to obtain it before getting into a lower bound computation. In addition, whether we12can obtain a good upper bound in the early stage of iterations has a signi cant e ect on theperformance of the branch-and-bound method.
First we note that any partial problem P as well as P0 reduces to an arc capacitatednetwork ow problem if we select p hubs. In general, it is more likely to be able to nd afeasible ow of the arc capacitated network ow problem by selecting hub candidates withlarger capacities. However, selecting such hubs is generally costly because hubs with large
capacities usually incur a high xed cost. On the other hand, to obtain a good lower bound,
it seems e ective to select those hub candidates which are expected to be selected at optimal
solutions. Consequently, we use the following strategy to obtain an upper bound. First we
select p hubs according to an appropriate criterion, then solve the resultant arc capacitatednetwork ow problem. If this problem is feasible, the obtained objective value is an upperbound. Here we introduce the following criterion as a likelihood of each hub candidate to be
selected as a hub at an optimal solution:Xh(k) = afk + ck ;k2k 2 H;

where > 0 is a parameter. The rst term represents the xed cost per unit capacity
multiplied by a positive parameter and the second term represents the total travel costwhen all passengers use the hub k. Note that we do not consider capacity constraints here.
By introducing this criterion, a hub candidate associated with small h(k) is likely to be
selected as a hub at an optimal solution. We sort h(k) in increasing order and express theordered indices of h(k) as 1; 2 ;    ; jH j. If the arc capacitated network ow problem inwhich hubs f1 ; 2;    ; p g are selected is feasible, we use the objective value as an initial
upper bound. Otherwise, we select p hubs with p largest ak 's and solve the resultant network
ow problem. If this problem is also infeasible, we proceed to branch-and-bound iterations
without knowing any upper bound at this stage.

3.4 Branching Strategies
The computational eciency of branch-and-bound method signi cantly depends on abranching strategy to select a branching variable in each stage to decompose the current
partial problem into new partial problems. Selecting a variable whose value is likely to be
1 at an optimal solution is one of the practical strategies. However, when we proceed tobranch-and-bound iterations without nding an upper bound for P0, we need to obtain atleast one feasible solution of P0 as fast as possible. In such a case, we therefore employ adi erent strategy.
When we do not have any feasible solution of P0, we select the index
k~ = arg kmaxa(22)
2H F k13

as the next branching variable with the hope to obtain a feasible solution. In this case, we
just decompose the current partial problem P by xing yk~ to 0 and 1, and generate two newpartial problems. We do not compute a lower bound because an upper bound is unknown.
If jH1 j = p, this partial problem P becomes an arc capacitated network ow problem andwe solve it exactly.
Once we nd a feasible solution of P0 , we select the index

k~ = arg max
k62K

X

2

xk

(23)

as the next branching variable. Note that when we solve a Lagrangian relaxation problem
RR [; ] of P , we obtain as a byproduct an exact solution yk as shown in (17). If somek 62 K violate the capacity constraints (11) of P , selecting yk that corresponds to one ofsuch hubs k as a branching variable may yield a good feasible solution to P . By (23), we
actually select yk such that hub k 62 K violates the capacity constraints (11) of P0 mostseriously.
Now we summarize the branch-and-bound method for solving P0, which uses the depthrst search rule. We let A denote the set of active partial problems, i.e., partial problems
that have been generated but not tested yet.
[Branch-and-Bound Method ]Step 0: (Initialize) Find an optimal solution ( x ; y ) of the arc capacitated network ow

problem, in which candidates 1 ; 2;    ; p are selected as hubs. Let the optimal valueof this problem be the incumbent value z . If this problem is infeasible, then nd an
optimal solution (x ; y ) of another arc capacitated network ow problem, in which
candidates with p largest ak 's are selected as hubs. Let the optimal value of thisproblem be the incumbent value z  . If this problem is infeasible, set z  = 1. Set
 := 0; HF := H; H1 := ;; H0 := ; and A := fP0g. Select an integer parameter S ，and real parameters 0 <  < 1;  > 0 and  > 0.

Step 1: (Search) If

A = ;, then go to Step 7. Otherwise, select an active partial problem

P +1 2 A according to the depth- rst search strategy. Set  :=  + 1. If H1 = p orH0 = n ; p, then go to Step 2. If a feasible solution of P0 is unknown, then go to Step
5. If H1 < p, then go to Step 3.

Step 2: (Compute upper bound) Find an optimal solution ( x ; y ) of P with the objective

value z. If z < z  , then let z  := z; x := x ; y := y and go to Step 5.

14

Step 3: (Subgradient method)

(0) := 0, (0) := 0, and l := 0. Let G^ be upper bound of P . If an upperbound is not known, set G^ := z .
(b) Set Gmax
 := ;1.
(c) Solve RP [(l) ; (l) ] to obtain an optimal solution (x; y), and optimal value G . Ifmax
 max := (l) and max := (l) . If G > z  ,G > Gmax
 , then update G := G , (a) Setthen go to Step 6.
(d) If mod(l + 1; S ) = 0, then go to Step 3(e). Otherwise, compute Lagrange multipliers (l+1) and (l+1) according to (18){(20). Set l := l + 1 and go to Step3(c).
(l;S +1) ; (l;S +1) )Gmax + (1 ;  )G^ ; G ((e) Update G^ as G^ :=. If G^ < z 
2 ; 2
(l+1) := max ,or j(G^ ; Gmax
 )=G^ j < , then go to Step 5. Otherwise, set (l+1) := max ; l := l + 1 and go to Step 3(b).

P0 is unknown, select the next branching variable k~ according to (22). Otherwise, select the next branching variable k~ according to

Step 5: (Branching) If a feasible solution of(23).
Set A := A [ f(H1 [ fk~g; H0 ; HF nfk~g); (H1; H0 [ fk~g; HF nfk~g)g and go to Step 1.
Step 6: (Backtrack) Set
Step 7: (Terminate)

A := AnfP g and go to Step 1.

x and y are optimal solution of P0, and z  is the optimal value. If

z  = 1, then P0 is infeasible.

4 Numerical Experiments
In this section, we report some computational results. All programs were coded in MATLAB R12.1 (version6.1) with optimization toolbox [6]. They were run on a DELL DIMENSION 4400 computer with Pentium4 2.0AGHz processor operated under Windows XP
professional with 512 Mb DDR-SDRAM memory. We used the function linprog includedin the optimization toolbox, which solves a linear programming problem using the interior
point method.
We prepared the demand data based on the well-known U.S. 25 cities data evaluatedin 1970 by CAB (Civil Aeronautics Board). The data set, however, does not contain nodecapacities and xed costs, so we generated them based on the data reported by FAA (Federal
Aviation Administration) [12]. They report the passenger enplanements in the scal year of
15

hub#
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

Table 1: Scaled xed cost f~k and capacity ~ak
City(a)FY2000 (b)FY2012 (b);(a) Fixed Cost Capacity
Enplanements Forecast
f~k
~ak
3
3
3
10
10
10
Atlanta
39,37559,353 19,978
11.69
8.98
Baltimore
10,61816,497
5,879
3.44
2.50
Boston13,816
18,695
4,879
2.86
2.83
Chicago34,153
46,178 12,025
7.04
6.99
Cincinnati
9,186
18,749
9,563
5.60
2.84
Cleveland
6,746
10,935
4,189
2.45
1.65
Dallas-FW
28,066
41,759 13,693
8.01
6.32
Denver
18,88427,541
8,657
5.07
4.17
Detroit
17,87428,287 10,413
6.09
4.28
Houston
16,564
26,510
9,944
5.82
4.01
Kansas City
5,889
7,978
2,089
1.22
1.21
Los Angeles
32,332
47,294 14,962
8.76
7.16
Memphis
6,234
8,488
2,254
1.32
1.28
Miami
16,71624,156
7,440
4.35
3.65
Minneapolis
17,20325,914
8,711
5.10
3.92
New Orleans
4,900
6,609
1,709
1.00
1.00
New York
16,081
22,991
6,910
4.04
3.48
Philadelphia
13,023
20,748
7,725
4.52
3.14
Phoenix
18,652
31,252 12,600
7.37
4.73
Pittsburgh
10,521
14,319
3,798
2.22
2.17
St. Louis
1,4553
21,631
7,078
4.14
3.27
San Francisco
18,49925,634
7,135
4.18
3.88
Seattle14,22520,854
6,629
3.88
3.16
Tampa
8,200
11,384
3,184
1.86
1.72
WashingtonDC
8,502
15,146
6,644
3.89
2.29

2000 and the passenger forecast in the scal year of 2012. In Table 1, these data are shownin the column labeled \(a) FY2000 Enplanements" and the column labeled \(b) FY2012
Forecast", respectively, and the growth forecast of passenger enplanements from 2000 to2012 is shown in the column labeled \(b);(a)". Moreover, f~k denotes a scaled xed cost
obtained by setting the xed cost of New Orleans to be 1.0. Similarly, a~k denotes a scaled
capacity obtained by setting the capacity of New Orleans to be 1.0.
Since the data of node capacities are not available in the published documents, we determined the node capacity ak for each candidate byP d~2  a~ ; 8k 2 H;
(24)
ak = p  mink
k2H akwhere d~ is a scaled trip demand between OD pair , that is d~ = d = min2 d . In a similar
16

way, we determined arc capacities as follows. First we calculate the scaled arc capacities
by assigning node capacities ak equally to all arcs connected to it. In the computational
experiments, we assume that H = N , i.e., all nodes are candidate hubs. So the number ofarcs connected to each hub candidate is 2(n ; 1). Consequently, ~bij are calculated as~bij = ai + aj ;(i; j ) 2 A;
2(n ; 1)from which we determined arc capacities bij as

bij = ~bij ;(i; j ) 2 A;

where   1 is a parameter. It is clear that a problem with small  has tight arc capacity
constraints. Similarly, we determined fk , the xed cost of hub candidate k, as

fk = f~k ;k 2 H;

where  > 0 is a parameter. To see how the arc capacity constraints and xed costs a ect
the total optimal cost and hub locations, we solved the problem with di erent values of and .
The computational results are given in Tables 2 and 3. Table 2 shows the results with
 = 1000, where the ratio of xed cost to the total cost is slightly less than 5%. Table 3 shows
the results with  = 10000, where the ratio of xed cost to the total cost is approximately
10%. The column labeled \optimal hubs" and \optimal cost" show the set of optimal hubs
and the optimal total cost, respectively. The column labeled \# problems solved" shows thenumber of partial problems P solved exactly, i.e., the partial problems with jH j = p. The
column labeled \B&B CPU time" shows the total CPU time required to solve the problem
by the branch-and-bound method. The column labeled \% cost reduced" shows how the
optimal cost reduces as  is increased from  = 1.
A problem with small  often contains many infeasible partial problems because of its tightarc capacity constraints. In such a case, it is generally dicult to obtain a good upper boundin the early stage of branch-and-bound iterations, and hence it takes much computation
time. In fact, Tables 2 and 3 indicate that the computational time and the number ofpartial problems solved increase as  decreases. The optimal cost decreases monotonicallyas  increases, because arc capacity constraints become loose as  increases. However, when
 = 1000, the decrease of the optimal cost is less signi cant. We can observe that thedi erence of the optimal cost between problems with  = 1 and  = 3 is less than 2%
regardless of the value of p. In particular, when p = 4, the di erence is much smaller thanthat in the cases of p = 2 and p = 3. The reason for this phenomenon may be resorted to
the fact that p is relatively large and the ratio of xed cost to total cost is relatively small.
17

Regardless of the value of , Kansas City(#11) is always selected as a hub when p = 2,
and St. Louis(#21) and San Francisco(#22) are always selected when p = 3. Anotherselected hub changes in the order of New York(#17), Philadelphia(#18), Baltimore(#2) as increases. Note that these three cities are located closely to each other. These results
indicate that the optimal hub location does not change substantially when  = 1000.
On the other hand, when  = 10000, hub locations are signi cantly a ected by the valueof . Table 3 shows that the maximum di erence of the optimal cost between the two cases
 = 1 and  = 4 is nearly 5%. Although the di erence of 5% seems small, it brings aconsiderable amount of cost reduction when we deal with a large network. Figures 2 〜 4
display the results for problems with n = 25 and p = 2; 3; 4. The thickness of each arc
represents ow congestion on that arc. These gures show that the optimal hub location
varies with the value of . In particular, since the ratio of xed costs to the total costis relatively large compared with the case of  = 1000, a candidate with lower xed costis likely to be selected. As  increases, the hub location changes from New York(#17) toPittsburgh(#20), from Memphis(#13) to New Orleans(#16) and from New York(#17) toBoston(#3). These changes occur in such a way that a candidate is replaced by anothercandidate with lower xed cost. However, when p = 4, the selected hub location is drasticallychanged from Boston(#3) to San Francisco(#22) in spite of the expensive xed cost of San
Francisco. Note that, when p  3, no candidate located in the western area is selected. It
becomes possible to select such a candidate when  = 4, because of capacity relaxation. By
selecting a candidate located in the western area, we can reduce the cost of travel which
originates from and destines for the western area. This travel cost reduction is large enoughto compensate for the increase of xed cost; thereby it induces the total cost reduction.
Figure 4 shows that, when p  3, the trips originating from Seattle, San Francisco and Los
Angeles are long and the links emanating from these cities are congested.
From the above observation, we see that the total cost is a ected by arc capacities,
especially for problems in which the ratio of the xed cost to the total cost is large, i.e., is large. Taking into account the fact that the arc capacities are likely to be changed
frequently by various surroundings, it is important to select a robust or stable candidate.
To select such hubs, airlines may expect to have a constant revenue. When  = 10000,
Kansas City is always selected regardless of the values of p and . Its good location in thecentral area and the cheap xed cost may be the main reason why it is selected constantly.
On the other hand, when  = 1000, there is no such stable hub candidate which is chosen
independently of the value of p and . In this case, almost all candidates except for somevery expensive ones are likely to become hubs, because the ratio of the xed cost to the total
cost is small. When p = 2, Kansas City(#11) is always selected regardless of the value of .
18

Similarly, when p = 3, St. Louis(#21) is always selected. Moreover, when p = 4, selectedhubs remain the same completely. We may regard such hubs as stable for each p.
50

50

23

23
15
9

4

Latitude

40

8

11 21

22

13

35

12

19

7

30

10

5

6

3
40

1

-110

-100

-90

8

11 21

22

13

35

12

19

7

30

24
14

-120

-80

9

4

17
20 18
2
25

16

25

20
-130

15

45

Latitude

45

10

20
-130

3
17
20 18
2
25

1

16
24
14

25

-70

5

6

-120

-110

-100

-90

Longitude

Longitude(a)  = 1, hubs = (11,17)(b)  = 2, hubs = (11,20)

-80

-70

Figure 2: Results for n = 25; p = 2;  = 10000

5 Conclusion
In this paper, we have considered a one-stop capacitated hub location problem in which
both hubs and arcs have capacity constraints. We formulated the problem as a mixed 0{1integer programming problem and solve the problem using branch-and-bound method with
Lagrangian relaxation bounding strategy. We also made some computational experiments
using actual aviation data based on the well-known CAB data set and the future enplanements forecast reported by FAA. From our computational results, we see that the total cost
is a ected by arc capacities, especially for problems in which the ratio of the xed cost to thetotal cost is large. We also observe that the maximum di erence in the optimal cost between
the case of loose capacity constraints and that of tight capacity constraints is nearly 5 %.
These results indicate that it is important to incorporate arc capacities into the model when
we deal with a large network. Taking into account the fact that the arc capacities are likely
to be changed frequently by various surroundings, it is also important to select hubs thatare insensitive to those changes. In our computational experiments, some hub candidates
are always chosen independently of the number of selected hubs and arc capacities. It is aninteresting future research to develop an e ective method of selecting such a stable hub set
by using sensitivity analysis.
19

n p 
25 2 12
3
25 3 12
3
25 4 12
3n p 
25 2 12
3
4
25 3 12
3
4
25 4 1
2
3
4

Table 2: Computational Results for  = 1000
Optimal Optimal #Problems B & B % Cost
hubs
cost(105 )solvedCPU time reduced
11,17
3.3289171102.31{11,18
3.3044
174
62.59
0.74
2,11
3.2746
172
61.28
1.63
17,21,22
2.9741890
848.16{
17,21,22
2.9541904
573.500.67
18,21,22
2.9448890
574.880.98
6,12,13,17 2.8762
3287
4433.56{
6,12,13,17 2.8571
3463
3400.27
0.66
6,12,13,17 2.8566
3272
3248.89
0.68

Table 3: Computational Results for  = 10000
OptimalOptimal #Problems B & B % Cost
hubs
cost(105 )solvedCPU time reduced
11,17
3.8023174
108.64{11,20
3.6991
166
59.86
2.72
11,20
3.6374
160
57.88
4.34
11,20
3.6373
157
57.45
4.34
11,13,17
3.82371047
1041.38{
11,16,17
3.7882867
550.580.93
11,16,20
3.6916796
501.55
3.45
11,16,20
3.6535772
489.38
4.45
11,16,17,20 3.8790
2272
3630.61{
3,11,16,20
3.85692027
2148.84
0.57
3,11,16,20
3.82861630
1714.691.30
11,16,20,22 3.7695
1593
1648.42
2.82

20

50

50

23

23
15
9

4
8

11 21

22

13

35

12

19

7

30

10

5

-120

-110

-100

40

16

-90

9

4

17
20 18
2
25

1

25

20
-130

6

3
8

11 21

22

13

35

12

19

7

30

24
14

10

5

24
14

20
-130

-70

-120

-110

-100

-90

Longitude(a)  = 1, hubs = (11,13,17)(b)  = 2, hubs = (11,16,17)

50

23
15

45

9

4
40

8

11 21

22

13

35

12

19

7

30

10

5

6

3
17
20 18
2
25

1

16
24
14

25

20
-130

-120

-110

-100

-90

-80

-70

Longitude(c)  = 3, hubs = (11,16,20)

Figure 3: Results for n = 25; p = 3;  = 10000

21

3
17
20 18
2
25

16

25

-80

6

1

Longitude

LatitudeLatitude

40

15

45

Latitude

45

-80

-70

50

50

23

23
15
9

4

Latitude

40

8

11 21

22

13

35

12

19

7

30

5

40

-110

-100

13

35

12

-90

-80

20
-130

-70

19

7

6

5

3
17
20 18
2
25

1

16

10

24
14

-120

-110

-100

-90

-80

Longitude

Longitude(a)  = 1, hubs =(11,16,17,20)(b)  = 2, hubs=(3,11,16,20)

-70

50

23

23
15
9

4
40

8

11 21

22

13

35

12

19

7

30

10

5

6

3
40

1

16

25

-110

-100

8

11 21

22

13

35

12

19

7

30

-90

-80

9

4

17
20 18
2
25

24
14

-120

15

45

Latitude45

Latitude

11 21

22

25

50

20
-130

8

30

24
14

-120

9

4

17
20 18
2
25

1

25

20
-130

6

3

16

10

15

45

Latitude

45

10

20
-130

5
1

24
14

-120

-110

-100

-90

-80

Longitude

Longitude(c)  = 3, hubs =(3,11,16,20)(d)  = 4, hubs =(11,16,20,22)

Figure 4: Results for n = 25; p = 4;  = 10000

22

3
17
20 18
2
25

16

25

-70

6

-70

References[1] T. Aykin: \Lagrangian relaxation based approaches to capacitated hub-and-spoke network design problem", European Journal of Operational Research, 79, 1994, pp. 501{523.
[2] D. Bryan: \Extensions to the hub location problem", Geographical Analysis, 30, 1998,
pp. 315{330.
[3] J. F. Campbell: \A survey of network hub location", Studies in Locational Analysis Issue,
6, 1994, pp. 31{49.
[4] J. F. Campbell: \Integer programming formulations of discrete hub location problems",
European Journal of Operational Research, 72, 1994, pp. 387{405.
[5] J. F. Campbell: \Hub location and the p-hub median problem", Operations Research,
44, 1996, pp. 923{935.
[6] T. F. Coleman, M. A. Branch and A. Grace: Optimization Toolbox User's Guide, 1999,
The Math Works, Inc., U.S.A.
[7] A. T. Ernst and M. Krishnamoorthy: \Ecient algorithms for the uncapacitated single
allocation p-hub median problem", Location Science, 4, 1996, pp. 139{154.
[8] A. T. Ernst and M. Krishnamoorthy: \Exact and heuristic algorithms for the uncapacitated multiple allocation p-hub median problem", European Journal of Operational
Research, 104, 1998, pp. 100{112.
[9] A. T. Ernst and M. Krishnamoorthy: \An exact solution approach based on shortestpaths for p-hub median problems", INFORMS Journal on Computing, 10, 1998, pp.
149{162.
[10] A. T. Ernst and M. Krishnamoorthy: \Solution algorithms for the capacitated single
allocation hub location problem", Annals of Operations Research, 86, 1999, pp. 141{159.
[11] J. Ebery, M. Krishnamoorthy, A. Ernst and N. Boland: \The capacitated multiple
allocation hub location problem: Formulations and algorithms", European Journal ofOperational Research, 120, 2000, pp. 614{631.
[12] FAA: \2001 Aviation Capacity Enhancement Plan", Federal Aviation Administration,
Washington DC, 2001.

23

[13] J. W. Herrmann, G. Ioannou, I. Minis and J. M. Proth: \A dual ascent approach to
the xed-charge capacitated network design problem", European Journal of Operational
Research, 95, 1996, pp. 476{490.
[14] K. Holmberg and D. Yuan: \Lagrangian heuristic based branch-and-bound approach forthe capacitated network design problem", Operations Research, 48, 2000, pp. 461{481.
[15] P. Jaillet, G. Song and G. Yu: \Airline network design and hub location problems",
Location Science, 4, 1996, pp. 195{212.
[16] M. Minoux: Mathematical Programming: Theory and Algorithms, John Wiley & Sons,
Chichester, 1986.
[17] M. E. O'Kelly: \A quadratic integer program for the location of interacting hub facilities", European Journal of Operational Research, 32, 1987, pp. 393{404.
[18] M. E. O'Kelly, D. Bryan, D. Skorin-Kapov and J. Skorin-Kapov: \Hub network design
with single and Multiple Allocation: A computational study", Location Science, 4, 1996,
pp. 125{138.
[19] M. E. O'Kelly and D. Bryan: \Hub location with ow economies of scale", Transportation Research, 32B, 1998, pp. 605{616.
[20] B. T. Polyak: \Minimization of unsmooth functionals", USSR Computational Mathematics and Mathematical Physics, 9, 1969, pp. 14{29.
[21] M. Sasaki, A. Suzuki and Z. Drezner: \On the selection of relay points in a logistics
system", Asia-Paci c Journal of Operational Research, 14, 1997, pp. 39{54.
[22] M. Sasaki, A. Suzuki and Z. Drezner: \On the selection of hub airports on the airline
hub-and-spoke system", Computers and Operations Research, 26, 1999, pp. 1411{1422.
[23] D. Skorin-Kapov, J. Skorin-Kapov and M. O'Kelly: \Tight linear programming relaxations of uncapacitated P -hub median problems", European Journal of Operational
Research, 94, 1996, pp. 582{593.

24