<figure><img src="image_1.png" /></figure>

HAL
open science

# Modeling of packet switching communication networks with finite buffer size at each node

Jacques Labetoulle, Guy Pujolle

► To cite this version:

Jacques Labetoulle, Guy Pujolle. Modeling of packet switching communication networks with finite buffer size at each node. [Research Report] IRIA-RR-239, IRIA. 1977, pp.22. hal-04716467

HAL Id: hal-04716467

https://inria.hal.science/hal-04716467v1

Submitted on 1 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
IRIA

laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B. P. 105 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

ARCHIVES

<figure><img src="image_3.png" /></figure>

MODELING OF
PACKET SWITCHING
COMMUNICATION NETWORKS
WITH
FINITE BUFFER SIZE
AT EACH NODE

Jacques LABETOULLE
Guy PUJOLLE

Rapport de Recherche N° 239

Juin 1977
MODELING OF PACKET SWITCHING COMMUNICATION NETWORKS
WITH FINITE BUFFER SIZE AT EACH NODE

Jacques Labetoulle, Guy Pujolle
IRIA/LABORIA

<figure><img src="image_4.png" /></figure>

Résumé :

Le sujet de ce travail est l'étude des réseaux à commutation de paquets, quand on considère des capacités de mémorisation limitées à chaque nœud. A notre connaissance il n'existe aucune méthode analytique exacte permettant de donner le temps de réponse de tels réseaux. Aussi nous proposons ici une solution numérique rapide qui donne les probabilités de saturation à chaque nœud et le temps de réponse moyen du réseau. Elle permet également de savoir si le réseau peut ou non accepter une charge donnée.

Cette méthode est appliquée au réseau CYCLADES et les résultats sont comparés à des simulations.

Abstract :

In this paper we study packet switching networks with finite storage capacity at each node. So far we do not know any method giving an analytical result to the problem of the response time of such networks. We propose here a very fast, thus very economical, numerical solution which gives the probabilities of saturation of each node, the mean response times of the network and permits to know if the network can accept a given load without saturation. This method is applied to the CYCLADES network and the results are compared with simulations.
# INTRODUCTION

All store-and-forward packet switching networks built until now have a finite storage capacity at each node. Until now, this choice has not been justified by convincing theoretical studies and has been frequently justified by considerations of simplicity or by experience.

In this paper we study a model of store-and-forward packet switching networks in which nodes will be of limited capacity. An approximate solution will be developed to compute the expected response times and the stability conditions of the model.

Today, two types of service are proposed for the switching of packets : the virtual circuit and the datagram service. The virtual circuit implies call procedures for the setting of a circuit and deliverance procedures (ARPA or EPSS). In this case we have an end traffic control, and our study will determine stability conditions. Expected response times will be only very approximate due to the procedures of call and deliverance.

The datagram service is a simple packet transmission service. There are neither procedures of call or of deliverance, nor end-to-end traffic control. Each node of the network carries out a control on the receipt of the packets. When there is an error or a non-receipt, the preceding node retransmits. Here, a node -to-node protocol intervenes, which we assume very simple in our study : a negative acknowledgment which returns immediately at the preceding node. More sophisticated protocols will be studied later.

It is the behavior of such a network that we study in the sequel (CYCLADE/CIGALE, EIN,...).

Till now, such systems have been modeled by networks of infinite capacity queues, assuming that the overflow probabilities are negligible. The queueing delay at each node with finite buffers size at low probability of overflow ($< 10^{-4}$) could
be approximated as the queueing delay at infinite buffer size as shown in [20].

To solve these models three types of methods exist :

1 - The classical Baskett, Chandy, Muntz and Palacios solution [1].

2 - Diffusion methods developed by Kobayashi [2] and Gelenbe [3], [4]. The method consists in representing the state of the system by a non discretized variable, or a vector. The equation which governs this variable is of a type found in classical problems of diffusion (heat for example).

3 - Iterative techniques based on heuristic considerations such as the method of Chandy et al. [5].

Recent studies of limited capacity networks have been made, but mainly for tandem queueing systems: Pennotti and Schwartz [6], Rudin [7] for exponentially distributed service times, Rubin [8], Labetoulle and Pujolle [9] for constant service times cases, Pujolle [10] for general service times.

Other studies compute the overflow probabilities at nodes : Lam [11], Schweitzer and Lam [12], or determine stability conditions : Lavenberg [13], Schweitzer [14].

Here, we develop a numerical solution technique which includes most of these results. Applications to the seven node Cyclade network are given in the last section.

## Description of the model

We consider a model with two levels (Figure 1) :

<figure><img src="image_5.png" /></figure>
- The store-and-forward packet switching itself. A fixed routing is assumed, i.e. all the links used by a packet depend only on the departure and destination nodes. Two nodes are linked by at most one link and each link has its own storage capacity of M packets. Packets are sent in FIFO order. When a packet is transmitted from a node it can :

- leave the network if the arrival node is the destination node

- enter the queue of the next link

- if this last queue is full, then it comes back to the preceding queue.

- The hosts : a host sends packets in the network and receives the packets leaving the network. We assume that the speed of the host is sufficiently high to avoid the rejection of packets when they arrive at their destination node. The packets received by a host to be sent in the network can be rejected if the buffers of the first following link are full. In this case the host keeps in its own buffer which must be sent again later. (See figures 1 and 2).

Therefore we model the host as :

- an infinite speed server when it is a receiver

- a set of infinite capacity queues when it is a sender.

**Remark** : It is possible to assume a network with limited capacity nodes (the total number of packets waiting at the node is limited). This case could be treated as an extension of our model with limited links.

## Method of solution

We assume in this section that interarrival times from the outside and service times (retransmission times) are exponentially distributed.

<figure><img src="image_6.png" /><figcaption>Figure 2</figcaption></figure>

It is obvious that the system precedingly described cannot be solved with an exact analytical method. We propose an approximate method based on the following idea : the arrival process to a queue of the network is the juxtaposition of several
flows : (see figure 2)

- the flow from outside (the host)

- the flows coming from neighbouring nodes

because of limited capacity each flow cannot be govern with exponential distributions. However we shall assume that the superposition of these flows is Poisson. If we assume an infinite number of joined traffic, then the assumption of exponential distribution would be true. Moreover this assumption has been verified to be very good even for tandem queueing networks [15]. That proves the validity of our assumptions for this model.

