Journal of the Operations ResearchSociety of Japan
2007, Vol. 50, No. 4, 428-443

MINIMUM UNDERFLOW ROUTING FOR STREAM-TYPE
COMMUNICATION SERVICES

Masaki Aida

Chisa Takano
Shunichi Osawa
Tokyo Metropolitan University(Received October 31, 2006; Revised July 12, 2007)
Abstract
The Internet of today supports various types of communication services including not only
conventional data communication services but also stream-type communication services. Typical examples
of stream-type services include interactive voice such as Voice over IP (VoIP) and live video delivery. Stream
services require real-time packet transmission, and the continuous playback of packets at the receiver side.
As a result, they are not only sensitive to the absolute value of the delay, but also sensitive to delay
variations. This paper addresses the problem of optimal routing for stream-type communication services.
Optimality is discussed in terms of the continuous playback of packets. From network observations, we
assume that the end-to-end delay statistics conform to a normal distribution. We model a network as aweighted graph with its link weights representing link delays. In the course of the analysis, we show that
this type of routing optimization problem can be formulated as a process of searching for a speciﬁc point
in a coordinate system deﬁned by the mean and variance of the end-to-end delay. This paper presents aneﬃcient algorithm for ﬁnding the optimal point in this coordinate system.

Keywords: Telecommunication, routing, underﬂow, stream traﬃc, VoIP

1. Introduction
The Internet of today has been regarded as a multimedia network that supports various typesof communication services. In particular, stream-type communication services, including
Voice over IP (VoIP) and live video delivery, are very diﬀerent from conventional data
communication services. Stream-type services require real-time packet transmission, andthe continuous playback of packets at the receiver side.
In the past, various ways have been oﬀered for routing in the Internet or, in general,
packet switched networks. The main goal is to minimize the summation of a certain quantity
associated with each link along a path, for example, the number of hops or the mean value of
end-to-end delay across the network [6]. This class of problem is the shortest path problem
and is quite solvable. It can be eﬃciently solved by applying the well-known Dijkstra’salgorithm [3].
Associated with each stream-type communication service, there exists the maximum
permissible delay (or allowable delay) that should be guaranteed to obtain satisfactory
quality of service (QoS) [1, 4]. Let us take interactive voice communication as an example[see Figure 1]. Voice packets received by the listener side terminal are stored temporarily
in a buﬀer before they are played back, to smooth out the packet delay jitter. The buﬀer
is called the playout buﬀer. Playback is performed under the control of a continuous clock
source. The required amount of storage must be large enough to make the probability of
buﬀer overﬂow/underﬂow reasonably small — a situation that may occur with excessive
packet delay variation. Furthermore, packets that experience delays longer than a certain
428

429

Minimum Underflow Routing
traffic with delay variations

stream traffic

continuous playback

network
senderplayout buffer

receivermaximum permissible delay

Figure 1: System model
value will be too late for play-back and must be discarded at the receiving buﬀer.
Thus, to obtain a satisfactory service quality, the routing should be designed to take the
distribution of delay into account. In this paper, we investigate end-to-end delay distribution
by means of observations of an actual network (Sec. 2). Based on these observations, weformulate the routing problem as the minimization of the probability that packets exceed
the maximum permissible delay (Sec. 3). We model a network as a graph with its linkweights representing delays and assume that the end-to-end delay statistics follow a normal
distribution. Our problem is regarded as a stochastic version of the shortest path problem
on a weighted graph. This class of problem requires a diﬀerent formulation to that which
may be used for optimal routing based simply on minimization of mean delay. As will be
described later, Dijkstra’s algorithm cannot be applied directly. Ichimori et.al.[5] consider
a related problem which can be dealt with using a minimum spanning tree. We will showthat the treatment proposed in their paper can be eﬃciently used as one tool to solve ourproblem (Sec. 4). However, their treatment is not suﬃcient to solve the stochastic shortest
path problem by itself. A shortest path problem is a more diﬃcult class of problem thanthe minimum spanning tree problem and requires some extensions. Next, we show that the
optimization of routing in our problem is equivalent to searching for a speciﬁc point in acoordinate system determined by the mean and variance of the end-to-end delay (Sec. 5).
Using this formulation, we construct a simple and eﬃcient algorithm for solving the problem(Sec. 6) and show that the residual estimation error is bounded (Sec. 7). Finally, we describe
the simulation results (Sec. 8).
2. Measurements of End-to-End One-Way Delay
2.1. Measured data
We conducted observations of the end-to-end delay of packets by using a real network. The
path from the PC of an asymmetric digital subscriber line (ADSL) customer to another PC
in a company’s LAN was used in the experiments, where both PCs were GPS-synchronized.
The LAN for the ADSL customer and the company’s LAN were connected via two ISPs
and the path between the LANs consisted of 15 hops. 64-byte UDP packets were sent fromthe ADSL customer to the PC in the company’s LAN according to a Poisson process andtheir one-way end-to-end delays were monitored for 900 seconds for each experiment period.
The experiment was conducted on a weekday in December 2003.
Five typical results are listed in Table 1. The load of the network was light in the ﬁrst
two cases and heavy in the last two cases. Measurements #1 and #2 were conducted atnight, and #4 and #5 were in busy hours. Measurement #3 was the result when there wasa medium load on the network.
c Operations Research Society of Japan JORSJ (2007) 50-4
430

