IEICE TRANS. COMMUN., VOL.E82–B, NO.12 DECEMBER 1999

2095

PAPER

Stability Analysis for Global Performance of Flow Control
in High-Speed Networks Based on Statistical Physics
Masaki AIDA† and Kenji HORIKAWA†† , Members

SUMMARY This paper focuses on ﬂow control in high-speed
and large-scale networks. Each node in the network handles its local traﬃc ﬂow only on the basis of the information it knows. It ispreferable, however, that the decision making of each node leads
to high performance of the whole network. To this end, the relationship between local decision making and global performance of
ﬂow control is the essential object. We propose phenomenological models of ﬂow control of high-speed and large-scale networks,
and investigate the stability of these models.
key words: high speed network, ow control, autonomous dis-tributed system, throughput, Fokker-Planck equation
1.

Introduction

This paper investigates performance and stability of
ﬂow control schemes in high-speed and large-scale networks. In order to clarify our motivation, this section
states the characteristics of high-speed and large-scale
networks and states issues concerning the frameworksof network control associated with them.

earth and the moon.
This means that it is impossible to exert timesensitive control based on collecting global information about the network. So, the frameworks of timesensitive control, in a large-scale and high-speed network, are inevitably autonomous distributed systems.
1.2 Issues in Large-Scale Networks
Let us consider an approach analyzing a large-scale network. If we model it based on traditional queueing network theory, a lot of degrees of freedom of the network
system are included in the model. For example, in order to derive the global performance from the state ofeach node, it is necessary to obtain the joint probability
distribution of the states of all nodes in the network. If
the distribution is easily determined, there is no problem. However, it is actually impossible to obtain the
distribution of a lot of nodes including their correlated

1.1 Issues in High-Speed Networks
In a high-speed network, propagation delay becomes adominant factor in the transmission delay. This is because light speed is a non-scaling factor, and is the samefor high-speed networks. Therefore, at a given time, alarge amount of data is being propagated on links inthe network (Fig. 1). The amount of this data is characterized by delay-bandwidth product, that is, the propagation distance times transmission rate. Therefore, inhigh-speed and/or long-distance transmission, there isa larger amount of data on links than in nodes. Figure 2
shows an example of how much data can be on a link.
Let us consider the situation involving data transmission between two nodes, with a distance between themof 1 km and a link speed of 1 Mbps. If transmission
speed increases to 1 Gbps, the data amount on the link
is equivalent to 103 km on a 1-Mbps link. In addition,
if its transmission speed increases to 1 Tbps, the dataamount is equivalent to 106 km on a 1-Mbps link. Thisdistance is about 2.5 times the distance between the
Manuscript received December 11, 1998.
Manuscript revised June 24, 1999.
†
The author is with Traﬃc Research Center, NTT Advanced Technology Corporation (NTT-AT), Musashino-shi,
180-0006 Japan.
††
The author is with NTT Information Sharing Platform
Laboratories, Musashino-shi, 180-8585 Japan.

Fig. 1

Fig. 2Large delay-bandwidth product.

An example of large delay-bandwidth product.

IEICE TRANS. COMMUN., VOL.E82–B, NO.12 DECEMBER 1999

2096

eﬀects.
When all nodes are mutually independent, as incases of BCMP theorem [2], the probability distribution is obtained in a simple product-form. These cases,
however, are strongly restricted situations. Conversely,
if we assume the independence among nodes, the interesting behavior of global performance can not be
described. Thus, the approach directly describing all
degrees of freedom is diﬃcult to handle.
In general, when a system is in equilibrium, thesystem is characterized by a few order parameters. The
kind of these order parameters is much smaller thanthe degrees of freedom of the whole system. Therefore
we adopt an approach which directly describes the behavior of a global performance as an order parameter.
Then the eﬀect of the behavior of each node is reﬂected
as ﬂuctuation of the global performance.
1.3 Framework of Flow Control in High-Speed andLarge-Scale Networks
This paper focuses on ﬂow control in networks in which
nodes handle their local traﬃc ﬂow themselves based
only on the information they know. It is, of course,
preferable that the decision making of each node leads
to high performance of the whole network. In ﬂow control, we adopt the whole throughput of a network as a
global performance measure.
To achieve such coherent ﬂow control, the relationship between the local decision and the global performance of ﬂow control is the essential object. We
propose phenomenological models of ﬂow control for alarge-scale and high-speed network, and investigate thestability of these models. These models assume thatthe whole throughput of a network has a Markovian
property, and they use the technique of Ω-expansion of
the master equation. Based on these models, we have
tried to connect the whole throughput and its ﬂuctuation caused by the local decision making of each node.
2.