We can equivalently consider an arriving flow $\lambda = \sum \lambda_i$ with packet recycling in case of saturation or a flow equal to $\lambda/(1-p)$ assuming packets are lost in case of saturation. ($p$ is the probability of saturation).

In the same way, we may consider an equivalent server which is a function of routing probabilities ($q_i$, $i=1,...,K+1$) and of saturation probabilities ($p_j$) of the following queues. This equivalent server is approximated by considering that a packet leaving the station and going towards station j has the stationary probability $p_j$ of being rejected.

Let $\mu$ denote the service rate and let $q_{K+1}$ denote the probability that a packet leaves the network.

This equivalent server has a service time distribution which is still exponential with rate $\tilde{\mu}$ where :

$$ \tilde{\mu}^{-1} = \mu^{-1} \left( q_{K+1} + \sum_{j=1}^{K} \frac{q_j}{1-p_j} \right) \quad (1) $$

This formula becomes obvious if we consider the expected time to serve a customer going to queue j. This quantity is equal to :

$$ \frac{1-p_j}{\mu} + \frac{2p_j(1-p_j)}{\mu} + \dots + \frac{np_j^{n-1}(1-p_j)}{\mu} + \dots = \frac{1}{\mu(1-p_j)} $$

The equivalent station is an M/M/1/N queue and the probability to have a customer rejected is equal to the probability that all the buffers are occupied, therefore :

$$ p = \frac{\tilde{\rho}^M \frac{1-\tilde{\rho}}{1-\tilde{\rho}^{M+1}}}{1-\tilde{\rho}} \quad (2) $$

$$ \text{with } \tilde{\rho} = \frac{\lambda}{\mu(1-p)}. \quad (3) $$

## Notations and system of equations governing the model

We consider a network with K transmission links and for each of them, we define the following notations :

- $A(i)$ represents the set of links arriving to the links i

- $S(i)$ represents the set of links following the link i

- $M_i$ is the maximal number of buffers of link i

- $\lambda_i$ is the traffic rate arriving to link i
- $\lambda = \sum_{i=1}^{K} \lambda_i$ is the total arriving flow

- $\mu_i$ is the retransmission rate of link i (speed of the link)

- $\tilde{\rho}_i$ is the equivalent traffic intensity of link i : $\tilde{\rho}_i = \frac{\rho_i}{1-\rho_i}$

- $p_i$ is the probability of rejection or the probability that the link is full

- $\lambda_i^{\text{ext}}$ is the rate of customers arriving from outside of the network (the host)

- $q_{ij}$ is the branching probability to link i

- $q_{i,K+1}$ is the probability to leave the network after link i

- $\rho_i$ is the traffic intensity of link i : $\rho_i = \lambda_i / \tilde{\mu}_i$

and

$$ \frac{1}{\tilde{\mu}_i} = \frac{1}{\mu_i} \left( q_{i,K+1} + \sum_{j=1}^{K} \frac{q_{ij}}{1-p_j} \right) \quad (4) $$

The model is governed by the following set of equations :

$$ p_i = \tilde{\rho}_i \frac{1 - \tilde{\rho}_i}{M_i + 1} \quad i = 1, \dots, K \quad (5) $$

## Conservation of flow

In this last system of equations we have :

- an equivalent server computed from the blocking probability

- an equivalent flow arriving at queue i, computed from external traffic flows.

This equivalent traffic flow will be the real flow only if the network is stable, in the sense of the ergodicity of queueing networks, i.e. if all the packets leave the network after a finite time.

Now, we shall prove that the steady state flows are the flows we use. Notice that

$$ \lambda_i = \sum_{j=1}^{K} \tilde{\mu}_j (1-p_j(0)) q_{ji} + \lambda_i^{\text{ext}} , \quad 1 \le i \le K $$

where $p_j(0)$ is the stationary probability that server j is idle. We must prove that the $\lambda_i$, $1 \le i \le K$, also satisfy the equations

$$ \lambda_i = \sum_{j=1}^{K} \lambda_j q_{ji} + \lambda_i^{\text{ext}} $$

which express the conservation of flow. It is sufficient to prove that

$$ \tilde{\mu}_j (1-p_j(0)) = \lambda_j $$

where we have
$$p_j(0) = \frac{1 - \tilde{\rho}_j}{1 - \tilde{\rho}_j^M}$$

Therefore

$$\tilde{u}_j(1-p_j(0)) = \tilde{u}_j \tilde{\rho}_j \left( \frac{1-\tilde{\rho}_j^M}{1-\tilde{\rho}_j^{M+1}} \right) \quad (6)$$

and, using (5) in (6) we write

$$\tilde{u}_j(1-p_j(0)) = \tilde{u}_j \tilde{\rho}_j(1-p_j)$$

This completes the proof, since $\tilde{\rho}_j = \lambda_j / \tilde{u}_j(1-p_j)$.

The system will be ergodic if the queues in the hosts remain finite. In the following section, we study conditions under which the system (5) has a solution with $p_i < 1$, $\forall i$. We shall show that these conditions guarantee the ergodicity of the network.

## Existence and uniqueness of the solution

First we must know if a solution to system (5) exists, different from the trivial solution $p_i = 1$, i. For this reason we begin to study the system of equations

(2) and (3) :

$$(2) \quad p = \frac{\tilde{\rho}^M}{1 - \tilde{\rho}^{M+1}}$$

$$(3) \quad \tilde{\rho} = \frac{\lambda}{\mathfrak{H}(1-p)}$$

$$\text{or} \quad 1-p = \frac{1-\tilde{\rho}^M}{1-\tilde{\rho}^{M+1}}$$

$$1-p = \frac{\lambda}{\mathfrak{H} \tilde{\rho}}$$

Therefore a positive solution $\tilde{\rho}$ to the equation

$$\frac{\lambda}{\mathfrak{H} \tilde{\rho}} = \frac{1-\tilde{\rho}^M}{1-\tilde{\rho}^{M+1}}$$

gives a unique solution to the system (2), (3).

So, let us study the function $f(\tilde{\rho}) = \frac{\tilde{\rho} - \tilde{\rho}^{M+1}}{1 - \tilde{\rho}^{M+1}}$. Notice that

$f(0) = 0, f(+\infty) = 1$. An elementary study of the function $f$ shows that

$$f'(\tilde{\rho}) = \frac{1-(M+1)\tilde{\rho}^M + M\tilde{\rho}^{M+1}}{(1-\tilde{\rho}^{M+1})^2}$$

