## Abstract

We present a possible way to study large-size packet-switched networks capable of accounting for interactions between adjacent queues. The interaction between queues arises, because of the influence of the routing protocol on each switching decision, and the stochastic nature of packet lengths and inter-arrival times.Both the methodology and the analysis tools are adaptations of methods of statistical mechanics. The justification for their use lies in recent experimental evidence indicating that aggregate, core-network IP traffic, exhibits quasi-Markovian properties when the network is heavily loaded.In this paper, we present a general methodology and introduce approximations that greatly simplify the analysis. These approximations, are owing to the quasi-Markovian nature of the traffic and the large size of the network.

## 1. Introduction

Next‐generation telecommunication networks are likely to rely on an IP core network infrastructure. As a consequence, unlike today's Internet, future networks will be subject to much more demanding requirements. In order to provide operational guarantees, network owners need tools which enable the dimensioning of the core of packet-switched networks. It is customary for large network design to be primarily based on simulations. However, replacing theory by simulations is non-tenable owing to the anticipated size of such future networks.

A theory of traffic in packet-switched networks must be capable of predicting a number of quantities of interest: end-to-end latency, packet loss rate, and so on. Based on such metrics of network performance, we would then wish to quantify the external loading point of a large system of interconnected queues, at which the network changes its behaviour. For example, if packet loss rate is of interest, we wish to know the loading point at which this rate exceeds a given threshold. One of the important features of such a theory is that it must model all sources of stochasticity in the system of interconnected queues that form the core network. This implies that correlations between the states of queues, which have a knock-on effect on each other, must be modelled explicitly as interactions. Conventional queueing theory does not account for these correlations, and its mathematical formalism becomes cumbersome to compute when the number of queues is very large.

Constructing such a theory has attracted a lot of attention (e.g. Kelly 1991; Schwartz 1987; Graham & Méléard 1994; Vvedenskaya *et al*. 1996; Delcoigne & Fayolle 1999; Fayolle *et al*. 2001), but no definitive theory has yet emerged. In this paper we take a somewhat different approach in attempting to formulate a phenomenological theory of packet-switched networks. The chosen approach can be scaled to large numbers of queues, which is essential in modelling a large-scale core IP network. The methodology adopted here is an adaptation of the methods of statistical mechanics, which are ideally suited to the study of large systems in the presence of sources of stochasticity.

## 2. General formulation of teletraffic theory

We consider a network consisting of *N*_{R} routers and *N*_{G} gateways, the latter considered as aggregate end-users. Each gateway is connected to several routers (*n*_{G}) and each router is connected to several gateways (*n*_{R}). These average connectivity numbers can be replaced by probability distributions in general.

We assume that a routing protocol operates in the network, which makes routing decisions, while attempting to minimize some metric (e.g. end-to-end latency, number of hops, etc.). To be specific, we assume, in what follows, that an OSPF-like protocol, operates in the network and makes its decision on routing an individual packet by trying to minimize the end-to-end latency for the packet. This can be modelled as follows: first, we introduce a ‘partition function’,(2.1)where *τ*_{P} is the latency of the packet along some path *P*, *τ*_{id} can be considered as a ‘temperature’ using a statistical physics analogy or an ‘ideally desired latency’, which one would like to keep as small as possible. The summation is over all possible paths connecting the source and the destination. As the ‘temperature’ goes to zero *τ*_{id}→0 (correspondingly, *β*→∞) the only contribution which survives in the sum in equation (2.1) is that from the path with the minimal latency:

We assume that all transmission lines are full duplex (transmitting in both directions independently). In what follows, we postulate the following packet flow-processing model for a router. First, a packet arrives to an input FIFO buffer; next, its header is processed and a forwarding decision is made, and the packet is segmented into cells of some fixed size; following that, the packet is sent to the switching fabric (we assume for simplicity that there is effectively only one switching fabric, and all input lines are served on a round‐robin basis); subsequent to switching, the packet is assembled and sent to the output buffer.

As an output buffer of one router is transmitting to a certain input buffer of another router, we are going to consider these pairs of buffers as one strongly correlated system of queues. Thus, the latency along some path *P*_{αβ} (from gateway *α* to gateway *β*) consists of four parts:(2.2)where is the sum of the transmission time from gateway *α* to the adjacent router *i*_{α} and the waiting time in the input buffer on the *i*_{α}th router; *τ*_{(ij)} is the cumulative time spent in the output buffer of the *i*th router and the input buffer of adjacent router *j* as well as the propagation time between these; *τ*_{i} is the time required for forwarding, processing and switching on the *i*th router; and is the waiting time in the output buffer of the *i*_{β}th router, plus the propagation time along the line. As we are going to study the behaviour of the network near a congested state, we expect the second term in equation (2.2) to be the most significant. This subdivision of delay terms in groups is not unique, but is guided (with some hindsight) by the need to consider strongly correlated subsystems of queues separately.