M. Aida, C. Takano & S. Osawa

Table 1: Results of measurement of end-to-end one-way delay (ms)reference measurement
mean
min.
max. usage ofID
perioddelay delay
delay network#1
00:00–00:15
5.319 4.942 11.948
light#2
02:30–02:45
5.279 4.944
9.948
light
#3
10:00–10:15
42.258 4.967 138.136 medium#4
10:45–11:00 130.590 28.516 187.936 heavy#5
11:45–12:00 116.767 27.746 191.342 heavy
2.0

2.0

1.6

1.6

1.2

1.2

0.8

0.8

0.4

0.4

0

0

1

2

3

4

5

6

7

8

9 [ms]

0

-2

-1.5

-1

-0.5

0

0.5

1

1.5

Figure 2: Distributions of end-to-end one-way delay for measurement #1 (left: distribution of measured delay on a linear scale, right: distribution for the logarithm of variable
component of delay)

2.2. Observed distribution of end-to-end one-way delay
First, we investigate the cases where the network load is light. The left-hand graph of Figure
2 shows a distribution of the measured delay for measurement #1. We subtracted 4.87 ms
as the constant component of the delay and show the histogram for the logarithm of residualcomponent of the delay in the right-hand graph of Figure 2. The right-hand graph implies
end-to-end delay consisting of a constant component and a variable component that follows
a lognormal distribution. Figure 3 shows a comparison between the measured, empirical
distributions of the variable component of end-to-end one-way delay for measurements #1
and #2, and lognormal distributions. These results imply that the variable component of
end-to-end delay follows a lognormal distribution when the network load is light.
Next, we examine the cases where the network load is heavy. Figure 4 shows a comparison
between the empirical distributions of end-to-end one-way delay for measurements #4 and#5, and normal distributions. These results imply that the end-to-end delay follows anormal distribution when the network load is heavy.
Finally, we consider the case where the load of the network is medium. Figure 5 shows
a comparison between the empirical distribution of the variable component of end-to-end
one-way delay for measurement #3, and theoretical distributions, normal and lognormal.
From this result, we cannot determine the distribution of the end-to-end delay when there
is a medium load on the network.
Table 2 summarizes the above results. Since controlling the QoS of streaming services isimportant when the network is congested, we should consider the situation where the loadof the network is high. So, we assume the end-to-end delay follows a normal distribution.
c Operations Research Society of Japan JORSJ (2007) 50-4


431

Minimum Underflow Routing
1

1
Lognormal distribution

Lognormal distribution

0.8

0.8
Empirical distrubution

0.6

Empirical distrubution

0.6

of the variable componentof end-to-end delay

0.4

of the variable componentof end-to-end delay

0.4

0.2

0.2

0
-3
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
Logarithm of the variable component of end-to-end delay

0
-3
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
Logarithm of the variable component of end-to-end delay