is always positive for $\tilde{\rho} \in [0, 1]$.

Consequently, the system (2), (3) has a unique solution with $0 < p < 1$ if and only if $\rho = \frac{\lambda}{\mathfrak{H}} < 1$.

If $\frac{\lambda}{\mathfrak{H}} \ge 1$, the only possible solution with $p \in [0, 1]$ is $p=1$. In this case, we see
that the quantity $\tilde{\rho}$ becomes infinite.

The system (5) may be represented by a matrix system :

$$ (7) \quad \vec{P} = \phi(\vec{P}) $$

where $\vec{P}$ is the vector with components $p_i$, $i = 1,...,K$ and $\phi$ is a vectorial function describe by (5).

The system (7) has at least one solution which is the vector $\vec{P} = \vec{i}$ or $p_i = 1, \forall i$.

Now, we can prove that if another solution exists for system (7), this solution is unique. Notice that $\phi$ is a monotone non-decreasing function for each component $p_i$ i.e. the reduction of $\phi(\vec{P})$ at each component is a non-decreasing function. This property is obvious from the system (5). We use the classical notation $\vec{P}_1 < \vec{P}_2$ if the inequality is satisfied for all the components.

Let us consider the vectorial series :

$$ \vec{P}^0 = \vec{0} \qquad (p_i^0 = 0, \forall i) $$

$$ \vec{P}^{i+1} = \phi(\vec{P}^i) \qquad (8) $$

The vector $\vec{0}$ not being solution of (7), we have $\vec{P}^1 > \vec{0}$. We deduce :

$$ \vec{P}^2 > \vec{P}^1 \text{ because } \vec{P}^1 > \vec{P}^0 \text{ and } \phi(\vec{P}^1) \ge \phi(\vec{P}^0). $$

and by induction :

$$ \vec{P}^{i+1} > \vec{P}^i $$

Each component $p_j^i$ of $\vec{P}^i$, $i = 1,...,\infty$ represents a non-decreasing series, bounded from above by 1.

Therefore this series has a limit which satisfies $p_j = \phi_j \vec{P}^\infty$ and $\vec{P}^i$ has a limit $\vec{P}$ which satisfies $\vec{P} = \phi(\vec{P})$.

We can see that, $\forall \vec{P}'$ solution of (7), $\vec{P}'$ satisfies $\vec{P}' \le \vec{P}'$. In fact :