The time spent in a queue is connected in a straightforward way to the length of the queue. So, the end-to-end transmission latency of a packet depends on the lengths of queues on all the routers along the path from source to destination. Lengths of queues on all routers (as an ensemble) are distributed according to some probability density function (PDF).

It is reasonable to expect that the system will exhibit the property of ‘self-averageness’ as the size of the network grows. This property means that observable quantities of the network do not depend on a particular configuration of queue lengths in the network, but on statistical characteristics of the whole set of queue lengths considered as an ensemble. As a result, we can average the latency over the queue lengths without changing the answer:(2.3)

The property of self-averageness can be explained as follows. If an observable can be represented as a sum of random variables *ξ*_{1}, *ξ*_{2}, …, *ξ*_{n} (not necessarily independent), for example, expression (2.2) for latency, then using the necessary and sufficient condition for the law of large numbers (cf. Gnedenko 1963) we can obtain the following restriction on the covariance matrix (this is a sufficient condition): In order to satisfy this relation, either the covariance matrix should have only a few non-zero elements in each row (e.g. a queue length should be correlated only with a few neighbouring queues), or the correlation should decay sufficiently fast with increased ‘distance’ (according to some metric, e.g. number of hops) between queues so that ∀*i*. The former case is realized in the framework of the mean-field approximation which we are going to introduce in §3.

In order to compute the averaging in equation (2.3) we need to know the PDF for queue lengths. Furthermore, equation (2.3) should also include an averaging over all possible network topologies from a given class, as well as an averaging of all external input traffic patterns for the network.

A number of the assumptions made in this section (e.g. protocol operation, router packet flow model, duplex lines, etc.), as well as the choice of metrics (e.g. latency, etc.) described, can be easily changed to consider the operation of a large number of variants of the network model adopted in this paper. The methodology, however, will remain the same.

## 3. P.D.F. of queue lengths

In this section we consider some possible ways of describing dynamical and static properties of the PDF of queue lengths.

Consider a joint PDF of all queue lengths of all routers at time *t*,(3.1)where *N*_{q} is the total number of queues in the network, which can range from tens of thousands to several millions. Methodologically, the most consistent way would be writing down a set of stochastic differential equations for interacting queue lengths (Langevin equations in the theoretical physics terminology). Then using standard field-theoretic methods (e.g. Zinn-Justin 1993) we could obtain the functional-integral representation of the joint PDF of all queue lengths. But this is a tremendously complex task and we postpone it for the future.

In this paper, we deploy a somewhat simplified and more phenomenological approach to the calculation of the PDF

It is natural to expect that some sets of queue lengths are more correlated than others. As it is quite difficult to write down a dynamical equation for equation (3.1), let alone solve it, we will try to exploit some kind of approximation scheme, the simplest possible one being the mean-field approximation of quantum field theory.

In our case the mean-field approximation means the following. We use for the PDF equation (3.1) the following representation:where is the PDF of a subsystem of strongly correlated queues. The choice of the set of subsystems {S_{j}} is made from heuristic considerations: types of correlations which are of most importance should be included in one subsystem. Other interactions (correlations) should be accounted for in a self-consistent way, that is, the dynamical equation for a subsystem PDF should account for the fact that it interacts with an ‘average field’ of all other subsystems (or ‘average’ network load).

Our next step is the derivation of a dynamical equation for a subsystem's PDF *p*(*ℓ*, *t*), **ℓ**≡{*ℓ*_{k}}_{k∈S}. In the core network, the relevant stochastic processes can be approximated by (quasi-)Markovian ones under heavy loading conditions. The experimental evidence for this is presented in recent papers by Cao *et al*. (2002) and Stepanenko *et al*. (2001). Here we restrict ourselves to giving some heuristic reasons. While traffic generated by individual processes might have self-similar properties and be highly correlated, in the network core a large number of such *incoherent* processes are mixed together and the theorem of large numbers must become applicable. Moreover, sequential packets corresponding to the same process are separated by a much bigger number of bits owing to the higher speed of the core network links than those at the edge of the network, thus decreasing the level of correlations significantly. Furthermore, in a *heavily loaded network* (*congested or pre-congested states*), the queueing and switching processes in the core network routers will significantly modify the underlying stochastic process. Essentially, the temporal gaps between successive packets will be determined by the lengths of the non-empty, interacting queues and the switching fabric operation. Deviations from Markovian properties are small, and getting smaller as the capacity of links is becoming larger and the load heavier. If necessary, these deviations can be accounted for in the framework of perturbation theory.