Background

2.1 Methods in Statistical Physics and Thermodynamics
Data in a network is in one of the following two states:in a node, or on a link. We deﬁne the throughput of
a network at time t as how much data is being propagated on the network [3], [6]. So, the throughput of a
network is the amount of data on the links at time t,
and is denoted by XE (t). Alternately, the average link
utilization of the whole network, that is, the normalized value of the amount of data on the links, may also
be thought of as the throughput of the network. We
denote it by XI (t).
These two quantities are categorized from the thermodynamical point of view as follows:• Extensive Quantity
Consider an equilibrium system made by combining n identical equilibrium subsystems. Whenthe thermodynamical quantity of the large system
equals n times the quantity of one small subsystem, the quantity is called an extensive quantity.
XE (t) is a quantity of this type.
• Intensive Quantity
In the same situation, when the thermodynamicalquantity of the large system equals the quantityof one small subsystem, the quantity is called an
intensive quantity. XI (t) is a quantity of this type.
These quantities are homogeneous functions of degreeΩ1 and Ω0 , respectively, where Ω is system size and,in this case, denotes the maximum amount of data onlinks in the whole network.
The behavior of throughput XE is determined by
the behaviors of all nodes in the network. However,when the network size Ω is large, there are many nodes
in the network and it is almost impossible to describe
their behavior, including the interactions among all
nodes. Therefore, we assume that the inﬂuence of the
behavior of a node is reﬂected in a small ﬂuctuation of
XE . So we regard XE as a random variable.
We deﬁne the probability density function PE (x, t)
such that the probability of the throughput XE (t) being
in x ≤ XE (t) < x + dx at time t is PE (x, t) dx. We
assume XE (t) as a Markovian random variable, anddescribe the temporal evolution of PE (x, t) using themaster equation:∂PE (x, t) =W (x − r, r, t) PE (x − r, t) dr∂t− W (x, r, t) PE (x, t) dr,(1)
where W (x, r, t) denotes the transition rate fromXE (t) = x to x + r at time t. Our interest is in thesituation for large networks, so we assume the following conditions:Ω  1, E(XE )  1.

(2)

2.2 Ω-Expansion of the Master Equation
Let us deﬁne moments of W (x, r, t) with respect totransition r asCn (x, t) := r n W (x, r, t) dr.
(3)
If all moments Cn (x, t) exist for all n, master equation(1) is written as
 n
∞

∂
(−1)n ∂PE (x, t) =Cn (x, t) PE (x, t).
∂t
n!
∂x
n=1
(4)

AIDA and HORIKAWA: STABILITY ANALYSIS FOR GLOBAL PERFORMANCE OF FLOW CONTROL

2097

This is called Kramers-Moyal expansion [4], [5].
In order to describe the relationship between thetemporal evolution of XE and its ﬂuctuation, we useΩ-expansion (size expansion).
When the network size Ω is very large and XE isnear equilibrium, we assume
E(XE ) = O(Ω), V ar[XE ] = O(Ω),(5)

where V ar[·] denotes the variance. Here, we introduce
normalized throughput XI (intensive quantity) as
XE (t),(6)Ω
which is independent of the network size. The varianceof XI is then expressed as
XI (t) :=V ar[XI ] = ε2 V ar[XE ],

(7)

−1

where ε := Ω . We deﬁne the probability that XI (t)is in x ≤ XI (t) < x + dx at time t as PI (x, t) dx. From
the normalization condition of PI (x, t), we obtain
PI (x, t) = Ω PE (Ωx, t),(8)

The master equation of Q(ξ, t) is then expressed as


 r
1
∂
Q(ξ, t) =dr dη
ds iξ eiεξs∂t
2π
0× Q(ξ − η, t) ω(η, r, t).
(14)
We assume that the solution Q(ξ, t) of masterequation (14) has the following form:∞ (iξ)nqn (t) . (15)Q(ξ, t) = exp q(ξ, t) = expn!
n=1Here, qn (t) is the n-th cumulant of XI .
From Eq. (14), we can derive the temporal evolution of cumulants qn (t)s. q1 (t) and q2 (t) respectively
correspond to the average and the variance of PI (x, t).
Expanding q1 (t) up to O(ε) with respect to the powerof ε, we deﬁne y(t), u(t) as
q1 (t) = y(t) + ε u(t) + O(ε2 ).

Similarly, expanding q2 (t) up to O(ε) with respect tothe power of ε, we deﬁne v(t) as
q2 (t) = ε v(t) + O(ε2 ).

under scaling law (6). The transition rate, w(x, r, t),from XI = x to x + r at time t isw(x, r, t) = ε W (Ωx, r, t).