$\vec{P}^0 < \vec{P}' \rightarrow \vec{P}^1 = \phi(\vec{P}^0) < \vec{P}' = \phi(\vec{P}') \rightarrow \vec{P}^j < \vec{P}' \rightarrow \vec{P} \le \vec{P}'$. So, if (7) has at least one solution, $\vec{P}$ may be called the smallest solution of (7). Now, we shall show that if $\vec{P} < \vec{i}$, then any other solution to (7) has at least one component equal to 1. If $\vec{P} < \vec{i}$, then $p'_i < 1$. The value $p_i$ is solution of equation (5) with $i = 1$, the other $p'_i$ is being the components of $\vec{P}$.

Therefore $p_i$ is a solution of (9) :

$$ p_i = \left( \frac{\rho_i}{1-p_i} \right)^M \frac{1 - \frac{\rho_i}{1-p_i}}{1 - \left( \frac{\rho_i}{1-p_i} \right)^{M+1}} = f_i(p_i) \qquad (9) $$

where $\rho_i$ is a constant.

If $\vec{P}_2$ is a solution of (7) with a first component different from $p_i$ then :
$$ \vec{P}_2 \ge \vec{Q}^0 = \{p_1 + \epsilon, p_2, \dots, p_K\} \qquad (10) $$

$\vec{Q}^0$ is the same as the vector $\vec{P}$ except for the first component, which is increased by an $\epsilon$ sufficiently small to satisfy (10).

Let us consider the numerical series $\{p_1^i\}$ defined by : $p_1^0 = p_1 + \epsilon, p_1^{i+1} = f_1(p_1^i)$.

This series defines a vectorial series $\{Q^i\}$ determined by : $Q^i = \{p_1^i, p_2, \dots, p_K^i\}$.

The series $\{Q^i\}$ is an increasing series and it is obvious by induction that : $\vec{Q}^i \le \vec{P}_2$. But equation (9) have only as solution $p_1$ and 1. Therefore, the first component of $\vec{P}_2$ equals 1.

We deduce that the system (7) possesses at most one solution $\vec{P}$ satisfying $\vec{P} < 1$.

## Numerical computation

Equation (7) can be solved only by a numerical method. In fact, we use the method described in the preceding paragraph which consists in computing the series of vectors $\vec{P}^i$. We stop the iteration when : $||\vec{P}_{i+1} - \vec{P}_i|| < \epsilon$. If at any step i, and for the queue j, we observe that $\rho_j = \lambda_j/\tilde{\mu}_j \ge 1$, then the only solution is $p_j = 1$. In this case we leave the domain of validity of the model.

When the saturation probabilities $p_i$ are known, we may compute the expected number of customers in the queue,

$$ N_i = \sum_{k=1}^{M_i} k \cdot \frac{\rho_i \cdot k \cdot \frac{1 - \frac{\rho_i}{1-\rho_i}}{1 - \left(\frac{\rho_i}{1-\rho_i}\right)}}{M_i + 1} $$

and also the expected response time $W_i : W_i = \frac{N_i}{\lambda_i}$. Hence it is possible to compute the expected response time of any host-to-host path.

## Study of the model outside its domain of validity

If the system is unstable, i.e. if the system (7) has no solution satisfying $\vec{P} < \vec{1}$, some queues in the hosts become infinite. It will be interesting in this case to know the behaviour of the network itself. The flow from a station to the following one cannot exceed the transmission rate of the links. Also, the fact that all the packets cannot be served modifies the probability $q_{ij}$ to go to the queue j from i.

The throughput of a station is equal to its service rate multiplied by the probability that the server is busy. Therefore the throughput of station i is exactly $\mu_i(1-\mu_i(0))$. Let $\tilde{\lambda}_i$ denote the real arrival rate to station i, then

$$ \tilde{\lambda}_i = \sum_{j \in A(i)} \mu_j q_{ji}^* (1 - \mu_j(0)) \qquad (12) $$

$$ \frac{1}{\mu_i} = \frac{1}{\mu_i} \sum_{k \in S(i)} \frac{q_{ik}}{1-p_k} \qquad (13) $$
$$ \tilde{\rho}_i = \frac{\tilde{\lambda}_i}{\tilde{\mu}_i} \qquad (14) $$

$$ p_i = \tilde{\rho}_i^{M_i} \frac{1 - \tilde{\rho}_i}{1 - \tilde{\rho}_i^{M_i}} \qquad (15) $$

$$ p_i(0) = \frac{p_i}{\tilde{\rho}_i^{M_i}} \qquad (16) $$

$$ q_{ji}^* = \frac{\sum_{k \in S(A(i))} (1-p_k)}{1-p_i} \qquad (17) $$

The iterative method of the previous section is used again, in which $p_i(0)$ and the $q_{ij}^*$'s intervene.

## Extension to state-dependent service rates

In the case where the overhead is added to the transmission time, the service time can depend on the state. When the network is congested, the overhead is higher than when traffic is light. This is due, for instance, to the waiting time for acknowledgements.

Let $\mu_i(n)$ be the service rate when there are n customers in the queue. Let us take:

$$ \frac{1}{\tilde{\mu}_i(n)} = \frac{1}{\mu_i(n)} \sum_{k \in S(i)} \frac{q_{ik}}{1-p_k} $$

$$ \rho_i(n) = \frac{\lambda_i}{\tilde{\mu}_i(n)} \quad \text{and} \quad \tilde{\rho}_i(n) = \frac{\lambda_i}{\tilde{\mu}_i(n)(1-p_i)} $$

$$ \text{We obtain : } p_i = \frac{\sum_{n=1}^{M_i} \tilde{\rho}_i(n)}{1 + \tilde{\rho}_i(1) + \tilde{\rho}_i(2) + \dots + \sum_{j=1}^{M_i} \tilde{\rho}_i(j)} \qquad (18) $$

An iterative method may be used as in the previous section.

## Extension of the model to diffusion processes

Up to this point, we have assumed that interarrival times in the network and service times were exponentially distributed. If we make no assumptions on these distributions, we do not know how to express in a simple way the probability of saturation of a finite capacity queue.

Approximations to diffusion processes allow us to obtain that quantity in a simple manner. Let us consider a queue and let $\lambda$, $\mu$, $K_a$, $K_s$ denote arrival rate, service rate, squared coefficient of variation of interarrival times and service times, respectively. With diffusion approximations, the probability of saturation is given by [3, 10].
$$p = \frac{\rho(1-\rho)}{\hat{\rho}^{-(M-1)} - \rho^2} \quad \text{with} \quad \rho = \frac{\lambda}{\mu} \quad \text{and} \quad \hat{\rho} = e^{\frac{2(\rho-1)}{\rho Ka+Ks}} \quad (19)$$

As in the preceding model, we introduce the notions of equivalent flow and equivalent service. The limits of validity of the model will be as before and from now, we shall assume that these conditions are satisfied. The following notations will be used :

$\tilde{K}_{i}^{\text{ext}}$ : squared coefficient of variation of arrival into queue i from the host if there is no feed-back

$K_{i}^{\text{ext}}$ : squared coefficient of variation of arrival into queue i from the host in the case of feed-back

$K_{i}$ : squared coefficient of variation of the equivalent flow arriving into queue i

$Ks_{i}$ : squared coefficient of variation of service time at queue i

$\tilde{K}s_{i}$ : squared coefficient of variation of the equivalent server.

Kobayashi [2] has established a formula giving the squared coefficient of variation of the arrival process in any station of a network. We shall use that formula in the sequel, assuming that the flow leaving a station is not regulated, but is lost. This with the condition to multiply the flow $1/(1-p_i)$. In that case, the formula becomes :

$$K_{i}^{\text{ext}} = \frac{1}{\lambda_{i}} \left\{ \tilde{K}_{i}^{\text{ext}} \lambda_{i}^{\text{ext}} + \sum_{j \in A(i)} [(\tilde{K}s_{j}-1)q_{ji} + 1] \lambda_{j}q_{ji} \right\}$$

$$\text{or} \quad K_{i}^{\text{ext}} = \tilde{K}_{i}^{\text{ext}} + \frac{1}{\lambda_{i}} \sum_{j \in A(i)} [(\tilde{K}s_{j}-1)q_{ji} + 1 - \tilde{K}_{i}^{\text{ext}}] \lambda_{j}q_{ji} \quad (20)$$

$\tilde{K}_{i}^{\text{ext}}$ can be expressed as a function of $K_{i}^{\text{ext}}$ by :

$$\tilde{K}_{i}^{\text{ext}} = P_{i} [K_{i}^{\text{ext}} - 1] + 1 \quad (21)$$

The squared coefficient of variation of service time at station i of a customer going to station j is : [10]

$$K_{j}^{i} = p_{j} + Ks_{i}(1-p_{j}) \quad \text{and its rate is :}$$

$$v_{j}^{i} = \mu_{i}(1-p_{j})$$

The variance of the equivalent service time is :

$$\text{var}_{i} = \sum_{j \in S(i)} \frac{q_{ij} [K_{j}^{i} + 1]}{(v_{j}^{i})^{2}} - \frac{1}{\mu_{i}}$$

Then we obtain the squared coefficient of variation of the equivalent service time at the station i as follows :
$$ \tilde{K}_{s_i} = -1 + \hat{\mu}_i^2 \sum_{j \in S(i)} \frac{q_{ij} [K_j^i + 1]}{\mu_i^2 (1-p_j)^2} $$

and finally substituting $K_j^i$ by its value yields :

$$ \tilde{K}_{s_i} = -1 + \frac{K_{s_i}}{\sum_{j \in S(i)} \frac{q_{ij}}{1-p_j}} + \frac{\sum_{j \in S(i)} \frac{q_{ij}(1+p_j)}{(1-p_j)^2}}{\left[ \sum_{j \in S(i)} \frac{q_{ij}}{1-p_j} \right]^2} \quad (22) $$

Let us denote that :

$$ \tilde{\rho}_i = \frac{\lambda_i}{(1-p_i)\tilde{\mu}_i} \quad (23) $$

$$ \hat{\rho}_i = e^{-2} \frac{1 - \tilde{\rho}_i}{\tilde{\rho}_i K a_i + \tilde{K} s_i} \quad (24) $$

and the feed-back probability can be written :

$$ P_i = \frac{\tilde{\rho}_i (1 - \tilde{\rho}_i)}{\hat{\rho}_i - (M_i - 1) - \tilde{\rho}_i^2} \quad (25) $$

## Numerical computation

The system consisting of equations (20) to (25) for $i = 1,...,K$ can be written in matrix form as :

$$ \vec{P} = \emptyset(\vec{P}), $$

where $P$ is the vector $(p_1,...,p_K)$. The existence and unicity of the solution of this system may be proved as in the exponential case and are treated in the appendix 1.

$\vec{P} = \emptyset(\vec{P})$ has a non-trivial solution ($\vec{P} = \vec{1}$) if $\rho_i = \frac{\lambda_i}{\mu_i} < 1$, $i = 1,...,K$.

An iterative method, identical to that of the exponential case is used : if at some iteration there exists $j$ such that $\rho_j \ge 1$ then the system is unstable. In that case, the behaviour of the network may be studied by a similar method as for the exponential case.

Let $h$ denote the function, solution of the diffusion equation [3] for the $G/G/1/N$ queue (equivalent to the closed two-server system with $N$ customers). The expected number of packets in queue $i$ is obtained by the integration of $h$ between $0$ and $M_i$. By removal of subscript $i$, we get :

$$ N = \frac{\lambda_D(0)}{b(1-p)} \left[ \frac{1}{\gamma} - \frac{1}{2} \right] + \frac{\lambda_D(0)}{b\gamma(1-p)} \left( M - 1 - \frac{1}{\gamma} \right) \hat{\rho}^{M-1} (1-\hat{\rho}^{-1}) \\ - \frac{\tilde{\mu}_p}{b} \left[ \frac{\hat{\rho}-1}{\gamma} \left( \frac{1}{\gamma} - M+1 \right) + \frac{1}{\gamma} \left( M - \frac{1}{\gamma} \right) \right] + \frac{\tilde{\mu}_p}{2b} (2M-1) + M_p \quad (26) $$
with

$$b = \frac{\lambda}{1-p} Ka + \tilde{\mu} \tilde{K}s, \quad a = \frac{\lambda}{1-p} \tilde{\mu}, \quad \gamma = \frac{2b}{a}, \quad \tilde{\rho} = e^{\gamma},$$

$$p = \frac{\tilde{\rho}(1-\tilde{\rho})}{\tilde{\rho}^{-(M-1)}-\tilde{\rho}^2}, \quad p(0) = \frac{1-\tilde{\rho}}{1-\tilde{\rho}^2 \tilde{\rho}^{M-1}}.$$

The expected response time is obtained by Little's formula : $W = N/\lambda$.

## Amelioration of the solution

An assumption that may be criticized in the previous model is the fact that a packet has always the same probability to be rejected independantly of the number of times it has been rejected. We can assume, and the simulations prove the validity of this assumption, that the probability of rejection for the first access is the probability of saturation of the next queue. For the second access and eventually the following ones the probability of rejection is certainly different.

Let us call $r_{ij}$ the probability that a packet going from i to j will be rejected after a first rejection. We show in appendix 2 that under exponential and Poisson assumptions :

- $r_{ij}$ is independent of the number of previous rejections

- $r_{ij}$ is a function of $\tilde{\mu}_i, \tilde{\mu}_j, \lambda_j, \lambda_i, q_{ij}$.

The equivalent service for a customer going from i to j can be represented by the scheme of figure 3.

<figure><img src="image_7.png" /></figure>

The new equivalent service has the following parameters :

$$\tilde{\mu}_i = \left[ \frac{1}{\mu_i} \sum_{j=1}^{K+1} q_{ij} \left( 1 + \frac{p_j}{1-r_{ij}} \right) \right]^{-1} \quad (27)$$

$$\tilde{K}s_i = -1 + \tilde{\mu}_i^2 \sum_j \frac{q_{ij}[K_i^j + 1]}{(\nu_i^j)^2} \quad (28)$$

with $K_i^j = p_j + Ks_i(1-p_j)$

$$\nu_i^j = \mu_i / \left( 1 + \frac{p_j}{1 - r_{ij}} \right)$$

The only assumption made here is that the value of $r_{ij}$ is the same as in the exponential case. All numerical results are obtained using formulas (27)(28). Let us
note that the new equivalent service is not exponential for exponentially distributed service times.

## Results and Conclusions

We apply the methods developed in previous section to compute the expected response time of the computer network Cyclade, with 7 nodes (Figure 4). We obtain a network with 23 queues. Sixteen of these queues denoted 1 through 16, are inside the store-end-forward packet switching network. We assume half-duplex links. The expected transmission time is shown by Figure 4. The transmission times corresponds to 4,8 and 9,6 kb/s links with an important overhead and with packets of 1000 bits. Packet routing is assumed fixed and is summarized in Table 1. The routes which are not indicated on Table 1 are chosen by the shortest path. Routing probabilities are computed from the routing data (given in Table 1 and in the text). We assume that external arrivals are identical at each host and equal to $\lambda$ and that all the queues have the same buffer size of customers M. The seven hosts represented by seven infinite queues have a service rate equal to 8. We also assume that a packet arriving at a host from the outside of the network has equal probability of having any of the other six nodes as a final destination.

We have summarized in these tables (2 to 4) some results obtained with the above methodology. In these tables we compare the analytical results with simulations. Confidence intervals (95 %) of these simulations have been obtained using the techniques described by Leroudier and Parent in [21].

Tables 2, 3 and 4 give the mean number in each queue (N), the probability of saturation of each buffer (P(M)) and the probability P(O) that each queue is empty under exponential assumptions. These tables compare the results obtained

1. by simulation

2. by the analytic method

3. by the classical BCMP method [1] assuming infinite buffers.

The parameters are

$$ \begin{array}{l} \underline{\text{Table 2}} \\ \underline{\text{Table 3}} \\ \lambda = 1 \\ \lambda = 1.5 \end{array} \qquad \left. \begin{array}{l} M = 4 \\ M = 4 \end{array} \right\} \text{exponential service times} $$