If the underlying stochastic processes are Markovian, we can write down the Chapman–Kolmogorov equation,(3.2)where *w*(*ℓ*, *t*|**ℓ**′, *t*′) is the transition probability density from the state with the queue length *ℓ*′ (occupancy) at time *t*′ to the state with the queue length *ℓ* at time *t*.

The next step is to derive a differential equation for this PDF, starting from the Chapman–Kolmogorov equation in the form of equation (3.2), and a set of conditions which slightly differs from the one normally used in derivation of the differential Chapman–Kolmogorov equation (or Fokker–Planck equation in the case of a continuous stochastic process; e.g. Gardiner 1994). Here, we will use a general notation, denoting by bold letters ** x**,

**, … a whole set of random variables.**

*y*As we are dealing with the dynamics of queue lengths, we consider time increments Δ*t*, which are much smaller than the characteristic time-scale on which the PDF changes significantly, but big enough so that a sufficiently large number of packets arrives to (and departs from) a queue in order to consider the stochastic process as a Markovian one (see Cao *et al*. 2002; Stepanenko *et al*. 2001) on the time-scale Δ*t*.

We assume the stochastic process to be continuous (i.e. without jumps) which means (see Gardiner 1994) that for any *ϵ*>0, we haveor, in a slightly different form (as Δ*t* is not infinitely small in the case at hand),(3.3)

This is certainly true in our case, because all transmission lines have finite capacities and all switching devices have finite throughput. This means that *w*(**x**,*t*+Δ*t*|**y**,*t*)=0 if |**x**−**y**|>*C*Δ*t*, where *C* is the maximum line capacity (or switching throughput). Hence, we have (putting *ϵ*(Δ*t*)=*C*Δ*t*)

The last relation is even stronger than equation (3.3).

Introduce the following moments (per unit time):(3.4)

The first moment *m*_{i}^{(1)}(**y**, *t*) is the average rate of changing of the random variable *y*_{i} under the condition that the system was in the state **y** at time *t*. The second moment *m*_{ij}^{(2)}(**y**, *t*) is the covariance matrix per unit time, and so on. For a Markovian process we have (Gardiner 1994)(3.5)(3.6)Because in our case, Δ*t* is not infinitely small, the higher moments might be non-zero (while probably small). The actual values of the moments need to be determined experimentally (or some models should be constructed based on experimental data).

We assume that the PDF *p*(**x**, *t*) and the transition PDF *w*(**x**, *t*′|**y**, *t*) are defined in some region *R* of the parameter space (**x**, **y**∈*R*), and consider the time evolution of the expectation of a function *f*(**y**), which is infinitely differentiable in *R*:(3.7)

Performing all the usual steps as in the derivation of the Fokker–Planck equation (Gardiner 1994) we obtain(3.8)where we have introduced the ‘slow’ derivative (on a time-scale Δ*t*),(3.9)and Δ*t* is chosen as described earlier. Integrating equation (3.8) by parts (as many times as required for each term) we obtainwhere ∂*R* is the boundary of the parameter region *R*, {*n*_{i}(**x**)}≡**n**(**x**) is outward pointing normal to ∂*R*, and we have introduced the following set of ‘currents’:

Finally, as the function *f*(**x**) is arbitrary, we arrive at the following differential equation:(3.10)where can be interpreted as a probability current, and the boundary conditions imposed on it are(3.11)

The relation (3.10) is a modified Fokker–Planck equation, and it can be rewritten in terms of the transitional PDF as well:(3.12)

Note the similarity of equation (3.10) with the Kramers–Moyal expansion (Kramers 1940; Moyal 1949; Gardiner 1994; Risken 1989). More precisely, equation (3.10) can be considered as a generalization of the Kramers–Moyal expansion: replacing the slow derivative with the normal one in equation (3.10), we obtain the latter.

Considering equations (3.5), (3.6) and (3.9) in the limiting case Δ*t*→0, we obtain the ordinary Fokker–Planck equation:

The remaining terms are of order o(Δ*t*), namely, these are the terms with *n*≥2 in equation (3.9) and the terms with *k*≥3 in equation (3.10). They correspond to deviations from the Markovian case and should be treated in the framework of perturbation theory. The former (i.e. the terms with higher derivatives in time) are small because of the choice of Δ*t*, and the latter (i.e. the higher moment's terms) are small because the process is approaching a Markovian process at the relevant time-scales. The last statement is based on experimental observations (Cao *et al*. 2002; Stepanenko *et al*. 2001). Note that if at least one of the higher moments is of order o(Δ*t*), then the Pawula theorem (Risken, 1989) guarantees that all higher moments are of order o(Δ*t*) as well.

Because of the same reasoning, only the first two of the boundary conditions (3.11) survive:(3.13)where(3.14)(3.15)

The first of the conditions in equation (3.13) is the standard one in the theory of the Fokker–Planck equation and means that the normal component of the probability current vanishes at the boundary. This translates very crudely to the intuitively obvious statement that the probability cannot ‘flow’ beyond the queue boundaries *ℓ*=0 and *ℓ*=*ℓ*_{max}. The second condition is quite peculiar, and to the best of our knowledge, has not been considered before.1 It is difficult to attribute the behaviour of the current to either the second-order moments or the PDF *p*(**x**, *t*) in the regions close to the borders, as we do not know the behaviour of the moments in this region and cannot easily quantify their dependence on the detailed packet dropping procedure.

This is a symptom of the failure of the Fokker–Planck equation to account for the behaviour of the dynamical system of queues in the small region around the border.

The preoccupation we have with the detailed boundary conditions for the PDF arises owing to the fact that some very important quantities of interest (e.g. the packet loss rate of the network) are almost entirely dependent on the behaviour of the PDF in the border region (not just *at* the boundary).

An alternative theoretical approach is to construct an integral equation describing the dynamical behaviour of the transitional PDF in this border region. As can be seen in the outline derivation given in Appendix A, this approach has the advantage that we can construct the boundary region behaviour of the PDF, by incorporating detailed information on the packet dropping procedure *ab initio*.

In this paper, we concentrate on calculating the PDF in the region away from the boundaries, which is suitable for quantifying observables such as the average end-to-end delay. To solve equation (3.10) or (3.12), in addition to the boundary conditions, we need the explicit form of the coefficients . Two examples of the derivation of these moments are presented in Appendices B and C, but in the main body of the work we only restrict ourselves to the first two moments alone.

### (a) Stationary case

In this section, we consider the steady-state (stationary) solution of equation (3.10). The steady-state solution is important for studying static properties of the system, when the external parameters (load) are unchanged. For example, the position of the congestion threshold as quantified by the average end-to-end delay falls in this category. Even when the load is changing in time, but slowly, so that the system quickly comes back to its equilibrium (stationary) state,2 we can use the steady-state solution with the external parameters dependent on time (this is the so-called adiabatic approximation). Many results presented in this section are well-known, but we include them for completeness.

If the process is homogeneous (the coefficients are time independent), then equation (3.10) takes the following formwhere *p*_{s}(**x**) is the stationary distribution. This equation can be simply written in terms of the probability current of equation (3.10):(3.16)

In the one-dimensional case (**x**=*x*), which describes the dynamics of an isolated queue, equation (3.16) takes the form(3.17)

Applying the boundary conditions (3.13) we find that the constant in equation (3.17) is equal to zero:

It should be kept in mind that terms with *k*≥3 are treated as a perturbation. Considering only terms *k*=1,2 corresponds to the stationary regime of the ordinary Fokker–Planck equation, and it is well known that the solution in this case can be represented in quadratures (e.g. Gardiner 1994).

Turning our attention to the multi-dimensional case, we note that one possibility to satisfy equation (3.16) is to assume that

In what follows, we restrict ourselves by considering a purely Markovian case (terms *k*=1, 2 only). Hence, from the equation above we obtain(3.18)

Provided that the matrix is invertible for all **x**, we have from equation (3.18)(3.19)

This equation cannot be satisfied for arbitrary and since the LHS is explicitly a gradient. Hence, *F*_{i} must also be a gradient, and a necessary and sufficient condition for this is(3.20)which is called the potential condition. If this is satisfied, the stationary solution can be obtained by simple integration of equation (3.19):

The value of the integral does not depend on a concrete choice of the contour of integration because *F*_{i} is defined as the gradient of a potential. The potential condition (3.20) coincides with the condition of the detailed balance. This means that condition (3.20) should be satisfied for the system to reach a true equilibrium state (not just a stationary one), which we will assume is the case in what follows.

The detailed balance condition can be interpreted as follows. In the absence of the routing protocol, the coefficients of the Fokker–Planck equation are not dependent on the queue lengths, and hence, the detailed balance condition is always satisfied. Thus, this condition is actually a restriction on the properties of the routing protocol. The protocols which satisfy the detailed balance condition can be considered as *stable* ones, as they do not generate vortices of the probability current inside a subsystem of several interconnected buffers in the stationary regime.

If the potential condition is not fulfilled, one should use a more complex method to derive the stationary solution (see Stratonovich 1989).

### (b) Linear Fokker–Planck equations

In this section, we consider a special type of Fokker–Planck equations: the linear Fokker–Planck equations (these can actually be related to the Ornstein–Uhlenbeck process; Gardiner 1994; Carmichael 1999; Risken 1989). For this type of system, the matrix is constant, and the coefficients are linear in *x*_{i} (e.g. buffer occupancies):(3.21)

This form of the coefficient accounts for the simplest feedback in the system: the inflow (or outflow) in a queue linearly depends on the lengths of other queues. The importance of equation (3.21) arises from the fact that this is simplest way to account for the presence of a routing protocol sensitive to congestion: the inflow to a queue depends on the state of all other queues in the subsystem (including itself) and the linearity is the simplest dependence (after none).3 This is in essence an idealized description of the manner in which the original routing protocols on the ARPANET were intended to function (,McQuillan *et al*. 1980). The coefficients *b*_{ij} can be considered as the corresponding protocol sensitivities (namely, *b*_{ij} is the sensitivity of the inflow of *i*th queue to the state, or the filling factor, of the *j*th queue). The potential condition of equation (3.20) then takes the following form(3.22)that is, the matrix *σ.b* should be symmetric. The detailed balance equation is(3.23)

We seek a solution of the form(3.24)where *A*_{i}, *B*_{ij} are quantities to be determined and *C* is the normalization constant. Inserting equation (3.24) into equation (3.23) and equating coefficients of powers of *x*_{i}, we obtain

Protocols satisfying equation (3.22) belong to the family of stable protocols.

## 4. Examples of subsystems

In this section, we consider three subsystems of the type of equation (3.21): two cascaded queues, *n* output queues connected at the output of the same switching device (this is a subsystem which is part of a router), and a system of *n* twin cascaded queues connected to the output of the same switching device. The first two have been considered in work by Stepanenko *et al*. (2002), and we only summarize them here. The third one is especially interesting because it can be used as the building block in implementing a non-trivial mean-field approximation, that is, to break a network consistently into a set of subsystems of the same type (see §3), thus enabling us to calculate an approximation of the partition function of an entire network. This latter calculation lies outside the scope of the present paper and will be presented separately.

### (a) Two cascaded queues

As an example of a subsystem, we first consider two cascaded queues consisting of an output buffer of a router connected to an input buffer of an adjacent router (see figure 1). The Fokker–Planck equation for this subsystem is of the type of equation (3.21) and has the following form:where

All quantities are normalized to the buffer length (which for simplicity is assumed the same for both buffers), and the boundary conditions of equation (3.13) are imposed,where

Parameters *a*′, *σ*′^{2} characterize the mean value and the variance per unit time of the traffic coming into the first buffer; *a*_{line}, *σ*_{line}^{2} characterize the mean value and the variance of the traffic on the line between the two buffers; and *a*″, *σ*″^{2} characterize the mean value and the variance of the switching capacity available to the second buffer (all quantities are per unit time). We do not discuss here the nature of these parameters, as they should be defined self-consistently with other parameters of a broader model for the overall network. Parameters *b*_{1}, *b*_{2} are sensitivities of the routing protocol to the congestion of the queues. The dependence of the Fokker–Planck equation on the protocol sensitivities can be explained as follows. Relative queue lengths *ℓ*_{1}, *ℓ*_{2} quantify the congestion level of the subsystem. The factor *a*_{1}−*b*_{1}*ℓ*_{1}−*b*_{2}*ℓ*_{2} determines the average amount of traffic diverted to the subsystem: the more it is congested the less traffic (on average) is (or should be) diverted to it.

We seek an equilibrium solution of the Fokker–Planck equation for which the detailed balance condition must be satisfied (Risken 1989), and this condition demands that the parameters *b*_{1}, *b*_{2} are constrained as follows:(4.1)with *b* being a single free parameter characterizing the sensitivity of the stable protocol to the congestion of the subsystem as a whole. The stationary solution iswhere

The normalization constant is a lengthy linear combination of error functions of different arguments, which is omitted here for the sake of brevity.

### (b) System of *n* queues connected to a switching device

Now, we analyse a system of *n* interacting queues connected to the same switching device as shown in figure 2. The Fokker–Planck equation for the PDF of this subsystem, *p*=*p*(ℓ_{1}, …, *ℓ*_{n}; *t*), is(4.2)where*b*_{ik} is the set of protocol sensitivities (similar to the ones described earlier), *a*_{serv}, *σ*_{serv}^{2} are the mean value and variance (per unit time) of the overall traffic coming from the central switching device, is the set of mean values of traffic coming towards of individual queues, *σ*′_{ik} is the corresponding covariance matrix, and *a*″_{i}, *σ*″_{i}^{2} are the mean value and variance of switching capacities available at the outputs of the individual queues. The equilibrium solution has the following form:(4.3)where is the normalization constant, and the following set of relations (owing to a detailed balance condition) should be imposed on the protocol sensitivities *b*_{ik}:(4.4)

A general solution to equation (4.4) allows a large number of free parameters. In order to arrive at a reasonable number of those, we impose the following restrictions:

The first condition means that the mean value of traffic arriving at an individual queue is proportional to the service rate capacity available to it, the second condition denotes that the variance of the traffic arriving at all the queues is the same, and the third is the statement that the correlation coefficients of the incoming traffic on any pair of queues are equal as well, in order to maintain symmetry (this implicitly assumes that the system is homogeneous). In addition to this, we are looking for a solution in the following class (the sensitivities of the idealized routing protocol for one particular queue, to congestion on all other queues, on the same switch are the same):

Then, we have the following parameterization for *σ*_{ik}, *b*_{ik}:where

Summarizing the resulting free parameters, we have: *a*″_{i}, *σ*″_{i}^{2} are the mean value and variance per unit time of the switching capacities available to each queue; *a*_{serv}, *σ*_{serv}^{2} are the mean value and variance per unit time of the overall traffic generated by the central switching fabric; *σ*′^{2} (which is variance of the traffic arriving at each queue) and are the ensemble average characteristics of the routing protocol.

### (c) *n*-port router

Finally, we consider a model of part of an *n*-port router: *n* twin cascaded queues connected to the same switching device; a combination of two previously considered subsystems. The Fokker–Planck equation for the PDF of this subsystem, *p*=*p*(ℓ_{11}, …, *ℓ*_{1n}, *ℓ*_{21}, …, ℓ_{2n}; *t*), is(4.5)where*b*_{ik}^{(1)}, *b*_{ik}^{(2)} is the set of protocol sensitivities similar to the ones described earlier, *a*_{serv}, *σ*_{serv}^{2} are again the mean value and variance (per unit time) of the overall traffic coming from the central switching device; is the set of mean values of the traffic switched towards individual output queues; *σ*′_{ik} is the corresponding covariance matrix; *a*_{line,i}, *σ*_{line,i}^{2} are the mean value and variance of the traffic flowing between the corresponding output and input buffers; and *a*″_{i}, *σ*″_{i}^{2} are the mean value and variance of the switching capacities available to the input queues of adjacent nodes.

The equilibrium solution has the form of equation (3.24) and is determined by the following set of equations (see (3.23)):(4.6)(4.7)

To satisfy the potential conditions of equation (3.22), we first impose restrictions such as equation (4.1) on the elements of the sensitivity matrices of the output and input queueswhere *b*_{ik} is a new set of protocol sensitivities (common for both output and input queues). Introducing a new set of variableswe obtain from the second equation in (4.7)Hence, the solution of (4.5) can be represented as follows:(4.8)

Substituting (4.8) into the first equation (4.6) we obtain(4.9)where differs from *J*_{1i} by the exponential factor in (4.8) and

From equation (4.9), it follows that the solution for is determined by the Fokker–Planck equation of the type of equation (4.2) with the solution similar to equation (4.3). The potential condition of equation (4.4) can be satisfied in the same way as in §4*b* with the appropriate symbol substitutions.

In this section, we have considered several examples of subsystems. As we have seen, the solution for the PDF of a subsystem may contain a large number of free parameters. However, the potential condition can be employed to decrease this number significantly. To further restrict the number of parameters in this analysis, we employed further simplifying assumptions on the operation of the protocol and the environment of the subsystem (those were mostly symmetry considerations). To arrive at a more realistic description of these subsystems, it is necessary to resort to experimental observations (which at the moment are hardly available), as they should involve traffic measurements on several ports of a router simultaneously.

## 5. Conclusion

In this paper, we present an outline phenomenological teletraffic theory of large-scale, core, packet-switched networks. All approximations made in the construction of this approach have preserved the description of the sources of stochasticity throughout, and maintain the discrete nature of packets without resorting to fluid-flow approximations. The careful selection of time-scales in the description of the underlying stochastic processes enables us to employ the Markovian approximation to the latter.

Furthermore, all approximations made take cognisance of the fact that the teletraffic theory we aim to construct only needs to have predictive powers in the regime of heavily loaded to congested network operation.

In order to make the analysis tractable, we identify strongly correlated subsystems of queues in the network which do not necessarily coincide with the distinct physical devices such as routers, and so on. We then proceed to build a probabilistic model which describes the ensemble average operation of each such subsystem in the presence of a routing protocol.

Subsystems are described through Fokker–Planck type differential equations governing the joint probability density function of the queue occupancies in the subsystem. Careful investigation of the boundary conditions pertinent to the problem at hand reveals that the behaviour of this PDF in a small region close to the boundary is determined by the packet dropping procedure and the manner in which this influences the resulting traffic statistics. Two levels of approximation are possible, one of which is approximately valid outside this border region and is thus not applicable to the calculation of quantities such as the packet loss rate, which by definition is significant only when buffers are nearly full. This form is described accurately by the Fokker–Planck equation with natural boundary conditions. The more exact form is described by an integral equation, which is only outlined in an appendix and forms the subject of a forthcoming paper.

A basic building block model of a switching fabric connected to a cascade of output queues, transmission lines and input queues of neighbouring network nodes is presented. This part of the analysis highlights the need for detailed correlated multipoint measurements in the core part of the Internet to support the full development of such a phenomenological theory.

This work is far from complete, but represents a first step towards the construction of a teletraffic theory for core, packet-switched networks. Work is in progress to integrate the partial router/transmission line building blocks to construct a simple overall network model.

## Acknowledgments

The authors wish to thank Nortel Networks for the initial funding of this work. In particular, we wish to thank Kevin Baughan and Colin Camilleri for their technological input. We also wish to thankfully acknowledge the current support of this project by an EPSRC ROPA award.

## Behaviour of PDF in the border region and the loss rate

In this appendix, we present an integral equation for the generating functional of conditional moments of the dropped traffic on time-scale *t* (under the condition that the queue is in state *ℓ* at the beginning of the observation interval). Details of this derivation will be published elsewhere. In this appendix, we will measure time in units equal to the average time of a cycle (cycle=packet+the following gap), that is, , , here is average length of a cycle, are the average packet and gap lengths and *c*_{in} is the capacity of the input line ( *c*_{in} are normalized to the buffer size *L*).

This functional can be expressed as follows (*ξ* is the Fourier conjugate variable that generates the moments of the dropped traffic):where *w*(ℓ′, *t*|ℓ; *ξ*) is the generating functional of the moments of the dropped traffic under the condition that the queue is in state *ℓ* at the beginning of the observation interval *t* and in state *ℓ*′ at the end of it. If we put *ξ*=0 then *w*(ℓ′, *t*|ℓ; 0) is simply the transitional probability density from {*ℓ*, 0} to {*ℓ*′, *t*}. The queue length *ℓ* is again normalized to the buffer size *L* (*ℓ* is a filling factor).

Here we adopt the following packet dropping procedure: if on arrival of a packet there is not enough space to accommodate it, the packet is dropped altogether. Considering the time-scale *t* covering, on average, a large number of packet arrivals, we can obtain the following expression for *w*(*ℓ*′, *t*|*ℓ;ξ*):where *δ*(*x*) is the Dirac *δ*-function, the angle brackets denote averaging over all packet (*p*_{k}) and gap (*g*_{k}) lengths and *w*_{n}(*ℓ*′|*ℓ*; *ξ*) is the conditional generating functional (to go from *ℓ* to *ℓ*′ in *n* cycles) of moments of the dropped traffic (*w*_{n}(*ℓ*′|*ℓ*; 0) is the transitional probability density to go from *ℓ* to *ℓ*′ in *n* cycles), for which we have the following recurrence relation:(A1)with the initial conditionwhere *w*_{s}(*ℓ*′, Δ*t*|ℓ) is the transitional probability density from *ℓ* to *ℓ*′ during time Δ*t* owing to service and ,∀*n*.

The term proportional to *θ*[−1+*ℓ*_{n}+*p*_{n+1}] in equation (A 1) corresponds to an event of packet dropping. Hence, we can express *w*(*ℓ*′, *t*|*ℓ*; *ξ*) as an expansion in the number of packet dropping occurrences:

For the zeroth-order term in this expansion *w*^{(0)}(*ℓ*′, *t*|*ℓ*; *ξ*)≡*w*^{(0)}(*ℓ*′, *t*|*ℓ*) we have(A2)

It can be shown that *w*(*ℓ*′, *t*|*ℓ*; *ξ*) satisfies the following integral equation:(A3)

The formal solution of equation (A 3) iswhere *ℓ*_{m}≡*ℓ*′.

In the same way, it can be shown that *w*^{(0)} (*ℓ*′, *t*|*ℓ*) defined by equation (A 2) satisfies the following equation:

The quantity *g*_{0}(*ℓ*′, *t*|*ℓ*) is determined by the solution of the following boundary value problem:(A4)with straightforwardly determined 4. The solution of equation ,(A 4) is

It can be shown that *w*^{(0)}(*ℓ*′, *t*|*ℓ*) can be represented as an expansion in powers of (which is a small parameter around the congestion threshold), and the zeroth-order term in this expansion is (assuming that (A5)

For *ℓ*′, *ℓ*<1 expression (A 5) coincides with the solution of the Fokker–Planck equation with the homogeneous boundary condition at *ℓ*′=1, and this is, what we would actually expect. If *ℓ*′>1, *ℓ*<1 or *ℓ*′<1, *ℓ*>1, we have *w*^{(0)} (*ℓ*′, *t*|ℓ)=0, which is again consistent with our initial assumptions.

## Model of packet flow in a transmission line

In this appendix, we consider the flow of packets in a line of capacity *C*. We will assume that the observation interval is large, that is, covers on average, a large number of packets. We are looking for a conditional PDF *w*_{line}(*ℓ*|*t*) for the amount of information *ℓ* (number of bits) in transit during the time-interval *t*.

The conditional PDF is defined by the following relation:(B1)where *p*_{line}(*ℓ*, *t*) is the joint PDF, and we impose a *Ct* upper integration limit in the second equality (B 1), since the amount of information in transit cannot exceed the amount permitted by the line capacity. The joint p.d.f's *p*_{line}(*ℓ*, *t*) can be represented as follows:(B2)where *p*_{k}(*g*_{k}) are consecutive packet (gap) lengths, the angled brackets denote averaging over all possible packet and gap lengths involved. The construction of equation (B 2) is as follows: the first *δ*-function selects configurations with an overall observation interval equal to *t*, and the second one ensures that the overall traffic is equal to *ℓ*. In principle, we should consider four cases, namely, that the observation interval starts inside a packet and ends inside a packet and the other three cases being ‘packet–gap’, ‘gap–packet’ and ‘gap–gap’. But it can be shown, that, for a large enough interval covering on average a large number of packets, they all give the same result up to multiplicative constants which are cancelled in equation (B 1). Hence, we consider only one these cases. Moreover, for the same reason we can assume that the observation interval starts at the beginning of a packet and ends at the end of a gap (this only simplifies formulae but does not have any impact on the final result).

Introducing the characteristic function of *w*_{line}(*ℓ*|*t*) (and recalling that it can be represented as a cumulant expansion)(B3)we obtain from equation (B 2)(B4)where are the characteristic functions for packet (*p*_{p}(*p*)) and gap (*p*_{g}(*g*)) length distributions, respectively.

Calculating the integral in equation (B 4), we obtain for the first two cumulants in expansion (B 3)where , *σ*_{p}^{2} (*σ*_{g}^{2}) are mean value and variance of packet (gap) lengths.

## Switching process

In this appendix, we consider the switching process in the router. We assume that there are *n* input queues on the router. Each packet is divided into cells of size *s*_{cell}. The last cell of a packet is used only partially (the remainder of the packet is probably less than the size of the cell), and this is the source of stochasticity in the switching process. All queues are served in turn—one cell a time in a round robin fashion. So, if *C*_{serv} is the capacity of the switch (in cells per second), then the rate of servicing of each individual queue is . This implicitly assumes that all queues are non-empty. As we are interested in studying a heavily loaded network, then this is a plausible assumption.

As was noted in Appendix B, if the observation time-interval is large (covering a large number of packets), we can choose the observation interval starting at the beginning of the first cell of a packet and finishing at the end of the last cell of another packet. We measure the observation interval by the number of processed cells *N*=*C*_{serv}*t*/*n*. As in the previous appendix, we seek the conditional PDF *w*_{serv}(*ℓ*|*t*) of number of bits *ℓ* served over the time-interval *t*. In complete analogy with equation (B 1), the conditional PDF is defined via the corresponding joint PDF *p*_{serv}(*ℓ*, *t*):

In complete analogy with equation (B 2), we have(C1)where the angled brackets denote averaging over all packet lengths and associated numbers of cells, for example, for each packet the averaging is of the formwhere *p*_{cell}(*p*, *m*) is a joint PDF for the length of packet *p* and the number of cells *m* it occupies, and *p*_{p}(*p*) is the packet length PDF.

Using equation (C 1) and introducing the following characteristic function for *p*_{cell}(*p*, *m*),we have for the characteristic function of *p*_{serv}(*ℓ*, *t*)(C2)

Calculating the integral in equation (C 2) we obtainwhere , are the average packet length and its variance, , are the average number of cells in a packet and its variance, is the average amount of data in a cell, and *c*^{(1,1)} is the covariance between the packet length and the number of cells it occupies.

The results of this appendix along with the results of Appendix B can be used to construct the coefficients *m*^{(k)} in the Fokker–Planck equation.

## Footnotes

As this paper exceeds the maximum length normally permitted, the authors have agreed to contribute to production costs

↵Normally the problem of boundary conditions such as equation (3.15) is avoided by considering observables

*f*(**x**) that tend to zero, with all their derivatives at the boundary (e.g. Gardiner 1994).↵This actually means that the relaxation time of the system is much smaller than the characteristic time on which the load changes significantly.

↵An ‘engineering’ approximation to the more useful case of nonlinear reaction of routing protocols to congestion can be dealt with by approximating this nonlinear dependence in a piecewise linear manner. The implementation of this approach is laborious, but straight-forward.

↵Parameters , in equation (A 4) are defined per unit time

*τ*and differ from their normally defined counterparts*a*,*σ*^{2}(per second) by an extra factor*τ*.- Received January 14, 2004.
- Accepted September 14, 2004.

- © 2005 The Royal Society