(9)∂
PI (x, t) =
εn−1
∂t
n!
n=1
 n
∂
×
cn (x, t) PI (x, t),
∂x

(10)

(17)

The temporal evolution equations of y(t), v(t), u(t) are,as shown in [5], expressed as

The physical meaning of this scaling is that the transition r in extensive quantity XE has Ω times the opportunity to occur as the same transition in normalized
intensive quantity XI . Using these scaling rules, werewrite master equation (1) in the following KramersMoyal expansion form:
∞

(−1)n

(16)

∂y(t)
= c1 (y, t),∂t∂v(t)= 2 c1 (y, t) v(t) + c2 (y, t),
∂t
1
∂u(t)
= c1 (y, t) u(t) + c1 (y, t) v(t),
∂t
2

(18)
(19)
(20)

where cn and cn denote the ﬁrst and second derivatives,
respectively, of cn with respect to y. Note that thetemporal evolution of the cumulants is determined by
the moments of the transition rate.

where cn (x, t) denotes the n-th moment of w(x, r, t)with respect to transition r, as
(11)cn (x, t) = r n w(x, r, t) dr.

3.

Equation (10) is called the Ω-expansion [4], [5] of masterequation (1).

Our network model has simple lattice topology and atorus boundary, that is, a closed Manhattan Street network (Fig. 3(a)). There are 400 nodes in the networkand they have a 20 × 20 lattice conﬁguration. All nodes
have two incoming links and two outgoing ones. Twolinks of a node correspond to the vertical and horizontal
directions. For simplicity, we assume an AsynchronousTransfer Mode (ATM) network in which the data unit,
the cell, has a ﬁxed length in bits.
All nodes have a switching function such that theincoming cells are switched into the vertical or horizontal outgoing direction with the probability 1/2 basedon a Bernoulli trial. Switched cells are stored in the

2.3 Cumulant Expansion and Temporal Evolution of
Cumulants
We deﬁne the characteristic function of PI (x, t) as
Q(ξ, t) := PI (x, t) eiξx dx,(12)and the characteristic function of the transition rate asω(ξ, r, t) := w(x, r, t) eiξx dx.
(13)

Flow Control Schemes and Simulation Results

3.1 Network and Node Models

IEICE TRANS. COMMUN., VOL.E82–B, NO.12 DECEMBER 1999

2098

ronment. We therefore give all nodes the same
rule. Note that although the rule is the same, eachnode is in a diﬀerent state at some time t and thus
the behavior of each node at t is diﬀerent. We areinterested in the behavior of global performance
under the situation that

Fig. 3

Network and node models.

output buﬀer of their direction (Fig. 3(b)).
All nodes are capable of receiving and sending rate
signals. When a node receives a rate signal from downstream, the node adjusts its transmission rate to downstream according to the rate signal. In addition, allnodes can signal the allowed transmission rate to the
upstream nodes by using a rate signal. Each node autonomously decides how much rate should be speciﬁedby using a rate signal. Diﬀerent rate control schemes
are applied in accordance with the diﬀerent processes
used to determine the speciﬁed rate.
Here, we state supplements of characteristics in the
network model.
• Network topology has high symmetry
Although the homogeneous topology of the network model may look un-realistic, the reason is asfollows. If we can identify the bottleneck points,
in performance, in the network model, then we canfocus on the part of the network model including
the bottleneck points. However, if the symmetryof the network topology is high and the bottleneck
points in performance are not explicitly identiﬁed,
the above model focusing on part of the network
can not be applied. Actually, any parts can be bottlenecks with the same probability. Let us consider
the situation when accidental congestion occurs at
a certain node. Then we are interested in the behavior of the local congestion, whether– it grows and causes deterioration of thewhole network performance through interaction among nodes, or– it remains a local phenomenon and diminisheswith time.
In order to evaluate the behavior without othercomplicated eﬀects, we adopt the above simple network model. This is the reason why we do not takehierarchical network models (or, for example, additional ﬂow control between terminals such as TCP)
into consideration.
• All nodes obey the same rule
There are no special nodes such as one to control
the whole network. More speciﬁcally, such nodes
can not exist under the high-speed network envi-– all nodes behave autonomously by the same
rule, and– there is no special node controlling the wholenetwork.
• Network size and boundary condition
There is no essential reason that the number ofnodes is 400 with a torus boundary. It is only necessary that there be as many nodes as possible, andthe nodes on the network boundary behave under
the same rule and conditions as the other ordinary
nodes. The reason for a closed network is that the
number of cells in the network should be invariant.
This is required, to compare the diﬀerent ﬂow control schemes under the same conditions, as shownin below. Note that if we focus on a part of the
network, it can be regarded as an open network inwhich the number of cells is not invariant.
3.2 Two Flow Control Schemes
In this subsection, we show two diﬀerent rate control
schemes for the diﬀerent processes used to decide thespeciﬁed rate.
One is for specifying rates using a rate signal fromdownstream, i.e., a rate-based ﬂow control. We callit Rate Driven Control (RDC). Ordinarily, each node
sends a rate signal at every deterministic time interval. Let D be the round-trip time between adjacentnodes. For simplicity, we set the time interval between
consecutive rate signals to D.
Consider a node downstream from link i. The nodecalculates the following quantity, at every time periodD:i (t) = Ni (t) + {Rrin (i, t − D) − Rrout (i, t − D)} D,N(21)
where Rrin (i, t) denotes the speciﬁed rate to send upstream at time t, Rrout (i, t) denotes the speciﬁed ratereceived from downstream at time t, and Ni (t) denotes
the queue length of the buﬀer. Then the node calculates the new speciﬁed rate as
Rrin (i, t) = Rrin (i, t − D) +i (t)
αL − N,D(22)