In table 4 we present results obtained with Erlang 2 and Hyperexponential 2 service times.

We can conclude that our results are very accurate for exponential service times, even for heavy traffic. In tables 2, and 3 the analytical results are always within the confidence interval of the simulations. Table 4 shows the limit of this study. It seems that results are underevaluate for Erlang distributions and overevaluate for hyperexponential ones.

Results obtained by increasing M are the same as results found using BCMP [1] (from M = 20 results are exactly identical).

Comparison between infinite buffer models and finite buffer model shows very significant difference and it is surely not possible to model a finite buffer system with an infinite model.
<figure><img src="image_8.png" /><figcaption>Figure 4</figcaption></figure>

<table><caption>Table 1</caption><thead><tr><td>Source</td><td>Destination</td><td>Via</td></tr></thead><tbody><tr><td>A'</td><td>F'</td><td>C</td></tr><tr><td>B'</td><td>F'</td><td>C</td></tr><tr><td>C'</td><td>G'</td><td>F</td></tr><tr><td>D'</td><td>G'</td><td>F</td></tr><tr><td>G'</td><td>C', D'</td><td>B</td></tr><tr><td>F'</td><td>A', B'</td><td>C</td></tr></tbody></table>
Paramètres M = 4 λ = 1

<table><thead><tr><th rowspan="2">File</th><th colspan="3">Simulation</th><th colspan="3">Analytique</th><th>BCMP</th></tr><tr><th>P(O) ±</th><th>P(M) ±</th><th>n̄ ±</th><th>P(O)</th><th>P(M)</th><th>n̄</th><th>n̄</th></tr></thead><tbody><tr><td>1</td><td>.684 .007</td><td>.011 0.002</td><td>.466 .018</td><td>0.677</td><td>0.015</td><td>.490</td><td>.428</td></tr><tr><td>2</td><td>.700 .005</td><td>.008 0.001</td><td>.425 .012</td><td>0.699</td><td>.009</td><td>.416</td><td>.428</td></tr><tr><td>3</td><td>.493 .008</td><td>.053 .004</td><td>.942 .027</td><td>.485</td><td>.051</td><td>.948</td><td>1.</td></tr><tr><td>4</td><td>.865 .004</td><td>.000 .000</td><td>.156 .006</td><td>.869</td><td>.001</td><td>.153</td><td>.154</td></tr><tr><td>5</td><td>.595 .006</td><td>.022 .002</td><td>.656 .017</td><td>.596</td><td>.021</td><td>.637</td><td>.667</td></tr><tr><td>6</td><td>.602 .008</td><td>.019 .003</td><td>.633 .019</td><td>.600</td><td>.020</td><td>.624</td><td>.667</td></tr><tr><td>7</td><td>.738 .006</td><td>.004 .001</td><td>.351 .012</td><td>.732</td><td>.006</td><td>.356</td><td>.464</td></tr><tr><td>8</td><td>.598 .008</td><td>.019 .003</td><td>.638 .020</td><td>.598</td><td>.021</td><td>.630</td><td>.667</td></tr><tr><td>9</td><td>.593 .008</td><td>.024 .003</td><td>.668 .024</td><td>.588</td><td>.027</td><td>.683</td><td>.667</td></tr><tr><td>10</td><td>.721 .006</td><td>.006 .001</td><td>.388 .013</td><td>.725</td><td>.008</td><td>.381</td><td>.364</td></tr><tr><td>11</td><td>.866 .004</td><td>.000 .000</td><td>.154 .006</td><td>.866</td><td>.001</td><td>.156</td><td>.154</td></tr><tr><td>12</td><td>.730 .007</td><td>.005 .001</td><td>.369 .014</td><td>.725</td><td>.008</td><td>.381</td><td>.364</td></tr><tr><td>13</td><td>.868 .005</td><td>.000 .000</td><td>.151 .006</td><td>.867</td><td>.001</td><td>.153</td><td>.154</td></tr><tr><td>14</td><td>.733 .006</td><td>.004 .001</td><td>.360 .012</td><td>.733</td><td>.006</td><td>.356</td><td>.364</td></tr><tr><td>15</td><td>.720 .008</td><td>.006 .001</td><td>.392 .016</td><td>.719</td><td>.009</td><td>.397</td><td>.364</td></tr><tr><td>16</td><td>.866 .005</td><td>.000 .000</td><td>.154 .007</td><td>.867</td><td>.001</td><td>.154</td><td>.154</td></tr></tbody></table>