Figure 3: Comparisons between empirical distributions of the variable component of end-toend one-way delay and lognormal distributions (left: measurement #1, right: measurement
#2)
1

1

Empirical distrubution

0.8

Empirical distrubution

0.8

of end-to-end delayof end-to-end delay

Normal0.6

distribution

Normal
0.6

0.4

0.4

0.2

0.2

0

0

20

40

60

80
100 120 140
End-to-end delay

160

180

0

0

distribution

20

40

60

80
100 120 140
End-to-end delay

160

180

Figure 4: Comparisons between empirical distributions of end-to-end one-way delay andnormal distributions (left: measurement #4, right: measurement #5)
Even if the assumption is not appropriate in all cases, we can assume that a light load on
the network does not have a crucial eﬀect on the QoS.
3. Optimal Routing Problem
We represent a network as a graph G = G(V, E) where V is the set of nodes and E is the setof links connecting the nodes. Each link denoted by i has a weight d(i) that represents thedelay experienced by packets on the link. In general, the value of d(i) ﬂuctuates dynamically.
We assume that d(i) is an independent random variable characterized by its mean m(i) andvariance v(i). We further assume that the mean m(i) and variance v(i) are stable over time.
We refer to d(i) as a stochastic weight. We treat a shortest path problem between given
two nodes in the following part of this paper.

Table 2: Distributions of end-to-end one-way delayload of the network
distribution
light
constant + lognormal
medium?
heavy
normalc Operations Research Society of Japan JORSJ (2007) 50-4
432

M. Aida, C. Takano & S. Osawa1
0.8
Empirical distrubution
0.6

of the variable componentof end-to-end delay

0.4
Normal distribution
0.2
Lognormal
distribution
0
-2
-1
0
1
2
3
4
5
Logarithm of the variable component of end-to-end delay

Figure 5: Comparisons between empirical distributions of end-to-end one-way delay andnormal distributions (for measurement #3)
Before considering the shortest path problem between two given nodes, we brieﬂy add
supplemental comments about the assumptions of independence and stationarity of d(i).
We do not assume stability over a long time period; we are interested only in the limited
period of the busy hour(s). Also in traditional telephone networks, stochastic properties oftraﬃc are assumed to be stable over a limited period even if they change dramatically overa longer time period. We also assume independence of ﬂows in the busy hour(s). In thisperiod, when there is a high volume of traﬃc, a single ﬂow has a relatively small impact on
the total amount of traﬃc. Even if the assumption of independence is not true for link-bylink delay statistics, then we can reconsider the link delay as the delay of a segment of a
network that is of a relatively large scale, e.g. the delay of each Internet Service Provider(ISP) network. This is because it is not easy to see how there could be short time dependenceof stochastic properties between diﬀerent large scale subnetworks.
Let π denote a path between two given nodes in G and Π = {π} be the set of all paths
between the two given nodes. Given any path π, based on the independence assumption,
the path weight D(π), its mean M (π), variance V (π), and standard deviation S(π) can be
expressed as:D(π) =
d(i), M (π) =
m(i),
i∈π

V (π) =



i∈π

v(i),

and S(π) =V (π).

(3.1)

i∈πHere, we assume the end-to-end delay D(π) follows a normal distribution
D(π) ∼ N (M (π), V (π)).

(3.2)

Conventionally, routing problems have been formulated in terms of ﬁnding a particularpath π between two nodes from the set of paths Π that minimizes the mean end-to-end
delay. In this formulation, the problem is expressed in the form:
min M (π).
π∈Π

(3.3)

This problem depends only upon the mean value of delay because (3.3) does not takethe variance of delay into consideration. Previous studies show that this type of problem
of ﬁnding the shortest path can be solved eﬃciently by applying dynamic programming
techniques; e.g. Dijkstra’s algorithm.
c Operations Research Society of Japan JORSJ (2007) 50-4


433

Minimum Underflow Routing

In contrast, our aim is to take the variance of delay into consideration, so we formulate
our minimum underﬂow routing problem as follows:
min P{D(π) > D0 },

(3.4)

π∈Π

where D0 is a given positive real constant representing the maximum permissible delay. The
purpose of our optimization is to minimize the probability that the end-to-end delay exceeds
a predeﬁned value. Here, we assume that D0 satisﬁes
D0 ≥ min M (π),

(3.5)

π∈Π

by considering real applications to communication networks. By standardization of D(π),
we can obtain an equivalent formulation:
D0 − M (π)
.
π∈Π
S(π)

(3.6)

max

Hereafter, we call the function to be maximized the “Evaluation function” of the path π:E(π, D0 ) :=D0 − M (π)
.
S(π)

(3.7)

4. Application of Dynamic Programing
Unfortunately, the evaluation function E(π, D0 ) still has properties that make it hard totreat. This section explains what this diﬃculty is and how it can be overcome.
Imagine a situation where there are three nodes (A, B, and C) which are connected bythree links as shown in Figure 6. Each link and each path have weights that are evaluatedby E(π, D0 ). We further assume that link c is the optimal path between nodes B and C.
It can easily be shown that, if the evaluation function E(π, D0 ) is used to evaluate thepath/link weights, the optimal path between node A and node C does not always include
link c which is the optimal path between node B and node C. We show an example using
Figure 6. We assume that link a has (m(a), v(a)) = (1, 1). Similarly, links b and c have(m(b), v(b)) = (7, 4) and (m(c), v(c)) = (8, 1), respectively. In addition, we set D0 = 10.
Then, the values of the evaluation function of the paths between nodes B and C areE(b, 10) =10 − 7√= 1.5,
4and E(c, 10) =10 − 8√= 2,
1respectively. Since E(b, 10) < E(c, 10), link c is better. On the other hand, values of theevaluation function of the paths between nodes A and C are
E(a → b, 10) =10 − (1 + 7)√
 0.894,
1+4

and E(a → c, 10) =10 − (1 + 8)√
= 0.707,
1+1

respectively. Since E(a → b, 10) > E(a → c, 10), the best path between nodes A and C
includes link b. It means that the principle of optimality [2], which is the basis for applying
dynamic programming techniques, is not available when this evaluation function is used.
Hence, the well-known, Dijkstra’s algorithm cannot be applied directly to our problem
because the algorithm employs the optimality principle. To avoid the above diﬃculty, we
introduce a new measure for the link and path weights. The following approach follows
c Operations Research Society of Japan JORSJ (2007) 50-4
434

M. Aida, C. Takano & S. Osawabacnode Anode Bnode C

Figure 6: Example of network
Ichimori et.al. [5] whose work deals with the stochastic minimum spanning tree problem.
This technique is also eﬀective in solving our problem. We introduce a positive parameter
x (which has the dimension of 1/(time)) and deﬁne the weight of the individual link i as:
wx (i) := m(i) + x · v(i).

(4.1)

Following this deﬁnition, the x-weight of a path π can be expressed as:
Wx (π) = M (π) + x · V (π).

(4.2)

With this formulation of link and path weights, an end-to-end path weight can be obtainedby simply summing up the weights of all links that comprise the path, i.e.,m(i) + xv(i) =
wx (i).
(4.3)
Wx (π) =
i∈π

i∈π

i∈π

Hereafter, we call wx (i) the x-weight of link i.
Next, we deﬁne the concept of the x-optimal path between a given pair of nodes. A path
πx-opt (x) is x-optimal if its x-weight satisﬁed the following minimization for a given value ofx ≥ 0.
min Wx (π).
π∈Π

(4.4)

In general, the conﬁguration of the path πx-opt varies with the value of x. Let Πx-opt denote a
set of such optimal paths and let πopt denote the real optimal path that satisﬁed the original
problem (3.6).
The x-weight thus deﬁned has a useful characteristic which allows our problem to besolved. We show this in the following theorem.
Theorem 4.1 The true optimal path πopt is included in the set Πx-opt of x-optimal paths.
πopt ∈ Πx-opt .

(4.5)

Proof
We prove this by contradiction. Suppose the real optimal path πopt is not included in theset Πx-opt . We then haveD0 − M (πopt ) = E(πopt , D0 ) S(πopt ),D0 − M (πx-opt (x)) < E(πopt , D0 ) S(πx-opt (x)),
from (3.7), which gives the following inequality.
M (πx-opt (x)) − M (πopt ) > E(πopt , D0 ) (S(πopt ) − S(πx-opt (x))) .
c Operations Research Society of Japan JORSJ (2007) 50-4


(4.6)

Minimum Underflow Routing

435

If we arbitrarily select the value of x as x = E(πopt , D0 )/(2S(πopt )), there should exist a
corresponding x-optimal path πx-opt (x ) in Πx-opt which should satisfy
Wx (πx-opt (x )) < Wx (πopt ).

(4.7)

However, it follows thatWx (πx-opt (x )) − Wx (πopt )= M (πx-opt (x )) + x · V (πx-opt (x )) − (M (πopt ) − x · V (πopt ))E(πopt , D0 )(S(πopt ) − S(πx-opt (x )))2 ≥ 0.
>2S(πopt )(4.8)

Thus, a contradiction is derived.
This theorem enables us to take the following steps to solve our problem: ﬁrst, we solve
for the set of x-optimal paths for the x-weight problem to which the optimality principle
holds and can be solved eﬃciently using Dijkstra’s algorithm; then we search for the true
optimal path from among the x-optimal paths.
5. Representation of Paths on M -V Plane
The main diﬃculty with the above derived approach is how to ﬁnd out the optimal pathfrom all the x-optimal paths for a given graph. To overcome this diﬃculty, we introduce, inthis section, a new coordinate system to analyze our problem. This oﬀers a useful means ofsystematically searching for the x-optimal paths and the true optimal path.
Speciﬁcally, our problem is formulated as searching for a speciﬁc point in the plane.
In the following, we show how all x-optimal paths are represented as points in the plane.
We then examine the geometrical properties of the points in the plane that can be used toshrink the domain where the target point should be searched for.
The x-weight deﬁned by (4.1) is a linear function of x. In general, any linear function
is completely determined by specifying two parameters, i.e., its gradient and the point ofintersection with the vertical axis. Here, we introduce a plane that is deﬁned by the M axis and V -axis, representing, respectively, the intersection and the gradient of the x-weightfunction, which in turn represent the mean delay of the corresponding path and its variance.
On this plane, referred to as the M -V plane, each path π is represented by a single point
P = (M (π), V (π)) as shown in Figure 7.
Next let us examine the contours of the original evaluation function, (3.7), on the M -V
plane. For a given value E of the evaluation function, the corresponding contour is expressed
as
V =

1
(D0 − M )2 .
2
E

(5.1)

The resulting shape of the contour is parabolic. Several examples of contours, for diﬀerent
values of E, are shown in Figure 7.
Following the above formulation, our purpose can be restated as ﬁnding the point withthe largest value of E. The point is expected to be found towards the lower left corner on
the M -V plane.
To systematize this searching process, we ﬁrst examine the basic relationship between
two points, P1 and P2 , on the M -V plane with the following two lemmas. In the following
discussion, a path πi is represented by the point Pi whose Cartesian coordinates are givenby (Mi , Vi ), where Mi = M (πi ) and Vi = V (πi ).
c Operations Research Society of Japan JORSJ (2007) 50-4


436

Variance

M. Aida, C. Takano & S. Osawa

O

D_0

Mean

Figure 7: Distribution of paths and contours for evaluation function
Lemma 5.1 Let π1 and π2 be diﬀerent paths belonging to Πx-opt . Let points P1 and P2 on
the M -V plane represent the paths π1 and π2 respectively. If there is some x = x0 > 0 thatsatisﬁesWx (π1 ) = Wx (π2 ) for x = x0 ,(5.2)

then the gradient of the line segment P1 P2 is negative and its value is −(1/x0 ).
ProofBy applying (4.2) to (5.2), x0 may be calculated asx0 = −

M1 − M2.
V1 − V2(5.3)

Therefore the gradient of the line segment P1 P2 is:gradient =1
V1 − V2= − < 0.
M1 − M2
x0

(5.4)


Lemma 5.2 Let π1 and π2 be diﬀerent paths belonging to Πx-opt , which are the x-optimal
paths for x = x1 and x = x2 , respectively. If x1 < x2 , there is some x = x0 that satisﬁes
Wx (π1 ) = Wx (π2 ),(x = x0 ),

(5.5)

and
x 1 < x0 < x2 .

(5.6)

Proof
From the assumption, we haveWx (π1 ) < Wx (π2 ),Wx (π1 ) > Wx (π2 ),(x = x1 ),(x = x2 ).

(5.7)

Since both Wx (π1 ) and Wx (π2 ) are linear functions of x, the two lines intersect, so there issome value of x such that Wx (π1 ) = Wx (π2 ), where x1 < x < x2 .
We next examine the geometrical distribution of points Pi , representing the collectionof x-optimal paths, over the M -V plane. This distribution has an outstanding property
described by the following theorem. Before entering into details, we should mention someproperties that can be understood intuitively:c Operations Research Society of Japan JORSJ (2007) 50-4
Minimum Underflow Routing

437

1. There exists only a ﬁnite number of x-optimal points.
2. Each point presents x-optimality for a certain range of x, i.e., for xi < x < xi+1 . So,
supposing that the complete set of x-optimal paths has been obtained, we can number
the points such that Pi is optimal for xi < x < xi+1 , where i = 1, 2 . . .
Theorem 5.1 Suppose that the complete set of the x-optimal points has been obtained andthey are numbered such that a point representing the x-optimality with a larger value of x is
given a greater number. Connect these Pi on the M -V plane by straight lines in ascendingorder and the resulting shape is convex and monotonically decreasing [see Figure 8].
Proof
Let Pi , Pj and Pk be diﬀerent points on the M -V plane, where they are determined tobe x-optimal for the value x of xi , xj and xk respectively, and 0 < xi < xj < xk . Let
gradient(Pi , Pj ) and gradient(Pj , Pk ) denote the gradients of line segments Pi Pj and Pj Pk ,respectively. Because these points correspond to paths in Πx-opt , we can always ﬁnd somereal values a, b for x such thatWx (πi ) = Wx (πj ),Wx (πj ) = Wx (πk ),(x = a),(x = b).

(5.8)

For these to be satisﬁed, from lemma 5.2, a and b must satisfy the relationship:
0 < x i < a < xj < b < xk .

(5.9)

From lemma 5.1, we have
1
gradient(Pi , Pj ) = − ,a1gradient(Pj , Pk ) = − ,b(5.10)

hence
gradient(Pi , Pj ) < gradient(Pj , Pk ) < 0.

(5.11)

That is, the gradient is less steep on the right side of the connected lines. This property
holds true for any combination of x-optimal points and hence the shape is convex.
Using the above results, we now develop an algorithm that locates the optimal point in
the M -V plane. Our algorithm starts with a known set of x-optimal points on the plane.
The main steps carried out in the algorithm consist of:
1. determining regions in which the candidate points are to be searched for,2. determining the region in which the search process is to be applied ﬁrst, and3. shrinking the search regions based on the newly obtained x-optimal points.
Assume that x-optimal points Pi , Pj , Pk and Pl are known on the M -V plane, where
xi < xj < xk < xl . We consider the problem of searching for a new point Pn for xj <xn < xk . The domain where Pn should be searched for is a triangle [see Figure 9] and its
location is determined in the following manner. Two vertices of its triangular domain arefound trivially to be Pj and Pk . We deﬁne the new vertex Qjk as the intersection point of
two extended lines, Pi Pj and Pk Pl . The reason for this is easily understood by considering
the fact that if Pn is to exist for xj < xn < xk , the connecting curve Pj Pn Pk must be convex
as dictated by theorem 5.1.
Within this triangle, the candidate point that has the potential to oﬀer the highest value
in terms of the evaluation function (3.7) is Qjk . If the evaluation with Qjk gives a larger
value than with the other two points, Qij and Qkl , the region is worth searching ﬁrst.
c Operations Research Society of Japan JORSJ (2007) 50-4


438

Variance

M. Aida, C. Takano & S. Osawa

O

D_0

Mean

Figure 8: Distribution of x-optimal paths
For a given set of known points, we can identify as many triangles as there are line
segments. Each newly found vertex must have a corresponding value of the Evaluation
function. The triangle with the new vertex which is evaluated to have the largest value isthe most promising place to search for a new x-optimal point, and this point is expected toyield the greatest value of the Evaluation function.
Once we ﬁnd a new x-optimal point, we can restart the algorithm with the increased
number of known points. In this case, following the same reasoning as the above, we canshrink the search domain as shown in Figure 10. Note that ﬁnding a single point results inthe shrinkage of not only the original triangle but also the two adjacent triangles. As thecycles of the algorithm are repeated, the search domain continues to shrink until the pointwith the highest evaluation, or the true optimal point, is determined.
For the given triangle (Pj , Pk , Qjk ), we can ﬁnd a new x-optimal point in the following
manner. We use
1x=−(5.12)
gradient(Pj , Pk )to solve for the x-optimal path by applying Dijkstra’s algorithm. The justiﬁcation for this isgiven in Appendix A. If we can ﬁnd a new path other that those corresponding to Pj or Pk ,this would be the one with the highest evaluation value. Otherwise, there are no x-optimal
points in this region, and this region should be excluded from further investigation.
6. Formal Description of Algorithm
In this section, based upon the previous discussion, we introduce a formal description of
our algorithm. As seen in Figures 9 and 10, each search region forms a triangle, two of thevertices or which are x-optimal paths and the other is the point that gives the upper-boundof the evaluation E(π, D0 ) in the triangle. In the following algorithm, Eopt takes the valueof the largest evaluation E(π, D0 ) among x-optimal paths that have already been searched
for, and Eub takes the value of the upper-bound of the evaluation E(π, D0 ) for all x-optimal
paths (the latter case means eﬀectively ﬁnding the shortest path using weight v(i) only).
The algorithm is as follows:
Step 1 Given two nodes on G, ﬁnd the initial x-optimal paths between them, ﬁrst with x
set to 0 and then with x set to a very large value (this means that ﬁnding the shortest
paths using weight v(i) only). The path ﬁnding processes can be done using Dijkstra’salgorithm. Calculate the corresponding evaluations E(π, D0 ) by applying (3.7) and set
c Operations Research Society of Japan JORSJ (2007) 50-4


439

Minimum Underflow Routing
P_i

Q_ij
P_j

P_n
P_k

Q_jk

P_l
Q_kl

Figure 9: Search domain I
P_i

Q_ij

P_j

P_n
Q_jn
P_k
Q_nk
P_l
Q_kl

Figure 10: Search domain II
Eopt to the larger value. In addition, set Eub to the largest evaluation for unsearchedpaths.
Step 2 Determine search domains from a set of known x-optimal points, and evaluateE(π, D0 ) for each new vertex Qjk . For the domain (Pj , Pk , Qjk ) with the highest evaluation, update Eub , and ﬁnd a new x-optimal path applying Dijkstra’s algorithm with x
set to −(1/gradient(Pj , Pk )). If there remains no candidate domain, go to step 4.
Step 3 If a new path exists, we plot the corresponding point Pn on the M -V plane. Calculate E(π, D0 ) for the new point Pn and if it is larger than Eopt , update Eopt to this value.
Exclude from the candidate list those domains whose evaluation values are smaller thanEopt . Otherwise no new path can be found, and we exclude the domain from the searchdomain list. In both cases, go back to step 2.
Step 4 Terminate the algorithm.
Figure 11 depicts the initial steps of the algorithm. The left-hand ﬁgure shows the M -V
plane for Step. 1 after two points have been calculated. If a third point is found in thetriangle triangle on the left, the new search region is restricted as shown in the right-handc Operations Research Society of Japan JORSJ (2007) 50-4


440

M. Aida, C. Takano & S. Osawa
P_1

Variance

Variance

P_1

Q_13

Q_12
O

Q_32

P_2
Mean

P_3

O

P_2
Mean

Figure 11: Initial steps in algorithm
ﬁgure.
7. Residual Error in the Algorithms
The main strength of our algorithm is that we can control the residual error caused bypremature termination of the algorithm. Let us consider the following two probabilities,
that end-to-end delay exceeds the maximum permissible delay;
 2
 +∞
1
x
popt = √
exp −
dx,
(7.1)
2
2π Eopt
 2
 +∞
1
x
exp −
dx,
(7.2)
pub = √
2
2π Eubwhere Eopt is current optimal value and Eub is the upper-bound of the optimum value
obtained from the algorithm. We deﬁne the error as the diﬀerence of the above probabilities,(popt − pub ). The residual error can be calculated as:(the maximum residual error) < popt − pub
 2
 Eub
1
x
=√
exp −
dx.
2
2π Eopt(7.3)

As the algorithm iterates, Eub is decreased and Eopt is increased (or remains unchanged).
Therefore the error decreases with each cycle in the algorithm and so we can stop thealgorithm when the error is suﬃciently small.
8. Simulation Results
This section describes the results of a simulation for investigating the eﬃciency of the minimum underﬂow routing algorithm and the number of iterations required for our algorithm
derived in Sec. 6. Figure 12 shows the network model used in the simulation. The modelhas 100 nodes and they are conﬁgured as a 10 × 10 lattice topology. We are searching forthe communication path is searched from the sender node 00 to the receiver node 99. The
minimum number of hops is 18, which is comparable in size to a real network environment.
Each link has two diﬀerent weights that correspond to the mean and variance of thelink delay. The values of the weights are set in advance as independent uniform random
numbers distributed in the range (0, 1).
We choose the maximum permissible delay as D0 = 12, 13, or 14. For each value ofD0 , we have conducted path searches 100 times and compared the paths obtained fromc Operations Research Society of Japan JORSJ (2007) 50-4


441

Minimum Underflow Routing
sender
00

01

02

08

09

10

11

12

18

19

20

21

22

28

29

80

81

82

88

89

90

91

92

98

99

receiverFigure 12: Simulation model
0.014
Dijkstra’s algorithm (optimization of (3.3))
0.012
Our algorithm (optimization of (3.4))
0.01

0.008

0.006

0.004

0.002

0
D_0 = 12

D_0 = 13

D_0 = 14

Figure 13: Simulation results: Probability that the end-to-end delay exceeds D0
optimizations (3.3) and (3.4). First, Figure 13 sshows a comparison of the probabilities thatthe end-to-end delay exceeds the maximum permissible delay D0 . The probability values
were obtained from the average of 100 path searches. From this ﬁgure it may be seen thatthe minimum underﬂow routing algorithm is most eﬀective for larger values of D0 .
In addition, Table 3 shows the number of iterations, which were required to search forthe optimal path using our algorithm. Here, the number of iterations means the number
of applications of Dijkstra’s algorithm before the optimal path is determined. The average
number of iterations is around 6, and the maximum is 9. From these results, we can expect
the algorithm to converge fairly rapidly in practical networks.

Table 3: Number of iterations
maximum
number of
number of iterations
permissible delay measurements average max. min.
D0 = 12
100
5.96
9
3
D0 = 13
100
6.28
9
4
D0 = 14
100
6.18
9
4

c Operations Research Society of Japan JORSJ (2007) 50-4
442

M. Aida, C. Takano & S. Osawa

9. Conclusion and Discussion
The QoS of stream-type communication services like interactive voice is sensitive to delay
variations. We have presented in this paper an eﬃcient algorithm for solving the stochastic
shortest path problem for which optimization function is given by (3.4). This is applicableto routing technologies for stream-type communication services to maintain their QoS. The
primary results of this paper may be summarized as follows. First, we investigated the endto-end delay distribution by measuring an actual network, and showed it follows a normal
distribution when the network is congested. Next, we introduced the notion of x-weight
to make the application of dynamic programming techniques possible, and we introduced
the M -V plane and re-state the problem in a graphical context, in which the problem isformulated as the process of searching for a speciﬁc point on the M -V plane. Based onthese developments, we proposed an eﬃcient algorithm to search for the point on the M -V
plane. In this process, we demonstrated that the remaining estimation error is bounded.
(Eub − Eopt ) decreases monotonically with the number of times the algorithm is applied andthe residual error (7.3) decreases similarly. Therefore, we can stop the algorithm taking into
account factors such as the network size, available processing power, etc., while controlling
the remaining estimation error caused by the untimely termination of the algorithm. Finally
we showed simulation results for this algorithm, which lead us to expect that the algorithm
will converge fairly rapidly in practical networks.
Acknowledgements
We would like to thank Mr. Teruyuki Kubo of NTT Software Corporation, Dr. Keisuke
Ishibashi and Mr. Ichizo Nakamura of NTT Corporation, for their enlightening discussions
and invaluable comments.
A. Determination of the Value of x
The x-weight function of each path P , Wx (P ), is a linear function of x and can be represented
on the x-Wx plane by a straight line. Suppose that two x-optimal paths, P1 and P2 havebeen found for x = x1 and x2 , respectively, and the two lines intersect at x = x3 . If anotherx-optimal path P3 is found for x1 < x < x2 , then the corresponding line should fall within
the shaded region in Figure 14. Thus, the corresponding x-weight function Wx (P3 ) shouldhave the following property.
Wx (P1 ) < Wx (P3 ) for x = x1 ,Wx (P2 ) < Wx (P3 ) for x = x2 ,Wx (P3 ) < Wx (P2 ) and Wx (P3 ) < Wx (P1 ) for x = x3 .

(A.1)
(A.2)
(A.3)

It can easily be understood that we can determine the existence of such a path by calculating
an x-optimal path with x set to x3 . The value of x3 is calculated asx3 = −

1.
gradient(P1 , P2 )(A.4)

References[1] M. Aida, I. Nakamura, and T. Kubo: Optimal routing in communication networks withdelay variations. IEEE Infocom ’92, (1992), 153–159.
[2] R. Bellman: Dynamic Programing (Princeton University Press, Princeton, NJ, 1957).
c Operations Research Society of Japan JORSJ (2007) 50-4


443

Minimum Underflow Routing
W_x
W_x(P_1)

W_x(P_2)

O

x_1

x_3

x_2

x

Figure 14: Path representation on (x, Wx )-plane[3] D. Bertsekas and R. Gallager: Data Networks (Prentice-Hall, Inc., Englewood Cliﬀs,
NJ, 1987), 322–323.
[4] G. Carle and E.W. Biersack: Survey of error recovery techniques for IP-based audiovisual multicast applications. IEEE Transactions on Networking, 11 (1997), 24–96.
[5] T. Ichimori, S. Shiode, H. Ishii, and T. Nishida: Minimal Spanning Tree with Normal
Variates as Weight. Journal of Operations Research Society of Japan, 24 (1981), 61–65.
[6] M. Schwartz and T. E. Stern: Routing Techniques Used in Computer Communication
Networks. IEEE Transactions on Communications, COM-28 (1980), 539–552.
Masaki Aida
Graduate School of System Design
Tokyo Metropolitan University
6-6 Asahigaoka, Hino-shi
Tokyo 191-0065, Japan
E-mail: maida@sd.tmu.ac.jp

c Operations Research Society of Japan JORSJ (2007) 50-4