where L denotes the buﬀer capacity, and α is a threshold parameter such that αL denotes the threshold of
the buﬀer. Ordinarily, each node sends a rate signal at
every deterministic time interval, D. However, in caseof emergency, that is, during buﬀer overﬂow/underﬂow,

AIDA and HORIKAWA: STABILITY ANALYSIS FOR GLOBAL PERFORMANCE OF FLOW CONTROL

2099

Fig. 4Behavior of normalized throughput.

each node can immediately send a rate signal.
The other scheme speciﬁes a rate based only on thequeue length of the nodes’ own buﬀers, i.e., a queuethreshold based ﬂow control. We call it Queue Driven
Control (QDC). Ordinarily, each node sends a rate signal at every deterministic time interval, D. The specifying rate is determined by the queue length of its buﬀers.
Each node makes an eﬀort to ensure that the queue
length is equal to the predeﬁned threshold in its buﬀers.
Let Rqin (i, t) be the speciﬁed rate to send upstream attime t:Rqin (i, t) = Rqin (i, t − D) +αL − Ni (t)
.
D

(23)

As with RDC, in case of emergency, that is, during
buﬀer overﬂow/underﬂow, each node can immediately
send a rate signal.
3.3 Experimental Evaluation of Stability
Let us investigate stabilities of the above two ﬂow control schemes, RDC and QDC, by using a Monte Carlosimulation. The environment of our simulation is asfollows:• bandwidth and distance between adjacentnodes
Each horizontal link in Fig. 3(a) is 300 Mbps and
20 km in length, and each vertical link is 600 Mbps
and 2 km in length. The links are assumed to be
WANs and LANs, respectively. Note that the horizontal links are the bottle-necks of the traﬃc ﬂow.
• buﬀer capacity
All nodes have the same capacity, 300 cell places.
The threshold parameter α is set to be 1/2.
• cell overﬂow
In order to maintain the network load, the discarded cells, when no buﬀer space is left, are again
added to network nodes chosen at random. Therefore, the total number of cells in the network isinvariant.
• total number of cells

There are 40,000 cells in the network. This corresponds to about 160% of the total capacity of alllinks. This implies that the load of the network is
high.
Based on the above environment, we compare thethroughput and queue length of the two networks, one
controlled by RDC and the other by QDC.
At the initial time t = 0, both networks are nearequilibrium states. Figure 4 shows the behavior of theaverage and the variance of the normalized throughputof RDC and QDC. These ﬁgures show that RDC isstable but QDC is not, and the variance is remarkably
increased for QDC.
Figure 5 shows the queue lengths of buﬀers in thenodes in the network. These data are from a typicalsample path of the simulation. Each pixel corresponds
to a node and the conﬁguration is the same as a network topology. The queue length of a node is denotedby the color density. “Step” denotes simulation time,“Link” denotes the normalized throughput, “Queue”
denotes the average queue length, and “Loss” denotes
the number of lost cells (running total).
The ﬁgure shows that the queue length remains homogeneously distributed for RDC. However, for QDC,
node queue lengths go to basically two states: very long
and almost empty. The nodes with very long queue
lengths are clustered. This clustering process corresponds to the process of decreasing normalized throughput, shown in Fig. 4.
3.4 Experimental Evaluation of Robustness
Next, let us investigate robustness of the performancewith respect to increasing network load. Simulation
environments are almost the same as the above model,
but the number of cells in the networks are increased
with time. The number of cells in the network is zero at
the initial state. Then, ﬁve cells are added to network
nodes chosen at random, every one simulation time.
The results obtained in the simulation are shownin Fig. 6. Note that the number of cells increases as