Table 2

Paramètres M = 4 λ = 1.5

<table><thead><tr><th rowspan="2">File</th><th colspan="3">Simulation</th><th colspan="3">Analytique</th><th>BCMP</th></tr><tr><th>P(O) ±</th><th>P(M) ±</th><th>n̄ ±</th><th>P(O)</th><th>P(M)</th><th>n̄</th><th>n̄</th></tr></thead><tbody><tr><td>1</td><td>.290 .018</td><td>.311 .028</td><td>1.995 .013</td><td>.297</td><td>.255</td><td>1.904</td><td>.818</td></tr><tr><td>2</td><td>.553 .007</td><td>.031 .003</td><td>.764 .021</td><td>.550</td><td>.030</td><td>.749</td><td>.818</td></tr><tr><td>3</td><td>.171 .010</td><td>.342 .016</td><td>2.376 .062</td><td>.165</td><td>.327</td><td>2.383</td><td>3.000</td></tr><tr><td>4</td><td>.797 .006</td><td>.003 .001</td><td>.263 .011</td><td>.800</td><td>.003</td><td>.250</td><td>.250</td></tr><tr><td>5</td><td>.386 .009</td><td>.101 .008</td><td>1.316 .038</td><td>.381</td><td>.093</td><td>1.291</td><td>1.5</td></tr><tr><td>6</td><td>.410 .009</td><td>.087 .007</td><td>1.213 .034</td><td>.400</td><td>.076</td><td>1.194</td><td>1.5</td></tr><tr><td>7</td><td>.602 .008</td><td>.023 .003</td><td>.641 .023</td><td>.600</td><td>.021</td><td>.629</td><td>.66</td></tr><tr><td>8</td><td>.395 .009</td><td>.100 .007</td><td>1.268 .038</td><td>.391</td><td>.086</td><td>1.247</td><td>1.5</td></tr><tr><td>9</td><td>.308 .012</td><td>.215 .017</td><td>1.745 .066</td><td>.318</td><td>.220</td><td>1.776</td><td>1.5</td></tr><tr><td>10</td><td>.531 .011</td><td>.055 .008</td><td>.883 .040</td><td>.544</td><td>.054</td><td>.866</td><td>.66</td></tr><tr><td>11</td><td>.793 .006</td><td>.002 .001</td><td>.263 .011</td><td>.795</td><td>.004</td><td>.266</td><td>.25</td></tr><tr><td>12</td><td>.542 .011</td><td>.050 .006</td><td>.835 .037</td><td>.544</td><td>.055</td><td>.869</td><td>.66</td></tr><tr><td>13</td><td>.799 .006</td><td>.002 .001</td><td>.255 .010</td><td>.800</td><td>.003</td><td>.249</td><td>.25</td></tr><tr><td>14</td><td>.597 .008</td><td>.021 .003</td><td>.649 .022</td><td>.600</td><td>.020</td><td>.629</td><td>.66</td></tr><tr><td>15</td><td>.414 .016</td><td>.157 .018</td><td>1.363 .076</td><td>.436</td><td>.137</td><td>1.324</td><td>.66</td></tr><tr><td>16</td><td>.799 .006</td><td>.003 .001</td><td>.261 .010</td><td>.800</td><td>.004</td><td>.262</td><td>.25</td></tr></tbody></table>

Table 3
<table><caption>Paramètres M = 4, λ = 1</caption><thead><tr><th rowspan="2">File</th><th colspan="2">Erlang 2</th><th colspan="2">Hyper 2</th></tr><tr><th>simulation<br>n ±</th><th>Analytique<br>n</th><th>Simulation<br>n ±</th><th>Analytique<br>n</th></tr></thead><tbody><tr><td>1</td><td>.334 .006</td><td>.302</td><td>.762 .049</td><td>.817</td></tr><tr><td>2</td><td>.375 .008</td><td>.308</td><td>.491 .019</td><td>.577</td></tr><tr><td>3</td><td>.756 .016</td><td>.671</td><td>1.230 .049</td><td>1.347</td></tr><tr><td>4</td><td>.143 .004</td><td>.115</td><td>.164 .009</td><td>0.220</td></tr><tr><td>5</td><td>.549 .010</td><td>.469</td><td>.788 .025</td><td>.889</td></tr><tr><td>6</td><td>.563 .013</td><td>.477</td><td>.713 .027</td><td>.826</td></tr><tr><td>7</td><td>.326 .008</td><td>.268</td><td>.432 .019</td><td>.492</td></tr><tr><td>8</td><td>.559 .013</td><td>.476</td><td>.750 .031</td><td>.857</td></tr><tr><td>9</td><td>.493 .010</td><td>.444</td><td>.872 .041</td><td>1.033</td></tr><tr><td>10</td><td>.305 .007</td><td>.260</td><td>.481 .024</td><td>.588</td></tr><tr><td>11</td><td>.138 .004</td><td>.114</td><td>.172 .009</td><td>.232</td></tr><tr><td>12</td><td>.302 .006</td><td>.260</td><td>.498 .028</td><td>.589</td></tr><tr><td>13</td><td>.146 .004</td><td>.114</td><td>.158 .008</td><td>.222</td></tr><tr><td>14</td><td>.325 .008</td><td>.264</td><td>.394 .016</td><td>.496</td></tr><tr><td>15</td><td>.301 .006</td><td>.263</td><td>.591 .042</td><td>.637</td></tr><tr><td>16</td><td>.140 .004</td><td>.114</td><td>.164 .010</td><td>.225</td></tr></tbody></table>

<caption>Table 4</caption>

We can estimate average network delay using Kleinrock's formula [16]. The average source-to-destination packet delay T is defined as follows :

$$T = \sum_{i=1}^{K} \lambda_{i} W_{i} / \quad \text{where} \quad \lambda_{j} = \lambda_{j}^{\text{ext}}$$

In figure 5 we compare in two cases the source-to-destination packet delay with $Ka = 1$. and $Ks = 1$, for different values of the $\lambda_j^{ext}$ s :

- end-to-end packet delay, response time of hosts included

- end-to-end packet delay, without response times of hosts.

We assume $M = 4$.

In figure 6 we compare the source-to-destination packet delay, with $Ka = 1$ and $Ks = 1$, for different values of the limited capacity M. We assume for each

$$j = 1,...,7, \lambda_j^{\text{ext}} = 1.5.$$

The method developed in this paper may be applied to study the different mechanisms available for controlling flows in a network. Such list of mechanisms used to control flow may be found in Pouzin [18] or Chou and Gerla [19]. This new analytical tool can be still improved to include protocols, acknowledgements, etc...

# APPENDIX 1

## Proof of the uniqueness of the solution in the diffusion case

As in the exponential case, we show that the system with only one queue has an unique solution if $\rho = \lambda/\tilde{\mu} < 1$.

Our system of equations is the following :

$$p = \frac{\tilde{\rho}(1-\tilde{\rho})}{\tilde{\rho}^{-(M-1)}-\tilde{\rho}^2} \quad (29)$$

$$\tilde{\rho} = \frac{\tilde{\rho}}{\tilde{\mu}(1-p)} = \frac{\rho}{1-p} \quad (30)$$
<figure><img src="image_9.png" /><figcaption>Figure 5</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 6</figcaption></figure>
$$ \hat{\rho} = e^{-2} \frac{1-\tilde{\rho}}{\tilde{\rho}Ka+Ks} \qquad (31) $$

From (29) and (30) we obtain:

$$ 1-p = \frac{1-\tilde{\rho}\hat{\rho}^{M-1}}{1-\tilde{\rho}^2 \hat{\rho}^{M-1}} \qquad (32) $$

$$ 1-p = \frac{\tilde{\rho}}{\hat{\rho}}. \qquad (33) $$

Aside from the trivial solution $p=1$, $\rho = +\infty$, there exists a solution only if

$$ f(\tilde{\rho}) = \frac{\tilde{\rho} - \tilde{\rho}^2 \hat{\rho}^{M-1}}{1 - \tilde{\rho}^2 \hat{\rho}^{M-1}} = \frac{\lambda}{\tilde{\mu}} $$