IEICE TRANS. COMMUN., VOL.E82–B, NO.12 DECEMBER 1999

2100

Fig. 5Behaviors of queue length in each node.

time progresses. As shown in Fig. 6(a), most of cells remain on the link between 0 and 4,000 cell times because
the network is not busy, and incoming cells are immediately sent out. After 5,000 cell times, the throughput
becomes constant at about 92%. This means that cells
begin to be stored in buﬀers.
After 15,000 cell times (after the number of cells in
the network exceeds 75,000), the throughput of the network controlled by the QDC begins to drop. Finally, its
throughput is lowered to nearly 0 and the control fails.

On the other hand, even at 25,000 cell times (numberof cells is 125,000), RDC maintains high throughput
without failure. In addition, RDC causes hardly any
cell loss, as shown in Fig. 6(b). The throughput RDC
begins to drop and the control is fails at 30,000 (150,000
cells). Since the total capacity of the buﬀer for the bottle neck links is 120,000 cells, however, it is safe to saythat RDC can endure in higher load conditions and canmaintain normal network operation.
Figure 7 shows the behaviors of the queue lengths

AIDA and HORIKAWA: STABILITY ANALYSIS FOR GLOBAL PERFORMANCE OF FLOW CONTROL

2101

Fig. 6

Behaviors of throughput and the number of lost cells.

of buﬀers in the nodes in the networks as the same in
Fig. 5. Since the number of cells in the network increases as the simulation advances, it becomes blackish
overall. It is possible to see that the areas of congested
nodes spread in the QDC. On the other hand, in the
RDC it can be understood that there are some congested nodes but they are not clustered, and the expansion of congestion is minimized to prolong control
without cell loss.
4.

Linear Relaxation Model of Rate

4.1 Linear Interaction Model
This subsection proposes a phenomenological model for
RDC. Let M (i, t) be the throughput of link i, that is,
the number of cells on link i at time t. We assume theutilization of the network is suﬃciently high.
When Rrin (i, t) = Rrin (i, t−D), no change occurs. If
in
Rr (i, t) = Rrin (i, t − D), M (i, t) changes to 12 DRrin (i, t).
The time for achieving 12 DRrin (i, t) is independent of its
value, and instead depends only on the round-trip timeD, that is,
M (i, t + D) − M (i, t)
1= D (Rrin (i, t) − Rrin (i, t − D)).
2

Let the normalized throughput be XI = y0 when the
network is in equilibrium. From Eqs. (6), (25) and (26),
we have
XI (t + ∆t) − XI (t)XI (t) − y0∆t + O(ε) ∆t,= −Ddetermines the transition rate of M (i, t). When M (i, t)
is near equilibrium M0 (i, t), and since we regard D asa time constant, we model∂
M (i, t) − M0 (i, t)
M (i, t) = −
,
(25)
∂t
D
as linear relaxation.
The throughput of the network is then denoted byXE (t) =M (i, t).
(26)

(27)

for small ∆t. Since y(t) in Eq. (16) means the leading
term (O(1) term) in the ﬁrst cumulant (the average) of
XI , we focus only on O(1) term in Eq. (27) and replace
XI (t) with y(t). Then we have∂y(t)y − y0
=−
.
∂t
D

(28)

From the temporal evolution equation (18), we havec1 (y, t) = −y − y0,D(29)

as shown in [1]. In addition, we assume
c2 (y, t) = const. = br ,(24)

This means the diﬀerence of Rrin (i, t) and Rrin (i, t − D)i

4.2 Gaussian Approximation near Stable Equilibrium
State

(30)

like Brownian motion. We set the initial distributionof XI asPI (x, 0) = δ(x − x0 ).

(31)

Then non-equilibrium distribution of XI near the equilibrium state is1
PI (x, t) = 
2πεv(t)× exp −(x − y0 − (x0 − y0 ) e−t/D )2,
2 ε v(t)(32)

IEICE TRANS. COMMUN., VOL.E82–B, NO.12 DECEMBER 1999

2102

Fig. 7

Behaviors of queue length in each node.

where, from v(0) = 0,v(t) =br D(1 − e−2t/D ),
2

= −
(33)

this is the solution of Eq. (19). The non-equilibrium
distribution (32) is the solution of the Fokker-Planck
equation [7]∂PI (x, t)
∂t

1 ∂2
∂
c1 (x, t) + εc2 (x, t) PI (x, t),
∂x
2 ∂x2
(34)

obtained by truncating Eq. (10) after O(ε2 ).
Figure 8 shows an example of temporal evolutionof the solution (32), where we set y0 = 0, x0 = −1,
D = 1, and ε = br = 1 for simplicity. Initially, thesolution is around x0 , and afterward it moves to y0 as
t increases. The movement is slowing down when the

AIDA and HORIKAWA: STABILITY ANALYSIS FOR GLOBAL PERFORMANCE OF FLOW CONTROL

2103

Fig. 9Fig. 8Behavior of the solution of Gaussian approximation.

average of the distribution (32) approaches the equilibrium y0 . The distribution (32) approaches, in accordance with increasing t,
lim PI (x, t) = √t→∞(x − y0 )2
1
exp −
,
ε br Dπεbr D(35)

and it is independent of the initial state XI = x0 .
5.

Nonlinear
LengthRelaxation

Model

of

Queue

5.1 Redeﬁnition of a Network Performance MeasureThis subsection shows a model of the QDC scheme.
For QDC, the behavior of the throughput is directly
controlled by the queue length in each node through(23). Since the time derivative of the throughput corresponds to the second time derivative of the queue
length, the equation describing temporal evolution of
the throughput is, therefore, the second diﬀerential
equation with respect to time. This implies that the solution of the equation contains some oscillation modes.
It is thought that such oscillation causes the deterioration of throughput shown in Fig. 5. Unfortunately, thetechnique based on the master equation is not applicable to this type of equation.
To avoid this, we want to focus on the trend of
the deterioration process of throughput without oscillation. To this end, we describe the behavior of the queue
length in each buﬀer directly using a phenomenologicalmodel. Since, in our network model, there is a constantnumber of cells in the closed network, if we know the
behavior of queue lengths, we can ﬁnd the throughput,that is, the number of cells on links. For simplicity,
we introduce the parameterized buﬀer shown in Fig. 9.
Let us assign 0 to the point of half the buﬀer capacity, that is, just at the threshold. The fully occupied

Buﬀer model.

and emptied points of the buﬀer are assigned + 12 L and− 12 L, respectively. Values of Ni (t) are denoted by thisbuﬀer.
We deﬁne the total number of cells in all buﬀers
in the network at time t, denoted by NE (t), and thenormalized value of NE (t) denoted by NI (t), that is,Ni (t),(36)
NE (t) :=i

NI (t) := NE (t)/ΩN ,(37)

where ΩN is system size and, in this case, denotes thenumber of buﬀers in the network. NE (t) and NI (t)are categorized into extensive and intensive quantity,
respectively, same as XE (t) and XI (t).
We deﬁne the probability density function ϕ(χ, t)such that the probability of the normalized queue
length NI (t) being in χ ≤ NI (t) < χ + dχ at time tis ϕ(χ, t) dχ. We assume NI (t) as a Markovian randomvariable, and describe the temporal evolution of ϕ(χ, t)using a master equation. Truncating the Ω-expansionof the master equation after O(+2 ), the temporal evolution equation of ϕ(χ, t) is obtained as the following
Fokker-Planck equation:∂ϕ(χ, t)
∂t
= −

∂
1 ∂2
c1 (χ, t) + +c2 (χ, t) ϕ(χ, t),
∂χ
2 ∂χ2
(38)

where + := Ω−1
N , and c1 and c2 are drift and diﬀusion
coeﬃcients, respectively.
5.2 Minimal Nonlinear Interaction ModelAs shown in Sect. 3.3, the QDC scheme is unstable inequilibrium at high throughput. In addition, complicated behavior such as clustering implies that there isnonlinear interaction among the nodes. It is thus necessary to consider a nonlinear diﬀerential equation.
Let us consider the minimal nonlinear interaction
describing the QDC scheme. We consider the FokkerPlanck equation (38) and assume the nonlinearity is
only in the drift coeﬃcient c1 (χ, t), and the diﬀusion
coeﬃcient c2 (χ, t) describes the simple Gaussian white
noise.

IEICE TRANS. COMMUN., VOL.E82–B, NO.12 DECEMBER 1999

2104

There are three equilibrium states of NI (t):1
1
NI (t) = 0, NI (t) = + L, NI (t) = − L.
2
2

(39)

As shown in Sect. 3.3, NI (t) = 0 is unstable.
We assume c1 (χ, t) is to have the following power
series expansion form:
c1 (χ, t) =∞an (t) χn .

(40)

n=0