We note that $f(0) = 0$, $f(+\infty) = 1$, and the function $f(\tilde{\rho})$ is an increasing function on $[0, +\infty[$. The system of equations (29) to (31) possesses solutions only if $\rho = \lambda/\tilde{\mu} < 1$.

The vectorial function $\emptyset$ satisfying $\vec{P} = \emptyset(\vec{P})$ is a monotone non-decreasing function.

The uniqueness of the solution of this system may be proved as in the exponential case.

# APPENDIX 2

## Computation of coefficients $r_{ij}$

Let us consider the scheme of figure 7.

<figure><img src="image_11.png" /><figcaption>Figure 7</figcaption></figure>

This scheme corresponds to the condition of queues i and j with

$$ - \lambda = \lambda_i q_{ij} $$

$$ - \lambda_1 = \lambda_j - \lambda_i q_{ij} $$

$$ - \nu = \mu_i $$

$$ - \mu = \tilde{\mu}_j $$

Let $\alpha_i$ = probability {after the next departure of service v the second queue will be full knowing that there are i customers in the second queue}.

Obviously $r_{ij} = \alpha_M$.
The $\alpha_i$ must satisfy the following equations :

$$ (1) \qquad \alpha_0 = \frac{\lambda_1}{\nu+\lambda_1} \alpha_1 $$

$$ (i-1) \qquad \alpha_i = \frac{\lambda_i}{\lambda_1 + \nu + \mu} \alpha_{i+1} + \frac{\mu}{\lambda_1 + \mu + \nu} \alpha_{i-1} $$

$$ (M) \qquad \alpha_{M-1} = \frac{\lambda_1}{\lambda_1 + \nu + \mu} \alpha_M + \frac{\mu}{\lambda_1 + \mu + \nu} \alpha_{M-2} $$

$$ (M+1) \qquad \alpha_M = \frac{\nu}{\nu+\mu} + \frac{\mu}{\nu+\mu} \alpha_{M-1} $$

To solve this system we can use a generating function satisfying equations (1) to (M) and use equation (M+1) as a normalizing equation.

Using $f(x) = \sum_{n=0}^{\infty} \alpha_i x^n$, we obtain :

$$ f(x) = [\lambda_1 - (\lambda_1 + \mu + \nu)x + \mu x^2] = \alpha_0 [\lambda_1 - \mu x] $$

with $\rho_1$ and $\rho_2$ being the solutions of $\lambda_1 - (\lambda_1 + \mu + \nu)x + \mu x^2 = 0$ we obtain :

$$ \frac{f(x)}{\alpha_0} = \frac{A}{\rho_1^{-x}} + \frac{B}{\rho_2^{-x}} $$

with

$$ A = \frac{-\rho_1 + \frac{\lambda_1}{\mu}}{\rho_2 - \rho_1} , \quad B = \frac{\rho_2 - \frac{\lambda_1}{\mu}}{\rho_2 - \rho_1} $$

This formula can be developed in a sum of power of x which gives $\alpha_M$ and $\alpha_{M-1}$ as function of $\alpha_0$. $\alpha_0$ is determined using equation (M+1). After computation we obtain:

$$ \alpha_M = \left\{ \frac{\nu+\mu}{\nu} - \frac{\lambda_1}{\nu} \frac{(\rho_2^M - \rho_1^M) - (\rho_2^{M-1} - \rho_1^{M-1})}{(\rho_2^{M+1} - \rho_1^{M+1}) - (\rho_2^M - \rho_1^M)} \right\}^{-1} $$

and we have shown that $r_{ij} = \alpha_M$.

We can remark that if $\lambda_1 = 0$, then we obtain the classical result : $\alpha_M = \frac{\nu}{\nu+\mu}$.
REFERENCES

[ 1 ] Baskett F., Chandy K.M., Muntz R.R., Palacios F.G. - "Open, closed, and mixed networks of queues with different classes of customers" - J.ACM 22, 248, 1975.

[ 2 ] Kobayashi H. - "Application of the diffusion approximation to queueing networks" - J.ACM 21, 316-328, 1974.

[ 3 ] Gelenbe E. - "On approximate computer models" - J.ACM 22, 261-269, 1975.

[ 4 ] Gelenbe E. and Pujolle G. - "Approximation to a single queue in a network" - Acta Informatica 7, 123-136, 1976.

[ 5 ] Chandy K.M., Herzog U., Woo L. - "Approximate analysis of general queueing network" - IBM Research Report RC 4931, 1974.

[ 6 ] Pennotti M. and Schwartz M. - "Congestion control in store and forward tandem links" - IEEE Trans. on Com. 23, 1434-1443, 1975.

[ 7 ] Rudin H. - "An introduction to flow control" - IBM Research Report RZ 756, Zurich, 1976.

[ 8 ] Rubin I. - "Path delays in communication networks" - Appl. Math. Optimization 1, 1974.

[ 9 ] Labetoulle J. and Pujolle G. - "A study of queueing networks with deterministic service and applications to computer networks" - Acta Informatica 7, 183-195, 1976.

[10] Pujolle G. - "Ergodicity conditions and congestion control in computer networks" - IRIA/LABORIA Research Report 190, 1976.

[11] Lam S.S. - "Store-and-forward buffer requirements in a packet switching network" - IEEE Trans. on Com. 24, 394, 1976.

[12] Schweitzer P. and Lam S.S. - "Buffer overflow in a store-and-forward network node" - IBM Research Report RC 5759, Yorktown Heights, 1975.

[13] Lavenberg S.S. - "Stability and maximum departure rate of certain open queueing networks having finite capacity constraints" - IBM Research Report RJ 1625, San Jose, 1975.

[14] Schweitzer P. - "Maximum throughput in finite capacity open queueing networks" - IBM Research Report RC 5996, Yorktown Heights, 1976.

[15] Irland M. and Pujolle G. - "Models of tandem networks of queues motivated by packet switching" - IRIA Research Report - 1977.

[16] Kleinrock L. - "Communication nets-stochastic message flow and delays" - Mc Graw-Hill - New York, 1964.

[17] Reiser M. and Kobayashi H. - "Accuracy of the diffusion approximation for some queueing systems" - IBM J. Res. and Develop. 18, 2, 110-124, 1974.

[18] Pouzin L. - "Flow control in data networks - methods and tools" - International Conference on Computer Communications - Toronto - Canada, 1976.
[19] Chou W. and Gerla M. - "A unified flow and congestion control model for packet networks" - *International Conference on Computer Communications* - Toronto - Canada, 1976.

[20] Chu W.W. - "Buffer behavior for poisson arrivals and multiple synchronous constant outputs" - *IEEE Trans. on Com.* 18, 530-534, 1970.

[21] Leroudier J. and Parent M. - "Discrete event simulation modelling of computer systems for performance evaluation" - *Research Report IRIA/LABORIA* 117, 1976.