Since the buﬀer threshold, that is, the target queue
length, is at half the buﬀer capacity, an (t) = 0 for all
even n. Thus, the terms of the even degrees of χ are outof consideration from the symmetry of the buﬀer. In
addition, if XI is near equilibrium, i.e., χ 0, the termsof the higher degrees of χ are negligible. Therefore, thedrift coeﬃcient (40) is approximately denoted asc1 (χ, t)

a1 (t) χ + a3 (t) χ3 ,(41)

near equilibrium χ 0. Since NI (t) = 0 is unstable, let
us choose a1 (t) > 0.
If we assume the simplest non-linear interaction
among nodes, c1 (χ, t) is assumed to have the form (41)for all χ. Since χ = 0, ±L/2 are equilibrium points andc1 (χ, t) = 0 at the equilibrium points, we havec1 (χ, t) = γ χ − g χ3 ,ϕ(χ, t)(42)

Behavior of the scaling solution.
−3/2

1
g
1 − χ2 (1 − e−2γt )γ2πτχ2× exp −,2 τ [1 − (g/γ) χ2 (1 − e−2γt )](45)√

where τ is the scaling time such that

where γ, g > 0 and
g = 4γ/L2 .

(43)

It is natural to consider that γ is determined by theround-trip time D.
On the other hand, since the diﬀusion coeﬃcientc2 (χ, t) is assumed to describe the simple Gaussian
white noise, we have c2 (χ, t) = const. = bq .
The minimal nonlinear Fokker-Planck equation isthen obtained by∂ϕ(χ, t)∂t∂+bq ∂ 2=(−γ χ + g χ3 ) +
∂χ
2 ∂χ2

Fig. 10ϕ(χ, t).
(44)

In addition, we denote the ﬁrst and second cumulantsof ϕ(χ, t) as Eqs. (16) and (17), after we replace ε with+. Then, the temporal evolution equations of the cumulants of ϕ(χ, t) are obtained as Eqs. (18)–(20).
5.3 Nonlinear Scaling and Scaling Solution
Let the initial distribution of the nonlinear FokkerPlanck equation (44) be ϕ(χ, 0) = δ(χ). This is theunstable equilibrium. According to the nonlinear scaling theory [8], the solution of (44) is obtained approximately asτ=+bq 2γt(e − 1).
2γ(46)

Solution (45) is called the scaling solution.
This scaling solution shows an interesting behavior. Figure 10 shows a temporal evolution of the scaling
solution, where we set L = 2, γ = g = 1, and + = bq = 1
for simplicity. Initially, the solution has a single peak,
and afterward it splits into double peaks in accordancewith increasing t. These peaks signify the very long andthe almost empty states of queue length. Note that we
can recognize that the double-peak state corresponds tothe ﬁnal state of Fig. 5 for QDC. Therefore, the emerging double peak means that the deteriorating processof throughput is caused by clustering of the busy nodes
and their spread.
Here, we state a supplemental comment about theabove interpretation. Note that the ﬁnal state in Fig. 10
does not mean that each node becomes busy or idle
alternatively in the ﬁnal state of QDC in Fig. 5. The
ﬁnal state in Fig. 10 means that state of the network
becomes one of the following states:• all nodes in the network are congested, or• all nodes in the network are idle.
Since our network model is closed and the number ofcells in the network is constant, it is impossible for thestate of the network to become one of the above two

AIDA and HORIKAWA: STABILITY ANALYSIS FOR GLOBAL PERFORMANCE OF FLOW CONTROL

2105

states. If we consider a very large closed network andwe focus on a part of the network as an open network,that partial network must become one of the above
states. We should therefore interpret Fig. 10 as a clustering process of busy state nodes.
5.4 Anomalous Fluctuation and Behavior of Throughput
We are interested in the time the throughput of the
network deteriorates. This is interpreted as the time
when double peaks of the solution emerge. Let thevariance of the solution be, initially, O(+). However, thevariance becomes anomalously large and reaches O(1)as t increases. This is called the Anomalous Fluctuation
or Anomalous Enhancement of Fluctuation, and often
appears in state transition from unstable equilibrium
to stable equilibrium. Temporal evolution (19) implies,
as shown in [8], [9],v(t) ∝

ẏ(t)
,
ẏ(0)

(47)

where ẏ(t) denotes the time derivative of y(t). Notethat y(t) and v(t) describe not the throughput but thequeue length in this QDC case. However, since the total
number of cells in the network is constant, the variancev(t) is the same as that of the throughput, and ẏ(t) is
the negative value. Therefore, we can also say that (47)is reasonable for describing the throughput.
Actually, the QDC case shown in Fig. 4 shows thetime when the maximum variance corresponds to thesteepest change in the average throughput.
The time the throughput of the network deteriorates corresponds to the time the variance becomeslarge and reaches O(1). The physical meaning of thisis that a small ﬂuctuation grows and causes the average value to change. The time for it, t0 , is called on-set
time [8], and is obtained as
t0
6.

1
2γ 2
log
.
2γ
+bq g

(48)

Conclusion

We have described two ﬂow control schemes, RDC and
QDC, in high-speed and large-scale networks as autonomous distributed systems. In both schemes, eachnode in the network handles its local traﬃc ﬂow only
on the basis of the information it knows. It is preferable, however, that the decision making of each nodeleads to high performance of the whole network. To
this end, the relationship between local decision making
and global performance of ﬂow control is the essential
object.
We are especially interested in the behavior of thelocal congestion, whether it grows and causes deterioration of the whole network performance through interaction among nodes, or it remains a local phenomenon

and diminishes with time. Experimental results showinteresting behaviors of throughput for both schemes.
Performance of RDC is stable, but that of QDC is unstable. For QDC, the congested nodes are clustered.
To describe the behaviors, we proposed phenomenological models based on the Ω-expansion technique and gave physical interpretations. Behavior ofthroughput for RDC is described by a linear relaxation model of the Fokker-Planck equation. The solution of the equation is a stable Gaussian distribution.
For QDC, we have assumed that the behavior of queue
length can be described by a nonlinear Fokker-Planck
equation, and have applied the simplest nonlinear relaxation model. The solution of the equation describes
the clustering process of the congested nodes.
Main residual issues are listed as follows:• To verify the minimal non-linear model Eq. (42).
(Is Eq. (42) valid for a large χ?)• To obtain the relationship between the round-trip
time, D, and coeﬃcients in Eq. (44).
• To obtain necessary and/or suﬃcient conditions
which enable us to design a stable ﬂow control
scheme in high-speed networks.
Acknowledgement
The authors would like to thank anonymous reviewers
for their valuable comments and instructions.
References[1] M. Aida and K. Horikawa, “A study on stability analysis
for high-speed networks based on statistical physics,” Proc.
Symp. on Performance Models for Information Communication Networks, pp.1–7, Jan. 1997.
[2] F. Baskett, K.M. Chandy, R.R. Muntz, and F.G. Palacios,“Open, closed, and mixed networks of queues with diﬀerent classes of customers,” Jour. Ass. Computing Machinery,vol.22, pp.248–260, 1975.
[3] K. Horikawa, M. Aida, and T. Sugawara, “Traﬃc control
scheme under the communication delay of high-speed networks,” Proc. 1996 Int. Conf. on Multi-Agent Systems (ICMAS’96), pp.111–117, 1996.
[4] N.G. van Kampen, “A power series expansion of the masterequation,” Can. J. Phys., vol.39, pp.551–567, 1961.
[5] R. Kubo, K. Matsuo, and K. Kitahara, “Fluctuation and relaxation of macrovariables,” Jour. Stat. Phys., vol.9, pp.51–
96, 1973.
[6] D. Mitra, “Optimal design of windows for high speed data
networks,” Proc. INFOCOM ’90, pp.1156–1163, 1990.
[7] H. Risken, Fokker-Planck Equation: Methods of Solution
and Applications, Springer-Verlag, Berlin, 1989.
[8] M. Suzuki, “Passage from an initial unstable state to a ﬁnal
stable state,” Adv. Chem. Phys., vol.46, pp.195–278, 1981.
[9] M. Suzuki, “Anomalous ﬂuctuation and relaxation in unstable systems,” Jour. Stat. Phys., vol.16, pp.477–503, 1976.

IEICE TRANS. COMMUN., VOL.E82–B, NO.12 DECEMBER 1999

2106

Masaki Aida
received the B.S. and
M.S. degrees in Theoretical Physics from
St. Paul’s University, Tokyo, Japan, in
1987 and 1989, respectively, and received
the Ph.D. degree in Telecommunication
Engineering from the University of Tokyo,
Japan, in 1999. Since he joined NTT
Laboratories in 1989, he had been mainly
engaged in research on traﬃc issues in
ATM networks and computer communication networks until March 1998. He iscurrently a manager at Traﬃc Research Center, NTT Advanced
Technology Corporation (NTT-AT). His current interests include
traﬃc issues in communication systems. He received the Young
Engineer Award of IEICE in 1996. Dr. Aida is a member of the
Operations Research Society of Japan.

Kenji Horikawa
received the B.S.
and M.S. degrees from Tokyo Instituteof Technology in 1992 and 1994, respectively. In 1994 he joined NTT (Nippon
Telegraph and Telephone Corporation).
He has been engaged in research and development of computer networks.