# Age of Information in a Decentralized Network of Parallel Queues with
Routing and Packets Losses

Josu Doncel and Mohamad Assaad J. Doncel is with the University of the Basque
Country, Spain, 48940 Leioa, Spain (e-mail:josu.doncel@ehu.eus)M. Assaad is
with CentraleSupélec, Université Paris-Saclay, Laboratoire des Signaux et
Systèmes 91190 Gif-sur-Yvette, France.

###### Abstract

The paper deals with Age of Information (AoI) in a network of multiple sources
and parallel queues with buffering capabilities, preemption in service and
losses in served packets. The queues do not communicate between each other and
the packets are dispatched through the queues according to a predefined
probabilistic routing. By making use of the Stochastic Hybrid System (SHS)
method, we provide a derivation of the average AoI of a system of two parallel
queues (with and without buffer capabilities) and compare the results with
those of a single queue. We show that known results of packets delay in
Queuing Theory do not hold for the AoI. Unfortunately, the complexity of
computing the average AoI using the SHS method increases highly with the
number of queues. We therefore provide an upper bound of the average AoI in a
system of an arbitrary number of $M/M/1/(N+1)*$ queues and show its tightness
in various regimes. This upper bound allows providing a tight approximation of
the average AoI with a very low complexity. We then provide a game framework
that allows each source to determine its best probabilistic routing decision.
By using Mean Field Games, we provide an analysis of the routing game
framework, propose an efficient iterative method to find the routing decision
of each source and prove its convergence to the desired equilibrium.

##  I Introduction

Age of Information (AoI) is a relatively new metric that measures the
freshness of information in the network. AoI is gaining interest in many areas
(e.g. control, communication networks, etc) due to the proliferation of
applications in which a monitor is interested in having timely updates about a
process of interest. As a typical example, AoI can capture the timeliness of
information in a sensor network where the status of a sensor is frequently
monitored. Since its introduction in the seminal papers [1, 2], AoI has
attracted the attention of many researchers in different fields.

A main part of the AoI literature focuses on the computation of the average
AoI and its minimization, where the channel/network in which the updates are
sent to the monitor is modeled as a queueing system. The computation of AoI in
various queueing models have therefore been investigated. For instance, the
authors in [2] considered an M/M/1 queue, an M/D/1 queue or a D/M/1 queue
model, and the authors in [3, 4] studied an M/M/2 queue model. Other queueing
models can also be found in [5, 6, 7, 8, 9]. While in the aforementioned work,
the status updates of the system are assumed to be sent under a predefined
transmission policy, the problem of the design of the update policy has been
considered in some papers as well, e.g. [10], [11],[12]. Furthermore, the
problem of scheduling and random access design with the aim of minimizing the
average age of the network has been considered recently in several papers [13,
14, 15, 16, 17, 18]. Besides, since single server queue models are not
representative of networks in which packets can be sent through multiple
paths, the average AoI has also been studied in networks with parallel servers
[19, 20, 21], or in more complex networks such as multihop systems [22, 23].
For instance, the scheduling of a single packet flow in multi-hop queueing
networks was studied in [22, 23]. It was shown that Preemptive Last Generated
First Served (P-LGFS) policy is age-optimal if service times are i.i.d.
exponentially distributed. In [24], a multihop scenario in which each node is
both a source and a monitor is considered. Fundamental age limits and near
optimal scheduling policies are provided in this work. Recently, the analysis
of AoI in a multihop multicast context has been studied in [25]. For more
detailed and comprehensive review of recent work on AoI, one can refer to
[26].

In this paper, the network model is different from the aforementioned previous
work since we consider multiple sources that can send their status updates
through a system of different parallel queues. The queues are assumed to be
decentralized in the sense that they cannot communicate between each other.
The incoming traffic from each source is dispatched through the parallel
queues according to a predefined probabilistic routing. Note that we also
develop a framework to optimize this probabilistic routing decision as we will
see later on in this paper. Furthermore, we consider a realistic assumption
that the transmissions through the parallel queues are not perfect and that
packets can be lost. This assumption is quite realistic in many scenarios,
e.g. when the transmission arises over wireless links (which induces errors
and hence packet losses) or even in wired networks when the service provider
breaks down. We aim to analyze the average AoI of this system composed of
parallel queues and, for this purpose, we use the stochastic Hybrid Systems
(SHS) method, which is introduced in [27] (we explain it in detail in Section
II). A related work to ours is [21], where the authors use the SHS method to
compute the average AoI for a system formed by multiple sources and an
arbitrary number of homogeneous M/M/1/1 queues (i.e. with no buffer) as well
as two heterogeneous M/M/1/1 queues, where in both cases preemption in service
is allowed. In our work, we compute the average AoI using the SHS method,
including a system formed by two heterogeneous M/M/1/2* queues with preemption
in service and packet losses. Due to the buffering capability at different
queues, the analysis becomes more challenging and complex as compared to [21].
In addition, we assume that queues are decentralized in the sense that they do
not communicate between each other. This makes our model different than [21],
where it is assumed that all the queues know where is the freshest update.
Besides, we provide an upper bound of the average AoI in a system of an
arbitrary number of $M/M/1/(N+1)$ queues. This allows obtaining an
approximation of the AoI with a low complexity. Finally, we provide in this
paper a game framework to optimize the probabilistic routing decision for each
source, which is to the best of our knowledge has not been considered before
in the AoI literature.

| Without losses | With losses
---|---|---
Single Source | (SERVER-ROUTING) and (SERVER-DOUBLE) | (SERVER-ROUTING) and (SERVER-DOUBLE).
| have equal age for $\lambda$ small and large. See Figure 4. | have equal age always. See Figure 6.
Multiple Sources | (SERVER-ROUTING) and (SERVER-DOUBLE) | (SERVER-ROUTING) and (SERVER-DOUBLE).
| have almost equal age. See Figure 5. | have equal age. See Figure 7.

TABLE I: Summary of Average AoI comparison of the systems without buffer (see
Section III-A)

.

| Without losses | With losses
---|---|---
Single Source | QUEUE-ROUTING and QUEUE-DOUBLE | QUEUE-ROUTING and QUEUE-DOUBLE
| have equal age for $\lambda$ small and large. See Figure 11. | have equal age. See Figure 13.
Multiple Sources | QUEUE-ROUTING and QUEUE-HALF | QUEUE-ROUTING and QUEUE-DOUBLE
| have equal age for $\lambda$ small. | have equal age always. See Figure 14.
| QUEUE-ROUTING and QUEUE-DOUBLE |
| have equal age for $\lambda$ large. See Figure 12. |
TABLE II: Summary of Average AoI comparison of the systems with buffer (see
Section III-B).

The main contributions of this work are twofold. First, we consider a system
with multiple sources where the packets in service can be lost and preemption
is allowed. The packets are sent to the parallel queues according to a
predefined probabilistic routing. We compute the average AoI of a system with
two parallel queues and we compare its average AoI with that of a single
queue. On one hand, in Table I, we present the obtained results where we
compare the following systems: (i) two parallel M/M/1/1 queues, each of them
with arrival rate $\lambda/2$, service rate $\mu$ and loss rate $\theta/2$,
(denoted by SERVER-ROUTING), (ii) one M/M/1/1 queue with arrival rate
$\lambda/2$, loss rate $\theta/2$ and service rate $\mu$ (denoted by SERVER-
HALF) and (iii) one M/M/1/1 queue with arrival rate $\lambda$, loss rate
$\theta$ and service rate $2\mu$ (denoted by SERVER-DOUBLE). On the other
hand, in Table II, we show results in which we compare the following systems:
(i) two parallel M/M/1/2* queues, each of them with arrival rate $\lambda/2$,
service rate $\mu$ and loss rate $\theta$ (denoted by QUEUE-ROUTING), (ii) one
M/M/1/3* queue with arrival rate $\lambda/2$, service rate $\mu$ and loss rate
$\theta/2$ (denoted by QUEUE-HALF) and (iii) one M/M/1/3* queue with arrival
rate $\lambda$, service rate $2\mu$ and loss rate $\theta$ (denoted by QUEUE-
DOUBLE). The description of M/M/1/3* and M/M/1/2* queues will be provided in
Section III. The main conclusion of this part of the work is that the known
results of packet delay in Queuing Theory do not hold for the average AoI. For
instance, we know that the delay of the systems SERVER-ROUTING and SERVER-HALF
is the same, whereas according to our results, the average AoI of SERVER-
ROUTING is smaller. This property also holds when we compare the systems
QUEUE-ROUTING and QUEUE-HALF. Besides, we also conclude that the average AoI
of SERVER-ROUTING is very close to the average AoI of SERVER-DOUBLE and also
that the average AoI of QUEUE-ROUTING is very close to the average AoI of
QUEUE-DOUBLE.
Besides, since the complexity of computing the exact average AoI with SHS
method increases hugely with the number of parallel queues, the second
contribution of this work consists of providing an upper-bound on the average
AoI of a system composed of multiple sources with an arbitrary number of
parallel M/M/1/(N+1)* queues. We also study numerically the accuracy of the
upper bound and we conclude that when the arrival rate is large or when there
are multiple sources, the upper bound is very tight. The interest of this
upper bound lies in the fact that it allows obtaining the average AoI with a
low complexity. The last contribution of this work consists in using the
derived upper bound of the average AoI in order to optimize the probabilistic
routing decision. For instance, we formulate a distributed framework where
each source optimizes its own routing decision using Game Theory. By using
Mean Field Games, we then provide a modification/simplification of the
framework and derived a simple iterative algorithm allowing each source to
find separately its own routing decision. We also provide a theoretical proof
of the convergence of the iterative algorithm to the desired fixed point
(Equilibrium) of the game.

The rest of the paper is organized as follows. In Section II, we formulate the
problem of calculating the average AoI and we present how the SHS can be used.
In Section III we focus on the average AoI derivation in the different systems
under consideration. We present the upper bound of the AoI in Section IV and,
finally, we provide the main conclusion of our work in Section V.

##  II AoI and SHS

We consider a transmitter sending status updates to a monitor. Packet $i$ is
generated at time $s_{i}$ and is received by the monitor at time
$s_{i}^{\prime}$. Hence, we define by $N(t)$ the index of the last received
update at time $t$, i.e., $N(t)=\max\\{i|s_{i}^{\prime}\leq t\\}$, and the
time stamp of the last received update at time $t$ as $U(t)=s_{N(t)}$. The
AoI, or simply the age, is defined as

$\Delta(t)=t-U(t).$

We are interested in calculating the average of the stochastic process
$\Delta(t)$, that is, the average age, which is defined as

$\Delta=\lim_{\tau\to\infty}\frac{1}{\tau}\int_{0}^{\tau}\Delta(t)dt.$

The computation of the average age in a general setting is known to be a
challenging task since the random variables of the interarrival times and of
the system times are dependent. To overcome this difficulty, the authors in
[27] introduce the SHS. For completeness, we describe hereinafter this method
and, for further details one can refer to [28].

In SHS, the system is modeled as a hybrid state $(q(t),{\mathbf{x}}(t))$,
where $q(t)$ is a state of a continuous time Markov Chain and
${\mathbf{x}}(t)$ is a vector whose component belong to $\mathbb{R}^{+}_{0}$
and captures the evolution of the age in the system.

A link $l$ of the Markov Chain represents a transition from two states $q$ and
$q^{\prime}$ with rate $\lambda^{l}$. The interest of SHS is that each
transition $l$ implies a reset mapping in the continuous process
${\mathbf{x}}$. In other words, in each transition $l$, the vector
${\mathbf{x}}$ is transformed to ${\mathbf{x}}^{\prime}$ using a linear
mapping where transformation matrix is given by ${\mathbf{A}}_{l}$, that is,
we have the following SHS transition for every $l$:
$\mathbf{x}^{\prime}={\mathbf{x}}{\mathbf{A}}_{l}$. Throughout this paper, we
denote by $x_{i}^{\prime}$ the i-th element of the vector
$\mathbf{x}^{\prime}$.

Furthermore, each state of the Markov Chain represents the elements of the
system whose age increases at unit rate. In other words, for each state $q$,
we define ${\mathbf{b}}_{q}$ as the vector whose elements are zero or one.
Besides, the evolution of the vector ${\mathbf{x}}(t)$ for state $q$ is given
by $\dot{\mathbf{x}}(t)={\mathbf{x}}{\mathbf{b}}_{q}$.

We assume the Markov Chain is ergodic and we denote by $\pi_{q}$ the
stationary distribution of state $q$. Let ${\mathcal{L}}_{q}$ the set of links
that get out of state $q$ and ${\mathcal{L}}^{\prime}_{q}$ the set of links
that get into state $q$. The following theorem allows us to characterize the
average AoI:

######  Theorem 1 ([27, Thm 4]).

Let $v_{q}(i)$ denote the $i$-th element of the vector ${\mathbf{v}}_{q}$. For
each state $q$, if ${\mathbf{v}}_{q}$ is a non-negative solution of the
following system of equations

${\mathbf{v}}_{q}\sum_{l\in{\mathcal{L}}_{q}}\lambda^{l}={\mathbf{b}}_{q}\pi_{q}+\sum_{l\in{\mathcal{L}}^{\prime}_{q}}\lambda^{l}{\mathbf{v}}_{q_{l}}{\mathbf{A}}_{l},$
(1)

then the average AoI is $\Delta=\sum_{q}v_{q}(0)$.

In the following section, we use the above result to characterize the average
AoI of several systems. In Section IV, we show that the method under
consideration can be also used to obtain an upper-bound on the average AoI of
very complex systems.

##  III Average AoI of Routing Systems Versus of a Single Queue

In this section, we aim to study the average AoI for different configurations
using the SHS method. We first focus on a system formed by queues without
buffer and then consider several cases of queues with buffer. Furthermore, we
consider in this section two main scenarios: i) system with single queue, and
ii) system with multiple parallel queues. In the latter, we consider that
multiple sources are dispatching their packets through the different parallel
queues according to a predefined probabilistic routing. This kind of routing
policies is used in practice and is widely considered in routing literature
since it can be implemented without knowing the instantaneous states of the
network or of the servers (this assumption is realistic as in real-life
networks such information cannot be known at the sources). In more detail, the
routing policy can be explained as follows. Each source $i$ dispatches its
packets according to the following policy: each job/packet of the source is
routed to queue $j$ with probability $p_{ij}$. We can see then that the
arrival rate from source $i$ to queue $j$ is $\lambda_{i}p_{ij}$. In addition,
it is obvious to see that $\sum_{j=1}^{K}p_{ij}=1$. Besides, we consider also
that the transmission through the queues is not reliable and that the loss
rate of each queue $j$ is denoted by $\theta_{j}$. We also assume that the
queues are decentralized in the sense that they do not communicate between
each other.

###  III-A Queues Without Buffer

In this section, we study the age when the queues do not have a buffer. We
first focus on a system formed by an M/M/1/1 queue and then in a system with
two parallel M/M/1/1 queues.

####  III-A1 The M/M/1/1 queue

We consider a system formed by one M/M/1/1 queue that receives traffic from
different $n$ sources when preemption of the packets that are in service is
permitted. We therefore consider a Poisson arrival process from each source
and hence the resulting arrival from all sources is a Poisson process and two
update arrivals cannot occur simultaneously. We are interested in calculating
the average AoI of any source. Without loss of generality, we focus on source
$1$. Thus, we consider that updates arrive to the system according to a
Poisson process, where, without loss of generality, the rate of updates of
source $1$ are denoted by $\lambda_{1}$ and of the rest of the sources
$\sum_{k>1}\lambda_{k}$. The total arrival rate in the system is denoted by
$\lambda$, i.e., $\lambda=\sum_{k=1}^{n}\lambda_{k}$. We assume that the
service time is exponentially distributed with rate $\mu$. We also assume that
an update that is in service is lost with an exponential time of rate
$\theta$.

We use the SHS method to compute the age of this system. The continuous state
of the SHS is ${\mathbf{x}}(t)=[x_{0}(t)\ x_{1}(t)]$, where the correspondence
between $x_{i}$(t) and the elements of the system is as follows: $x_{0}$ is
the age at the monitor and $x_{1}$ the generation time of the packet in
service. One can notice that $x_{0}(t)=x_{1}(t)$ if the packet/update in
service is delivered. The discrete state of the SHS is a two-state Markov
Chain, where $0$ represents that the system is empty and $1$ that an update is
being executed. As explained in the previous section, in each transition in
the Markov chain, the continuous state of the SHS ${\mathbf{x}}$ changes to
${\mathbf{x}}^{\prime}$. The Markov Chain is represented in Figure 1 and the
SHS transitions are given in Table III in which we state explicitly the
transitions from ${\mathbf{x}}$ to ${\mathbf{x}}^{\prime}$.

$l$ | $q_{l}\to q_{l^{\prime}}$ | $\lambda^{l}$ | $x^{\prime}=xA_{l}$ | $\bar{v}_{q_{l}}A_{l}$
---|---|---|---|---
0 | $0\to 1$ | $\lambda_{1}$ | $[x_{0}\ 0]$ | $[v_{0}(0)\ 0]$
1 | $0\to 1$ | $\sum_{k>1}\lambda_{k}$ | $[x_{0}\ x_{0}]$ | $[v_{0}(0)\ v_{0}(0)]$
2 | $1\to 0$ | $\mu$ | $[x_{1}\ 0]$ | $[v_{1}(1)\ 0]$
3 | $1\to 0$ | $\theta$ | $[x_{0}\ 0]$ | $[v_{1}(0)\ 0]$
4 | $1\to 1$ | $\lambda_{1}$ | $[x_{0}\ 0]$ | $[v_{1}(0)\ 0]$
5 | $1\to 1$ | $\sum_{k>1}\lambda_{k}$ | $[x_{0}\ x_{0}]$ | $[v_{1}(0)\ v_{1}(0)]$
TABLE III: Table of transitions of Figure 1 $0$$1$$4$$5$$0$$1$$2$$3$ Figure 1:
The SHS Markov Chain for the one M/M/1/1 queue system with multiple sources
and losses of packets in service.

We now explain each transition $l$:

  *     * $l=0$

There is an arrival of source $1$ and the queue is idle. Therefore, the age of
the monitor does not change, i.e., $x_{0}^{\prime}=x_{0}$ and the age of the
packet in service is zero since there is a fresh update arrived.

    * $l=1$

There is an arrival of one of the others sources when the queue is idle. Since
we are interested in the age of source 1, an arrival from the rest of the
sources does not bring a fresh update of the status of source $1$ and hence it
changes the value of $x_{1}$ to the age of the monitor, that is, the age of
the update in service satisfies $x_{1}^{\prime}=x_{0}$, where $x_{0}$ is the
age of the monitor. Therefore, when this update ends its service, the age of
the monitor remains unchanged.

    * $l=2$

The update under execution ends its service and the age of the monitor is
updated by that of this update, i.e., $x_{0}^{\prime}=x_{1}$.

    * $l=3$

The update that is in service is lost and, therefore, the age of the monitor
does not change.

    * $l=4$

There is an arrival of source $1$ when the queue is in service. For this case,
the update in service is replaced by the fresh one and, therefore, the age of
the monitor does not change, i.e., $x_{0}^{\prime}=x_{0}$, but the age of the
update in service changes to zero.

    * $l=5$

There is an arrival of another source when the queue is in service. For this
case, the update in service is replaced by the fresh one and the age of the
update in service changes to that of the monitor, i.e.,
$x_{1}^{\prime}=x_{0}$.

The stationary distribution of the Markov Chain of Figure 1 is

$\pi_{0}=\frac{\mu+\theta}{\lambda+\mu+\theta},\
\pi_{1}=\frac{\lambda}{\lambda+\mu+\theta}.$

Besides, for the state $q=0$, we have that $\mathbf{b}_{0}=[1,0]$ since the
age of the monitor is the only one that grows at unit rate and the age of the
update in service is irrelevant, whereas for the state $q=1$ we have that
$\mathbf{b}_{1}=[1,1]$ and the age of the monitor and of the update in service
grow at a unit rate. On the other hand, we have that
$\mathbf{v}_{0}=[v_{0}(0)\ v_{0}(1)]$ and $\mathbf{v}_{1}=[v_{1}(0)\
v_{1}(1)]$. From Theorem 4 of [27], we know that the age of this system is
$v_{0}(0)+v_{1}(0)$, where

$\displaystyle\lambda\mathbf{v}_{0}=$ $\displaystyle[\pi_{0}\ 0]+\mu[v_{1}(1)\
0]+\theta[v_{1}(0)\ 0],$ $\displaystyle(\lambda+\mu+\theta)\mathbf{v}_{1}=$
$\displaystyle[\pi_{1}\ \pi_{1}]+\lambda_{1}[v_{1}(0)\
0]+\sum_{k>1}\lambda_{k}[v_{1}(0)\ v_{1}(0)]$
$\displaystyle+\lambda_{1}[v_{0}(0)\ 0]+\sum_{k>1}\lambda_{k}[v_{0}(0)\
v_{0}(0)],$

The above expressions can be written as the following system of equations:

$\displaystyle\lambda v_{0}(0)$ $\displaystyle=\pi_{0}+\mu v_{1}(1)+\theta
v_{1}(0),$ (2a) $\displaystyle(\mu+\theta)v_{1}(0)$
$\displaystyle=\pi_{1}+\lambda v_{0}(0),$ (2b)
$\displaystyle(\lambda+\mu+\theta)v_{1}(1)$
$\displaystyle=\pi_{1}+\sum_{k>1}\lambda_{k}v_{1}(0)+\sum_{k>1}\lambda_{k}v_{0}(0).$
(2c)

The solution to the above system of equations is

$v_{0}(0)=\frac{1}{\lambda_{1}}+\frac{\theta}{\lambda_{1}\mu}-\frac{\pi_{1}}{\lambda+\mu+\theta}$
$v_{1}(0)=\frac{\lambda}{\lambda_{1}\mu}+\frac{\pi_{1}}{\lambda+\mu+\theta}$
$v_{1}(1)=\frac{\sum_{k>1}\lambda_{k}}{\lambda_{1}\mu}+\frac{\pi_{1}}{\lambda+\mu+\theta}.$

Therefore, since, from (1), the average AoI for this case is
$v_{0}(0)+v_{1}(0)$, the next result follows:

######  Proposition 1.

The average AoI of source $1$ in the aforementioned system is given by

$\frac{1}{\lambda_{1}}+\frac{\theta}{\lambda_{1}\mu}+\frac{\lambda}{\lambda_{1}\mu}.$

.

######  Remark 1.

We remark that, when $\theta=0$, (2) coincides with the result of Theorem 2a)
in [27]. In their model, they consider a Markov Chain with a single state, but
when there are updates that are lost their model cannot be considered.
Therefore, we believe that the model presented above is the simplest one to
study the average AoI using the SHS method when there are update losses.

####  III-A2 Two parallel M/M/1/1 queues

We now consider a system formed by two parallel M/M/1/1 queues receiving
traffic from different $n$ sources and where preemption of packets in service
is permitted. We aim to calculate the average age of information of source
$1$. As in the previous case, the arrivals are Poisson and the rate of source
$1$ is denoted by $\lambda_{1}$ and that of the rest of the sources is denoted
by $\sum_{k>1}\lambda_{k}$. Hence, $\lambda=\sum_{k=1}^{n}\lambda_{k}.$ The
packets are dispatched according to a predefined probabilistic routing, where
$p_{1j}$ (resp. $p_{kj}$) is the probability that a job of source $1$ (resp.
of another source $k\neq 1$) is routed to queue $j$. We denote by
$\lambda_{1}p_{11}$ and by $\lambda_{1}p_{12}$ the arrival rates from source
$1$ to queue 1 and to queue 2, respectively. Likewise,
$\sum_{k>1}\lambda_{k}p_{k1}$ and $\sum_{k>1}\lambda_{k}p_{k2}$ denotes the
arrival rates from the rest of the sources to queue 1 and to queue 2,
respectively. The service rate and the loss rate in queue $j$ are respectively
$\mu_{j}$ and $\theta_{j}$, where $j=1,2$. We assume that the queues are
decentralized in the sense that they do not communicate between each other.

######  Remark 2.

The latter assumption makes the model under study here different than [21],
where it is assumed that the servers know where is the freshest update.

We compute the average AoI using the SHS method. First, we define the
continuous state as $\mathbf{x}(t)=[x_{0}(t)\ x_{1}(t)\ x_{2}(t)]$, where the
correspondence between $x_{i}$(t) and each element in the system is as
follows: $x_{0}$ is the age of the monitor and $x_{1}$ (resp. $x_{2}$) is the
age if an update of queue 1 (resp. of queue 2) is delivered. The discrete
state is a Markov Chain with four states, represented in Figure 2 and where
state $k_{1}k_{2}$ represents that in queue 1 there are $k_{1}$ updates, with
$k_{1}\in\\{0,1\\}$, and in queue 2 there are $k_{2}$ updates, with
$k_{2}\in\\{0,1\\}$. We also represent the SHS transitions in Table IV.

$00$$01$$10$$11$$4$$5$$10$$11$$16$$17$$22$$23$$0$$1$$2$$3$$6$$7$$8$$9$$12$$13$$14$$15$$18$$19$$20$$21$ Figure 2: The SHS Markov Chain for system with two parallel M/M/1/1 queues with multiple sources and losses of packets in service. $l$ | $q_{l}\to q_{l^{\prime}}$ | $\lambda^{l}$ | $x^{\prime}=xA_{l}$ | $\bar{v}_{q_{l}}A_{l}$
---|---|---|---|---
0 | $00\to 01$ | $\lambda_{1}p_{12}$ | $[x_{0}\ 0\ 0]$ | $[v_{00}(0)\ 0\ 0]$
1 | $00\to 01$ | $\sum_{k>1}\lambda_{k}p_{k2}$ | $[x_{0}\ 0\ x_{0}]$ | $[v_{00}(0)\ 0\ v_{00}(0)]$
2 | $01\to 00$ | $\mu_{2}$ | $[x_{2}\ 0\ 0]$ | $[v_{01}(2)\ 0\ 0]$
3 | $01\to 00$ | $\theta_{2}$ | $[x_{0}\ 0\ 0]$ | $[v_{01}(0)\ 0\ 0]$
4 | $01\to 01$ | $\lambda_{1}p_{12}$ | $[x_{0}\ 0\ 0]$ | $[v_{01}(0)\ 0\ 0]$
5 | $01\to 01$ | $\sum_{k>1}\lambda_{k}p_{k2}$ | $[x_{0}\ 0\ x_{0}]$ | $[v_{01}(0)\ 0\ v_{01}(0)]$
6 | $00\to 10$ | $\lambda_{1}p_{11}$ | $[x_{0}\ 0\ 0]$ | $[v_{00}(0)\ 0\ 0]$
7 | $00\to 10$ | $\sum_{k>1}\lambda_{k}p_{k1}$ | $[x_{0}\ x_{0}\ 0]$ | $[v_{00}(0)\ v_{00}(0)\ 0]$
8 | $10\to 00$ | $\mu_{1}$ | $[x_{1}\ 0\ 0]$ | $[v_{10}(1)\ 0\ 0]$
9 | $10\to 00$ | $\theta_{1}$ | $[x_{0}\ 0\ 0]$ | $[v_{10}(0)\ 0\ 0]$
10 | $10\to 10$ | $\lambda_{1}p_{11}$ | $[x_{0}\ 0\ 0]$ | $[v_{10}(0)\ 0\ 0]$
11 | $10\to 10$ | $\sum_{k>1}\lambda_{k}p_{k1}$ | $[x_{0}\ x_{0}\ 0]$ | $[v_{10}(0)\ v_{10}(0)\ 0]$
12 | $01\to 11$ | $\lambda_{1}p_{11}$ | $[x_{0}\ 0\ x_{2}]$ | $[v_{01}(0)\ 0\ v_{01}(2)]$
13 | $01\to 11$ | $\sum_{k>1}\lambda_{k}p_{k1}$ | $[x_{0}\ x_{0}\ x_{2}]$ | $[v_{01}(0)\ v_{01}(0)\ v_{01}(2)]$
14 | $11\to 01$ | $\mu_{1}$ | $[x_{1}\ 0\ x_{2}]$ | $[v_{11}(1)\ 0\ v_{11}(2)]$
15 | $11\to 01$ | $\theta_{1}$ | $[x_{0}\ 0\ x_{2}]$ | $[v_{11}(0)\ 0\ v_{11}(2)]$
16 | $11\to 11$ | $\lambda_{1}p_{11}$ | $[x_{0}\ 0\ x_{2}]$ | $[v_{11}(0)\ 0\ v_{11}(2)]$
17 | $11\to 11$ | $\sum_{k>1}\lambda_{k}p_{k1}$ | $[x_{0}\ x_{0}\ x_{2}]$ | $[v_{11}(0)\ v_{11}(0)\ v_{11}(2)]$
18 | $10\to 11$ | $\lambda_{1}p_{12}$ | $[x_{0}\ x_{1}\ 0]$ | $[v_{10}(0)\ v_{10}(1)\ 0]$
19 | $10\to 11$ | $\sum_{k>1}\lambda_{k}p_{k2}$ | $[x_{0}\ x_{1}\ x_{0}]$ | $[v_{10}(0)\ v_{10}(1)\ v_{10}(0)]$
20 | $11\to 10$ | $\mu_{2}$ | $[x_{2}\ x_{1}\ 0]$ | $[v_{11}(2)\ v_{11}(1)\ 0]$
21 | $11\to 10$ | $\theta_{2}$ | $[x_{0}\ x_{1}\ 0]$ | $[v_{11}(0)\ v_{11}(1)\ 0]$
22 | $11\to 11$ | $\lambda_{1}p_{12}$ | $[x_{0}\ x_{1}\ 0]$ | $[v_{11}(0)\ v_{11}(1)\ 0]$
23 | $11\to 11$ | $\sum_{k>1}\lambda_{k}p_{k2}$ | $[x_{0}\ x_{1}\ x_{0}]$ | $[v_{11}(0)\ v_{11}(1)\ v_{11}(0)]$
TABLE IV: Table of transitions of Figure 2.

We now explain each transition $l$:

  *     * $l=0$

There is an arrival of source $1$ to queue 2, when the system is empty.
Therefore, $x_{0}$ and $x_{1}$ do not change and the age of the update in
service in queue 2 is zero due to a fresh update arrival.

    * $l=1$

There is an arrival of an update of the other sources to queue 2 when it is
idle. Therefore, we set $x_{2}^{\prime}=x_{0}$, which means that, when this
update ends its service, the age of the monitor is again $x_{0}$.

    * $l=2$

The update under execution in queue 2 is delivered and the age of the monitor
is updated by that of this update, i.e., $x_{0}^{\prime}=x_{2}$.

    * $l=3$

The update that is in service in queue 2 is lost and, therefore, the age of
the monitor does not change and queue 2 has no updates.

    * $l=4$

There is an arrival of source $1$ to queue 2 when it is in service. For this
case, since preemption is permitted, the update in service is replaced by the
fresh one and, therefore, the age of the update in service in queue 2 changes
to zero.

    * $l=5$

There is an arrival of another source to queue 2 when it has an update in
service. For this case, the update in service is replaced by the fresh one and
the age of the update in service in queue 1 changes to that of the monitor,
i.e., $x_{2}^{\prime}=x_{0}$.

  * The transitions 6-11 are symmetric to 0-5, respectively.

    * $l=12$

When there is an update in queue 2, if an update of source $1$ arrives to
queue 1, the age of the monitor and of the update in queue 2 do not change,
whereas that of queue 1 is set to zero, that is, $x_{1}^{\prime}=0$.

    * $l=13$

When there is an update in queue 2, if there is an arrival of one of the other
sources in queue 1, the age of the update in queue 1 changes to the age of the
monitor, i.e,. $x_{1}^{\prime}=x_{0}$, whereas $x_{0}$ and $x_{2}$ do not
change.

    * $l=14$

When there are updates in both queues, an update of queue 1 is delivered and
the age of the monitor changes to $x_{1}$, i.e., $x_{0}^{\prime}=x_{1}$.

    * $l=15$

When there are updates in both queues, if an update of queue 1 is lost, the
age of the monitor does not change and queue 1 is idle.

    * $l=16$

When both queues have updates in service, if an update of source $1$ arrives
to queue 1, we set $x_{1}^{\prime}$ to zero and the rest does not change.

    * $l=17$

When both queues have updates in service, if an update of the other sources
arrives to queue 1, we set $x_{1}^{\prime}$ to the same as the monitor.

  * The transitions 18-23 are symmetric to 12-17, respectively.

The stationary distribution of the Markov Chain in Figure 2 is given by

$\pi_{k_{1}k_{2}}=\frac{\rho_{1}^{k_{1}}\rho_{2}^{k_{2}}}{(1+\rho_{1})(1+\rho_{2})},\
\text{ for }k_{1},k_{2}=0,1,$

where
$\rho_{j}=\frac{\lambda_{1}p_{1j}+\sum_{k>1}\lambda_{k}p_{kj}}{\mu_{j}+\theta_{j}},j=1,2$.
Moreover, we define the value of $b_{q}$ for each state
$q\in\\{00,10,01,11\\}$ as follows: $\mathbf{b}_{00}=[1\ 0\ 0]$,
$\mathbf{b}_{10}=[1\ 1\ 0]$, $\mathbf{b}_{01}=[1\ 0\ 1]$ and
$\mathbf{b}_{11}=[1\ 1\ 1]$. We also define, for $q\in\\{00,10,01,11\\}$, the
following vector $\mathbf{v}_{q}=[v_{q}(0)\ v_{q}(1)\ v_{q}(2)]$.

Let $\hat{\mu}=\mu_{1}+\mu_{2}$ and $\hat{\theta}=\theta_{1}+\theta_{2}$. The
SHS method says that the average AoI of this system is $\sum_{q}v_{q}(0)$,
where $v_{q}(0)$ is the solution of the following system of equations:

$\displaystyle\lambda\mathbf{v}_{00}$ $\displaystyle=[\pi_{00}\ 0\
0]+\mu_{1}[v_{10}(1)\ 0\ 0]$ $\displaystyle+\theta_{1}[v_{10}(0)\ 0\ 0]$
$\displaystyle+\mu_{2}[v_{01}(2)\ 0\ 0]+\theta_{2}[v_{01}(0)\ 0\ 0]$ (3)
$\displaystyle(\lambda+\mu_{1}+\theta_{1})\mathbf{v}_{10}$
$\displaystyle=[\pi_{10}\ \pi_{10}\ 0]+\lambda_{1}p_{11}[v_{00}(0)\ 0\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}p_{1k}[v_{00}(0)\ v_{00}(0)\ 0]$
$\displaystyle+\mu_{2}[v_{11}(2)\ v_{11}(1)\ 0]$
$\displaystyle+\theta_{2}[v_{11}(0)\ v_{11}(1)\
0]+\lambda_{1}p_{11}[v_{10}(0)\ 0\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}p_{1k}[v_{10}(0)\ v_{10}(0)\ 0]$ (4)
$\displaystyle(\lambda+\mu_{2}+\theta_{2})\mathbf{v}_{01}$
$\displaystyle=[\pi_{01}\ 0\ \pi_{01}]+\lambda_{1}p_{12}[v_{00}(0)\ 0\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}p_{k2}[v_{00}(0)\ 0\ v_{00}(0)]$
$\displaystyle+\mu_{1}[v_{11}(1)\ 0\ v_{11}(2)]$
$\displaystyle+\theta_{1}[v_{11}(0)\ 0\ v_{11}(2)]$
$\displaystyle+\lambda_{1}p_{12}[v_{01}(0)\ 0\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}p_{k2}[v_{01}(0)\ 0\ v_{01}(0)]$ (5)
$\displaystyle(\lambda+\hat{\mu}+\hat{\theta})\mathbf{v}_{11}$
$\displaystyle=[\pi_{11}\ \pi_{11}\ \pi_{11}]+\lambda_{1}p_{11}[v_{01}(0)\ 0\
v_{01}(2)]$ $\displaystyle+\sum_{k>1}\lambda_{k}p_{k1}[v_{01}(0)\ v_{01}(0)\
v_{01}(2)]$ $\displaystyle+\lambda_{1}p_{11}[v_{11}(0)\ 0\ v_{11}(2)]$
$\displaystyle+\sum_{k>1}\lambda_{k}p_{k1}[v_{11}(0)\ v_{11}(0)\ v_{11}(2)]$
$\displaystyle+\lambda_{1}p_{12}[v_{10}(0)\ v_{10}(1)\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}p_{k2}[v_{10}(0)\ v_{10}(1)\ v_{10}(0)]$
$\displaystyle+\lambda_{1}p_{12}[v_{11}(0)\ v_{11}(1)\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}p_{k2}[v_{11}(0)\ v_{11}(1)\ v_{11}(0)].$
(6)

Since the first equation has two irrelevant variables and the second and third
ones have one irrelevant variable, the above expression can be written
alternatively as a system of 8 equations.

######  Proposition 2.

The average AoI of source $1$ in the aforementioned system is given by
$v_{00}(0)+v_{10}(0)+v_{01}(0)+v_{11}(0)$, where for $q\in\\{00,10,01,11\\}$,
$v_{q}(0)$ is the solution of (3)-(6).

####  III-A3 Average AoI Comparison

(a) M/M/1/1 queue with arrival rate $\lambda/2$.

(b) M/M/1/1 queue with service rate $2\mu$.

(c) Two parallel M/M/1/1 queues.

Figure 3: Representation of the models under comparison in Section III-A3 for
two sources. Figure 4: Average AoI of source $1$ when $\lambda_{1}$ varies
from $0.1$ to $10^{3}$ with a single source and without losses
($\lambda_{k}=0$ for all $k>1$ and $\theta=0$). $\mu=1$. Figure 5:  Average
AoI of source $1$ when $\lambda_{1}$ varies from $0.1$ to $10^{3}$ with
multiple sources and without losses ($\sum_{k>1}\lambda_{k}=10$ and
$\theta=0$). $\mu=1$. Figure 6:  Average AoI of source $1$ when $\lambda_{1}$
varies from $0.1$ to $10^{3}$ with a single source and losses ($\lambda_{k}=0$
for all $k>1$ and $\theta=10$). $\mu=1$.  Figure 7:  Average AoI of source $1$
when $\lambda_{1}$ varies from $0.1$ to $10^{3}$ with multiple sources and
losses ($\sum_{k>1}\lambda_{k}=10$ and $\theta=10$). $\mu=1$.

We now compare the average AoI of source $1$ for the models we have studied in
this section. For this purpose, we consider three systems. The first one
consists of a single M/M/1/1 queue with arrival rate of source $1$
$\lambda_{1}/2$ and of the rest of the sources $(1/2)\sum_{k>1}\lambda_{k}$,
loss rate $\theta/2$ and service rate $\mu$ (see Figure 3(a)). The average AoI
of this model is represented with a solid line. The second system we consider
is a single M/M/1/1 queue with arrival rate of source $1$ $\lambda_{1}$ and of
the rest of the sources $\sum_{k>1}\lambda_{k}$, loss rate $\theta$ and
service rate $2\mu$ (see Figure 3(b)). The average AoI of this model is
represented with a dotted line. Finally, we consider a system with two
parallel M/M/1/1 queues with arrival rate of source $1$ equal to $\lambda_{1}$
and of the rest of the sources $\sum_{k>1}\lambda_{k}$. Besides, we consider
that $p_{kj}=1/2$ for all $k=1,\dots,n$ and $j=1,2$, the loss rate and the
service rate are in both servers $\theta/2$ and $\mu$, respectively (see
Figure 3(c)). The average AoI of this model is represented with a dashed line.
Our goal is to determine which system has the smallest average AoI when
$\lambda_{1}$ varies and the rest of the parameters are fixed. To this end, we
have solved numerically the systems of equations in (2) and in (3)-(6). We set
$\mu=1$ in these simulations. When we study the system with multiple sources,
we consider that $\sum_{k>1}\lambda_{k}=10$, and in the case of losses in the
packets in service, we set $\theta=10$.
In Figure 4, we plot the average AoI of these systems as a function of
$\lambda_{1}$ when there is a single source and there are no losses in the
queues. We observe that the smallest age is achieved for the single M/M/1/1
queue system with service rate $2\mu$. We also observe that the age of the two
parallel M/M/1/1 queues is the same as the latter when $\lambda_{1}$ is either
very small or very large.

In Figures 5-7, we show that the average AoI of a system with two parallel
M/M/1/1 queues is equal to that of a single server with service rate $2\mu$.

In queueing theory, it is known that, among the systems under consideration in
this section, the one that minimizes the delay is the single M/M/1/1 queue
with service rate $2\mu$. Therefore, these illustrations show that the AoI
also verifies this property. On the other hand, for classical queueing theory
metrics such as delay, the performance of two parallel M/M/1/1 queues
coincides with that of a single M/M/1/1 queue with arrival rate $\lambda/2$
and loss rate $\hat{\theta}/2$. However, as far as average AoI is concerned,
one can see that, according to the figures we present in this section, this is
not the case (solid and dashed lines do not coincide on these figures).

###  III-B Queues With Buffer

We now focus on queues with buffer. For this case, an update starts getting
served upon its arrival to a queue, if the queue is idle. However, if the
queue is busy, the incoming update is put in the last position of the queue
and, if the queue is full, the last update of the buffer is replaced by the
new one. In this section, we aim to compare the average AoI of a system with
one queue and buffer size two with that of two parallel queues with buffer
size one. We first compute the average AoI of the former system and then of
the latter one.

####  III-B1 The M/M/1/3* queue

We concentrate on a system formed by a queue with a buffer of size two. When
an update arrives and the system is empty, it gets served immediately.
However, if an update arrives when another packet is in service, it replaces
the last update in the queue if it is full and, otherwise, it is put in the
last position of the queue. This system will be denoted in the remainder of
the paper as the M/M/1/3* queue.

When traffic comes from $n$ different sources, we are interested in computing
the average AoI of source $1$. Updates of source $1$ arrive to the queue
according to a Poisson process of rate $\lambda_{1}$ and of those of the rest
of the sources with rate $\sum_{k>1}\lambda_{k}$. We assume that the updates
that are waiting in the queue are served according to the FCFS discipline and
that the service time is exponentially distributed with rate $\mu$, as well as
the update that is in service is lost with exponentially distributed time with
rate $\theta$.

$0$$1$$2$$3$$12$$13$$0$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$11$ Figure 8: The SHS
Markov Chain for the M/M/1/3* queue with multiple sources and losses of
packets in service.

We employ the SHS method to calculate the average AoI of this system. The
continuous state is given by $\mathbf{x}(t)=[x_{0}(t)\ x_{1}(t)\ x_{2}(t)\
x_{3}(t)]$, where the correspondence between $x_{i}$(t) and each element in
the system is as follows: $x_{0}$ is the age of the monitor, $x_{1}$ is the
age if the update in service is delivered and $x_{2}$ and $x_{3}$ is
respectively the age if the update in the first and second positions of the
queue are delivered. The discrete state is a four state Markov Chain, where
state $k$ represents that there are $k$ updates present in the system, with
$k=0,1,2,3$. The Markov Chain under consideration and the SHS transition are
represented respectively in Figure 8 and Table V.

$l$ | $q_{l}\to q_{l^{\prime}}$ | $\lambda^{l}$ | $x^{\prime}=xA_{l}$ | $\bar{v}_{q_{l}}A_{l}$
---|---|---|---|---
0 | $0\to 1$ | $\lambda_{1}$ | $[x_{0}\ 0\ 0\ 0]$ | $[v_{0}(0)\ 0\ 0\ 0]$
1 | $0\to 1$ | $\sum_{k>1}\lambda_{k}$ | $[x_{0}\ x_{0}\ 0\ 0]$ | $[v_{0}(0)\ v_{0}(0)\ 0\ 0]$
2 | $1\to 0$ | $\mu$ | $[x_{1}\ 0\ 0\ 0]$ | $[v_{1}(1)\ 0\ 0\ 0]$
3 | $1\to 0$ | $\theta$ | $[x_{0}\ 0\ 0\ 0]$ | $[v_{1}(0)\ 0\ 0\ 0]$
4 | $1\to 2$ | $\lambda_{1}$ | $[x_{0}\ x_{1}\ 0\ 0]$ | $[v_{1}(0)\ v_{1}(1)\ 0\ 0]$
5 | $1\to 2$ | $\sum_{k>1}\lambda_{k}$ | $[x_{0}\ x_{1}\ x_{1}\ 0]$ | $[v_{1}(0)\ v_{1}(1)\ v_{1}(1)\ 0]$
6 | $2\to 1$ | $\mu$ | $[x_{1}\ x_{2}\ 0\ 0]$ | $[v_{2}(1)\ v_{2}(2)\ 0\ 0]$
7 | $2\to 1$ | $\theta$ | $[x_{0}\ x_{2}\ 0\ 0]$ | $[v_{2}(0)\ v_{2}(2)\ 0\ 0]$
8 | $2\to 3$ | $\lambda_{1}$ | $[x_{0}\ x_{1}\ x_{2}\ 0]$ | $[v_{2}(0)\ v_{2}(1)\ v_{2}(2)\ 0]$
9 | $2\to 3$ | $\sum_{k>1}\lambda_{k}$ | $[x_{0}\ x_{1}\ x_{2}\ x_{2}]$ | $[v_{2}(0)\ v_{2}(1)\ v_{2}(2)\ v_{2}(2)]$
10 | $3\to 2$ | $\mu$ | $[x_{1}\ x_{2}\ x_{3}\ 0]$ | $[v_{3}(1)\ v_{3}(2)\ v_{3}(3)\ 0]$
11 | $3\to 2$ | $\theta$ | $[x_{0}\ x_{2}\ x_{3}\ 0]$ | $[v_{3}(0)\ v_{3}(2)\ v_{3}(3)\ 0]$
12 | $3\to 3$ | $\lambda_{1}$ | $[x_{0}\ x_{1}\ x_{2}\ 0]$ | $[v_{3}(0)\ v_{3}(1)\ v_{3}(2)\ 0]$
13 | $3\to 3$ | $\sum_{k>1}\lambda_{k}$ | $[x_{0}\ x_{1}\ x_{2}\ x_{2}]$ | $[v_{3}(0)\ v_{3}(1)\ v_{3}(2)\ v_{3}(2)]$
TABLE V: Table of transitions of Figure 8

We now explain each transition $l$:

  *     * $l=0$

The system is empty and an update of source $1$ arrives. The age of the
monitor is not modified and we set $x_{1}^{\prime}=0$.

    * $l=1$

The system is empty and an update of another source arrives. The age of the
monitor is not modified and the age $x_{1}$ changes to $x_{0}$, i.e.,
$x_{1}^{\prime}=x_{0}$.

    * $l=2$

When there is an update getting in service and the queue is empty. If the
update in service is delivered, the age of the monitor changes to $x_{1}$,
i.e., $x_{0}^{\prime}=x_{1}$.

    * $l=3$

The queue is empty and the update in service is lost. For this case, the age
of the monitor does not change and the age of $x_{1}$ is replaced by zero.

    * $l=4$

The queue is busy and an update of source $i$ arrives. The age of the monitor
and that of $x_{1}$ are not modified and we set $x_{2}^{\prime}=0$.

    * $l=5$

The queue is busy and an update of source $i$ arrives. The age of the monitor
and that of $x_{1}$ are not modified and the age $x_{2}$ changes to $x_{1}$,
i.e., $x_{2}^{\prime}=x_{1}$.

    * $l=6$

There are two updates in the system and the update in service is delivered
and, therefore, the age of the monitor changes to $x_{1}$ and the age $x_{2}$
to $x_{1}$, i.e., $x_{0}^{\prime}=x_{1}$ and $x_{1}^{\prime}=x_{2}$.

    * $l=7$

There are two updates in the system and the update in service is lost. For
this case, the age of the monitor does not change, but the age $x_{1}$ is
replaced by $x_{2}$, i.e., $x_{1}^{\prime}=x_{2}$ since the update that was
waiting start getting served.

    * $l=8$

There are two updates in the system and an update of source $1$ arrives. The
ages of the updates that are present in the system do not change and we set
$x_{3}^{\prime}=0$.

    * $l=9$

There are two updates in the system and an update of another source arrives.
The ages of the updates that are present in the system do not change and the
age $x_{3}$ changes to $x_{2}$, i.e., $x_{3}^{\prime}=x_{2}$.

    * $l=10$

The system is full and the update in service is delivered. For this case, the
age of the monitor changes to $x_{1}$, the age of $x_{1}$ to $x_{2}$ and the
age of $x_{2}$ to $x_{3}$.

    * $l=11$

The system is full and the update in service is lost. For this case, the age
of the monitor does not change, but the age of $x_{1}$ changes to $x_{2}$ and
the age of $x_{2}$ to $x_{3}$.

    * $l=12$

The system is full and an update of source $1$ arrives. The age of the monitor
and of that of $x_{1}$ and $x_{2}$ are not modified and we set
$x_{3}^{\prime}=0$.

    * $l=13$

The system is full and an update of another source arrives. The ages of the
monitor, of $x_{1}$ and of $x_{2}$ do not change, but the age of $x_{3}$ is
set to $x_{2}$, i.e., $x_{3}^{\prime}=x_{2}$.

Let $\lambda=\sum_{k=1}^{n}\lambda_{k}$ and $\rho=\frac{\lambda}{\mu+\theta}$.
The stationary distribution of the Markov Chain in Figure 8 is

$\pi_{j}=\frac{\rho^{j}}{1+\rho+\rho^{2}+\rho^{3}},\ j=0,1,2,3.$

We now define the vector $\mathbf{b}_{q}$ for all state $q$ of the Markov
Chain of Figure 8: $\mathbf{b}_{0}=[1\ 0\ 0\ 0]$, $\mathbf{b}_{1}=[1\ 1\ 0\
0]$, $\mathbf{b}_{2}=[1\ 1\ 1\ 0]$ and $\mathbf{b}_{3}=[1\ 1\ 1\ 1]$. Besides,
for all state $q\in\\{0,1,2,3\\}$, $\mathbf{v}_{q}=[v_{q}(0)\ v_{q}(1)\
v_{q}(2)\ v_{q}(3)]$. From Theorem 4 in [27], we have that the average AoI of
the M/M/1/3* queue is $v_{0}(0)+v_{1}(0)+v_{2}(0)+v_{3}(0)$, where $v_{q}(0)$
is the solution of the following system of equations:

$\displaystyle\lambda\mathbf{v}_{0}=$ $\displaystyle[\pi_{0}\ 0\ 0\
0]+\mu[v_{1}(1)\ 0\ 0\ 0]+\theta[v_{1}(0)\ 0\ 0\ 0]$ (7)
$\displaystyle(\lambda+\mu+\theta)\mathbf{v}_{1}=$ $\displaystyle[\pi_{1}\
\pi_{1}\ 0\ 0]+\lambda_{1}[v_{0}(0)\ 0\ 0\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}[v_{0}(0)\ v_{0}(0)\ 0\ 0]$
$\displaystyle+\mu[v_{2}(1)\ v_{2}(2)\ 0\ 0]$ $\displaystyle+\theta[v_{2}(0)\
v_{2}(2)\ 0\ 0].$ (8) $\displaystyle(\lambda+\mu+\theta)\mathbf{v}_{2}=$
$\displaystyle[\pi_{2}\ \pi_{2}\ \pi_{2}\ 0]+\lambda_{1}[v_{1}(0)\ v_{1}(1)\
0\ 0]$ $\displaystyle+\sum_{k>1}\lambda_{k}[v_{1}(0)\ v_{1}(1)\ v_{1}(1)\ 0]$
$\displaystyle+\mu[v_{3}(1)\ v_{3}(2)\ v_{3}(3)\ 0]$
$\displaystyle+\theta[v_{3}(0)\ v_{3}(2)\ v_{3}(3)\ 0].$ (9)
$\displaystyle(\lambda+\mu+\theta)\mathbf{v}_{3}=$ $\displaystyle[\pi_{3}\
\pi_{3}\ \pi_{3}\ \pi_{3}]+\lambda_{1}[v_{2}(0)\ v_{2}(1)\ v_{2}(2)\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}[v_{2}(0)\ v_{2}(1)\ v_{2}(2)\ v_{2}(2)]$
$\displaystyle+\lambda_{1}[v_{3}(0)\ v_{3}(1)\ v_{3}(2)\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}[v_{3}(0)\ v_{3}(1)\ v_{3}(2)\ v_{3}(2)]$
(10)

Since the first equation has three irrelevant variables and the second and
third equations have respectively two and one irrelevant variables, the above
expression can be alternatively written as a system of 10 equations.

######  Proposition 3.

The average AoI of source $1$ in the aforementioned system is given by
$v_{0}(0)+v_{1}(0)+v_{2}(0)+v_{3}(0)$, where for $q\in\\{0,1,2,3\\}$,
$v_{q}(0)$ is the solution of (7)-(10).

####  III-B2 Two parallel M/M/1/2* queues

We consider a system formed by two parallel queues with buffer size equal to
one and $n$ different sources. The packets are dispatched according to a
predefined probabilistic routing. An update of source $1$ arrives to the
system with rate $\lambda_{1}$ and it is sent to queue $j$ with probability
$p_{1j}$. Therefore, the arrival rate of source $1$ to queue $j$ is
$\lambda_{1}p_{1j}$. Besides, an update of the rest of the sources arrives to
the system with rate $\sum_{k>1}\lambda_{k}$ and an update of source $k\neq i$
is sent to queue $j$ with probability $p_{kj}$. Therefore, the arrival rate of
the rest of the sources to queue $j$ is $\sum_{k>1}\lambda_{k}p_{kj}$. If an
update finds the queue full, it replaces the update that is waiting in the
queue, whereas when the queue is idle, it starts being served immediately.
This system will be denoted as two parallel M/M/1/2* queues.

We assume that the service time of queue $j$ is exponentially distributed with
rate $\mu_{j}$ and that updates/packets in service in queue $j$ are lost with
exponential time of rate $\theta_{j}$, $j=1,2$. We assume that the queues are
decentralized in the sense that they do not communicate between each other.

We seek to compute the average AoI of this system using the SHS method. The
continuous state is $\mathbf{x}(t)=[x_{0}(t)\ x_{11}(t)\ x_{12}(t)\ x_{21}(t)\
x_{22}(t)]$, where the correspondence between $x_{i}$(t) and each element is
as follows: $x_{0}$ is the age of the monitor, $x_{j1}$ is the age if the
update/packet in service in queue $j$ is delivered and $x_{j2}$ the age if the
update that is waiting for service in queue $j$ is delivered. The discrete
state is described by a Markov Chain, where the state $k_{1}k_{2}$ denotes
that there are $k_{1}$ updates in queue 1 and $k_{2}$ in queue 2, with
$k\in\\{0,1,2\\}$. The Markov Chain we study is depicted in Figure 9. We note
that some of the links are unified to avoid heavy notation. The SHS
transitions for this model are reported in Appendix C.

$00$$01$$02$$10$$20$$11$$21$$12$$22$$4$$15$$22$$25$$28$$29$$0$$1$$2$$3$$5$$6$$7$$8$$9$$10$$11$$12$$13$$14$$16$$17$$18$$19$$20$$21$$24$$23$$27$$26$
Figure 9: The SHS Markov Chain for system with two parallel M/M/1/2* queues
with multiple sources and losses of packets in service.

Let
$\rho_{1}=\frac{\lambda_{1}p_{11}+\sum_{k>1}\lambda_{k}p_{k1}}{\mu_{1}+\theta_{1}}$
and $\rho_{2}=\frac{\lambda_{1}p_{12}+\lambda_{k}p_{k2}}{\mu_{2}+\theta_{2}}.$
The stationary distribution of the Markov Chain in Figure 9 is

$\pi_{k_{1}k_{2}}=\frac{\rho_{1}^{k_{1}}\rho_{2}^{k_{2}}}{(1+\rho_{1}+\rho_{1}^{2})(1+\rho_{2}+\rho_{2}^{2})},k_{1},k_{2}\in\\{0,1,2\\}.$

Let $\hat{\mu}=\mu_{1}+\mu_{2}$, $\hat{\theta}=\theta_{1}+\theta_{2}$ and
$\mathcal{Q}=\\{00,10,20,01,11,21,02,12,22\\}$. For every $q\in\mathcal{Q}$,
we define $\mathbf{v}_{q}=[v_{q}(0)\ v_{q}(1)\ v_{q}(2)\ v_{q}(3)\ v_{q}(4)]$
and the vector $\mathbf{b}_{q}$ as follows: $\mathbf{b}_{00}=[1\ 0\ 0\ 0\ 0]$,
$\mathbf{b}_{10}=[1\ 1\ 0\ 0\ 0]$, $\mathbf{b}_{20}=[1\ 1\ 1\ 0\ 0]$,
$\mathbf{b}_{01}=[1\ 0\ 0\ 1\ 0]$, $\mathbf{b}_{11}=[1\ 1\ 0\ 1\ 0]$,
$\mathbf{b}_{21}=[1\ 1\ 1\ 1\ 0]$, $\mathbf{b}_{02}=[1\ 0\ 0\ 1\ 1]$,
$\mathbf{b}_{12}=[1\ 1\ 0\ 1\ 1]$ and $\mathbf{b}_{22}=[1\ 1\ 1\ 1\ 1]$. We
use the result of Theorem 4 in [27] that shows that the average AoI of this
system is given by $\sum_{q\in\mathcal{Q}}v_{q}(0),$ where $v_{q}(0)$ is the
solution to the following system of equations:

$\displaystyle\lambda\mathbf{v}_{00}=$ $\displaystyle[\pi_{00}\ 0\ 0\ 0\
0]+\mu_{1}[v_{10}(1)\ 0\ 0\ 0\ 0]$ $\displaystyle+\theta_{1}[v_{10}(0)\ 0\ 0\
0\ 0]+\mu_{2}[v_{01}(3)\ 0\ 0\ 0\ 0]$ $\displaystyle+\theta_{2}[v_{01}(0)\ 0\
0\ 0\ 0]$ (11) $\displaystyle(\lambda+\mu_{1}+\theta_{1})\mathbf{v}_{10}=$
$\displaystyle[\pi_{10}\ \pi_{10}\ 0\ 0\ 0]$
$\displaystyle+\lambda_{1}p_{11}[v_{00}(0)\ 0\ 0\ 0\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}p_{k1}[v_{00}(0)\ v_{00}(0)\ 0\ 0\ 0]$
$\displaystyle+\mu_{2}[v_{11}(3)\ v_{11}(1)\ 0\ 0\ 0]$
$\displaystyle+\theta_{2}[v_{11}(0)\ v_{11}(1)\ 0\ 0\ 0]$
$\displaystyle+\mu_{1}[v_{20}(1)\ v_{20}(2)\ 0\ 0\ 0]$
$\displaystyle+\theta_{1}[v_{20}(0)\ v_{20}(2)\ 0\ 0\ 0]$
$\displaystyle+\mu_{2}[v_{11}(3)\ v_{11}(1)\ 0\ 0\ 0]$
$\displaystyle+\theta_{2}[v_{11}(0)\ v_{11}(1)\ 0\ 0\ 0]$ (12)
$\displaystyle(\lambda+\mu_{1}+\theta_{1})\mathbf{v}_{20}=$
$\displaystyle[\pi_{20}\ \pi_{20}\ \pi_{20}\ 0\ 0]$
$\displaystyle+\mu_{2}[v_{21}(3)\ v_{21}(1)\ v_{21}(2)\ 0\ 0]$
$\displaystyle+\theta_{2}[v_{21}(0)\ v_{21}(1)\ v_{21}(2)\ 0\ 0]$
$\displaystyle+\lambda_{1}p_{11}[v_{10}(0)\ v_{10}(1)\ 0\ 0\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}p_{k1}[v_{10}(0)\ v_{10}(1)\ v_{10}(1)\ 0\
0]$ $\displaystyle+\lambda_{1}p_{11}[v_{20}(0)\ v_{20}(1)\ 0\ 0\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}p_{k1}[v_{20}(0)\ v_{20}(1)\ v_{20}(1)\ 0\
0]$ (13) $\displaystyle(\lambda+\mu_{2}+\theta_{2})\mathbf{v}_{01}=$
$\displaystyle[\pi_{01}\ 0\ 0\pi_{01}\ 0]$
$\displaystyle+\lambda_{1}p_{12}[v_{00}(0)\ 0\ 0\ 0\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}p_{k2}[v_{00}(0)\ 0\ 0v_{00}(0)\ 0]$
$\displaystyle+\mu_{2}[v_{02}(3)\ 0\ 0\ v_{20}(4)\ 0]$
$\displaystyle+\theta_{2}[v_{02}(0)\ 0\ 0\ v_{20}(4)\ 0]$
$\displaystyle+\mu_{1}[v_{11}(1)\ 0\ 0\ v_{11}(3)\ 0]$
$\displaystyle+\theta_{1}[v_{11}(0)\ 0\ 0\ v_{11}(3)\ 0]$ (14)
$\displaystyle(\lambda+\hat{\mu}+\hat{\theta})\mathbf{v}_{11}=$
$\displaystyle[\pi_{11}\ \pi_{11}\ 0\ \pi_{11}\ 0]$
$\displaystyle+\lambda_{1}p_{11}[v_{01}(0)\ 0\ 0\ v_{01}(3)\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}p_{k1}[v_{01}(0)\ v_{01}(0)\ 0\ v_{01}(3)\
0]$ $\displaystyle+\lambda_{1}p_{12}[v_{10}(0)\ v_{10}(1)\ 0\ 0\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}p_{k2}[v_{10}(0)\ v_{10}(1)\ 0v_{10}(0)\
0]$ $\displaystyle+\mu_{2}[v_{12}(3)\ v_{12}(1)\ 0\ v_{21}(4)\ 0]$
$\displaystyle+\theta_{2}[v_{12}(0)\ v_{12}(1)\ 0\ v_{21}(4)\ 0]$
$\displaystyle+\mu_{1}[v_{21}(1)\ v_{21}(2)\ 0\ v_{21}(3)\ 0]$
$\displaystyle+\theta_{1}[v_{21}(0)\ v_{21}(2)\ 0\ v_{21}(3)\ 0]$ (15)
$\displaystyle(\lambda+\hat{\mu}+\hat{\theta})\mathbf{v}_{21}=$
$\displaystyle[\pi_{21}\ \pi_{21}\ \pi_{21}\pi_{21}\ 0]$ $\displaystyle+$
$\displaystyle\lambda_{1}p_{11}[v_{11}(0)\ v_{11}(1)\ 0\ v_{11}(3)\ 0]$
$\displaystyle+$ $\displaystyle\sum_{k>1}\lambda_{k}p_{k1}[v_{11}(0)\
v_{11}(1)\ v_{11}(1)\ v_{11}(3)\ 0]$ $\displaystyle+$
$\displaystyle\lambda_{1}p_{11}[v_{21}(0)\ v_{21}(1)\ 0\ v_{21}(3)\ 0]$
$\displaystyle+$ $\displaystyle\sum_{k>1}\lambda_{k}p_{k1}[v_{21}(0)\
v_{21}(1)\ v_{21}(1)\ v_{21}(3)\ 0]$ $\displaystyle+$
$\displaystyle\lambda_{1}p_{12}[v_{20}(0)\ v_{20}(1)\ v_{20}(2)\ 0\ 0]$
$\displaystyle+$ $\displaystyle\sum_{k>1}\lambda_{k}p_{k2}[v_{20}(0)\
v_{20}(1)\ v_{20}(2)v_{20}(0)\ 0]$ $\displaystyle+$
$\displaystyle\mu_{2}[v_{22}(3)\ v_{22}(1)\ v_{22}(2)\ v_{22}(4)\ 0]$
$\displaystyle+$ $\displaystyle\theta_{2}[v_{22}(0)\ v_{22}(1)\ v_{22}(2)\
v_{22}(4)\ 0]$ (16)
$\displaystyle(\lambda+\mu_{2}+\theta_{2})\mathbf{v}_{02}=$
$\displaystyle[\pi_{02}\ 0\ 0\ \pi_{02}\ \pi_{02}]$
$\displaystyle+\lambda_{1}p_{12}[v_{01}(0)\ 0\ 0\ v_{01}(3)\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}p_{k2}[v_{01}(0)\ 0\ 0\ v_{01}(3)\
v_{01}(3)]$ $\displaystyle+\lambda_{1}p_{12}[v_{02}(0)\ 0\ 0\ v_{02}(3)\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}p_{k2}[v_{02}(0)\ 0\ 0\ v_{02}(3)\
v_{02}(3)]$ $\displaystyle+\mu_{1}[v_{12}(1)\ 0\ 0\ v_{12}(3)\ v_{12}(4)]$
$\displaystyle+\theta_{1}[v_{12}(0)\ 0\ 0\ v_{12}(3)\ v_{12}(4)]$ (17)
$\displaystyle(\lambda+\hat{\mu}+\hat{\theta})\mathbf{v}_{12}$
$\displaystyle=[\pi_{12}\ \pi_{12}\ 0\ \pi_{12}\ \pi_{12}]$
$\displaystyle+\lambda_{1}p_{11}[v_{11}(0)\ v_{11}(1)\ 0\ v_{11}(3)\ 0]$
$\displaystyle+\sum_{k>1}\lambda_{k}p_{k1}[v_{11}(0)\ v_{11}(1)\ 0\ v_{11}(3)\
v_{11}(3)]$ $\displaystyle+\lambda_{1}p_{12}[v_{20}(0)\ v_{20}(1)\ v_{20}(2)\
0\ 0]$ $\displaystyle+\sum_{k>1}\lambda_{k}p_{k2}[v_{20}(0)\ v_{20}(1)\
v_{20}(2)\ v_{20}(0)\ 0]$ $\displaystyle+\mu_{1}[v_{22}(1)\ v_{22}(2)\ 0\
v_{22}(3)\ v_{22}(4)]$ $\displaystyle+\theta_{1}[v_{22}(0)\ v_{22}(2)\ 0\
v_{22}(3)\ v_{22}(4)]$ (18)
$\displaystyle(\lambda+\hat{\mu}+\hat{\theta})\mathbf{v}_{22}$
$\displaystyle=[\pi_{22}\ \pi_{22}\ \pi_{22}\ \pi_{22}\ \pi_{22}]$
$\displaystyle+$ $\displaystyle\lambda_{1}p_{11}[v_{12}(0)\ v_{12}(1)\ 0\
v_{12}(3)\ v_{12}(4)]$ $\displaystyle+$
$\displaystyle\sum_{k>1}\lambda_{k}p_{k1}[v_{12}(0)\ v_{12}(1)\ v_{12}(1)\
v_{12}(3)\ v_{12}(4)]$ $\displaystyle+$
$\displaystyle\lambda_{1}p_{11}[v_{22}(0)\ v_{22}(1)\ 0\ v_{22}(3)\
v_{22}(4)]$ $\displaystyle+$
$\displaystyle\sum_{k>1}\lambda_{k}p_{k1}[v_{22}(0)\ v_{22}(1)\ v_{22}(1)\
v_{22}(3)\ v_{22}(4)]$ $\displaystyle+$
$\displaystyle\lambda_{1}p_{12}[v_{21}(0)\ v_{21}(1)\ v_{21}(2)\ v_{21}(3)\
0]$ $\displaystyle+$ $\displaystyle\sum_{k>1}\lambda_{k}p_{k2}[v_{21}(0)\
v_{21}(1)\ v_{21}(2)\ v_{21}(3)\ v_{21}(3)]$ $\displaystyle+$
$\displaystyle\lambda_{1}p_{12}[v_{22}(0)\ v_{22}(1)\ v_{22}(2)\ v_{22}(3)\
0]$ $\displaystyle+$ $\displaystyle\sum_{k>1}\lambda_{k}p_{k2}[v_{22}(0)\
v_{22}(1)\ v_{22}(2)\ v_{22}(3)\ v_{22}(3)]$ (19)

From the above expressions, if we remove the irrelevant variables, we obtain a
system of 27 equations.

######  Proposition 4.

The average AoI of source $1$ in the aforementioned system is given by
$\sum_{q\in\mathcal{Q}}v_{q}(0)$, where $v_{q}(0)$ is the solution of
(11)-(19).

####  III-B3 Age Comparison

(a) M/M/1/2 queue with arrival rate $\lambda/2$.

(b) M/M/1/2 queue with service rate $2\mu$.

(c) Two parallel M/M/1/2 queues.

Figure 10: Representation of the models under comparison in Section III-B3 for
two sources.

We compare the average AoI of the models presented in this section. We focus
on the following three systems. First, we consider an M/M/1/3* queue with
arrival rate of source $1$ equal to $\lambda_{1}/2$ and that of the rest of
the sources $(1/2)\sum_{k>1}\lambda_{k}$, loss rate $\theta/2$ and service
rate $\mu$ (see Figure 10(a)). The average AoI of this model is represented
with a solid line. The second system we consider is an M/M/1/3* queue with
arrival rate of source $1$ equal to $\lambda_{1}$ and that of the rest of the
sources $\sum_{k>1}\lambda_{k}$, loss rate $\theta$ and service rate $2\mu$
(see Figure 10(b)). The average AoI of this model is represented with a dotted
line. We also consider a system with two parallel M/M/1/2* queues with arrival
rate of source $1$ $\lambda_{1}$ and that of the rest of the sources
$\sum_{k>1}\lambda_{k}$. Each of the servers satisfies that
$p_{k1}=p_{k2}=1/2$ for all $k=1,\dots,n$, has a loss rate equal to $\theta/2$
and a service rate $\mu$ (see Figure 10(c)). The average AoI of the latter
model is represented with a dashed line. We aim to investigate which system
has the smallest average AoI when $\lambda_{1}$ varies. Thus, we have solved
numerically the systems of equations in (7)-(10) and of (11)-(19). We set
$\mu=1$ in these simulations. When we study the system with multiple sources,
we consider that $\sum_{k>1}\lambda_{k}=10$ and in the case of packet losses,
we set $\theta=10$.

We first focus on the average AoI for a single source and when there are no
packet losses. The evolution of the AoI of source $1$ with respect to
$\lambda_{1}$ is represented in Figure 11. We observe that for the M/M/1/3*
queue with service rate $2\mu$, the average AoI coincides with that of the two
parallel M/M/1/2* queues when $\lambda_{1}$ is either very small or very
large. Another interesting property obtained from these simulations is that
the average AoI of the M/M/1/3* queue is not monotone in $\lambda_{1}$. This
phenomenon is due to the FCFS discipline and the presence of the buffer of
size 2 and can be interpreted as follows. For small arrival rates, all the
packets will be delivered directly without staying too much time in the buffer
and the system will behave like an M/M/1/1 queue. Therefore, when the arrival
rate increases and on average there is only one packet (or less) in the
server, the AoI will keep decreasing since more fresh packets improves the
AoI. We can observe in this figure that the minimum AoI is achieved when
$\lambda_{1}=\mu$ (which means that on average we have one packet in the
server as explained above). Then, when the arrival rate keeps increasing,
there will be always packets in the buffer (in addition to the packet in the
server) and the arrived packets will be delayed, which will increase the AoI.
When the arrival rate grows very large, the packet in the second place in the
buffer will be constantly replaced by the new arrived packet. However, there
will be always a delay due the fact that the packet in the first place of the
buffer should wait until the packet in the server is delivered. In other
words, the AoI will converge to an asymptotic value. However, this asymptotic
value is greater than the AoI when $\lambda=\mu$ since in that case there is
on average one packet in the system (i.e. the packet is directly served by the
server) and hence the packets are not delayed by the buffer. In addition to
these results, we provide in Appendix B some results obtained by simulations
for single M/M/1/3* queue. As expected The results assess the accuracy of the
SHS method used to evaluate the average AoI.

In Figure 12, we study the average AoI for a system with multiple sources and
without losses. We see that the AoI of the two parallel M/M/1/2* queues
coincides with that of the M/M/1/3* queue with half traffic rate when
$\lambda$ is small, whereas it coincides with that of the M/M/1/3* queue with
double service rate when $\lambda$ is large. It is worth mentioning that we
consider here that the arrival rate of the rest of the sources is
$\sum_{k>1}\lambda_{k}=10$, which implies that there will be always packets in
the server and in the buffer and the system cannot behave as an M/M/1/1 by
changing $\lambda_{1}$ of the first source. This explains why the average AoI
decreases with $\lambda_{1}$ until reaching a limiting value and the average
AoI does not have the same shape as in Figure 11.

We also study the average AoI with a single source and losses in Figure 13.
For this case, the AoI of the system with two parallel M/M/1/2* queues and of
an M/M/1/3* queue with double service rate coincide when $\lambda_{1}$ is
either small or large. In this case, we can see that the average AoI decreases
with the arrival rate $\lambda_{1}$. This can be explained by the fact that,
since the packets can get lost, it is better from AoI perspective to have more
arrived packets (even if these packets are delayed in the buffer).

Finally, in Figure 14, we show the average age of information for different
values of $\lambda_{1}$ when there are multiple sources and losses. This
illustration presents that, depending on the value of $\lambda_{1}$, the AoI
approaches that of an M/M/1/3* with half arrival rate and half loss rate or
that of an M/M/1/3* with double service rate, as in Figure 14.

The main conclusion of these illustrations is that, from an AoI perspective,
the M/M/1/3* queue with double service rate is the optimal one among the
systems under consideration. Besides, we characterize the instances where the
AoI of the two parallel queues coincides with the optimal AoI.

Figure 11: Average AoI comparison when $\lambda_{1}$ varies from $0.1$ to
$10^{3}$ with a single source and without losses ($\lambda_{k}=0$ for all
$k>1$ and $\theta=0$). $\mu=1$.  Figure 12:  Average AoI comparison when
$\lambda_{1}$ varies from $0.1$ to $10^{3}$ with multiple sources and without
losses ($\sum_{k>1}\lambda_{k}=10$ and $\theta=0$). $\mu=1$. Figure 13:
Average AoI comparison when $\lambda_{1}$ varies from $0.1$ to $10^{3}$ with a
single source and losses ($\lambda_{k}=0$ for all $k>1$ and $\theta=10$).
$\mu=1$. Figure 14:  Average AoI comparison when $\lambda_{1}$ varies from
$0.1$ to $10^{3}$ with multiple sources and losses ($\sum_{k>1}\lambda_{k}=10$
and $\theta=10$). $\mu=1$.

###  III-C Parallel Queues With Buffer Size $N>1$

We now focus on the study of the average AoI in a system with $K$ parallel
queues with buffer size $N>1$. We notice that using the SHS method leads to
the analysis of a Markov Chain with a number of states equal to $K\cdot(N+2)$.
This implies that the number of SHS transitions increases at a very high rate
with the number of queues and with the buffer size. As a result, according to
(1), the number of equations to be solved so as to obtain the AoI suffers from
the curse of dimensionality. Thus, providing an analytical expression of the
AoI of source $1$ in a system with an arbitrary routing system (with an
arbitrary number of queues and an arbitrary buffer size) seems to be
intractable using the considered method. However, as we will see in the next
section, it is possible to provide an upper-bound on the AoI.

##  IV Upper-bound on the Average AoI for an Arbitrary Routing System

We study the average AoI of a system with $K>2$ parallel queues with $N>1$
buffer size. In this section, we provide an upper bound on the age of
information using the SHS method in a system with a single and multiple
sources.

We now explain the system we study here. We consider a system with $n$ sources
where, for all i, the updates of source $i$ and of the rest of the sources
arrive to the system with rate $\lambda_{i}$. We denote by $p_{ij}$ the
probability that a job of source $i$ is routed to server $j$. Hence,
$\lambda_{i}=\sum_{j=1}^{K}\lambda_{i}p_{ij}$. Besides, the total incoming
traffic to the system is denoted by $\lambda$, i.e.,
$\lambda=\sum_{i=1}^{n}\lambda_{i}$. We assume that the service rate in queue
$j$ is exponentially distributed with rate $\mu_{j}$. In the following result,
we provide an upper bound of the average AoI of the system under study here.
Without loss of generality, the result is provided for source $1$, however one
can see that the result can be obtained for any source $i$. The proof of this
result is reported in Appendix A.

######  Theorem 2.

For the aforementioned system, the average AoI of source $1$ is upper bounded
by

$\frac{1}{\sum_{j=1}^{K}\mu_{j}}\left(1+KN+\sum_{j=1}^{K}\frac{\sum_{k>1}\lambda_{k}p_{kj}+\mu_{j}}{\lambda_{j}p_{1j}}\right).$
(20)

It is also important to note that to obtain the above result we consider that
when an update completes the service, we create a _fake_ update to keep the
system full of packets. This is the reason why, unlike in the previous
section, we are not able to study the influence of the packet losses on the
upper bound of the average AoI we provide in Theorem 20. In fact, let us
consider that $N=0$ as an example. In this case, when a packet in service is
lost, we put a false/fake update with the same age of the lost one in service.
Therefore, this fake update modifies the age of the monitor when it is served
and delivered to the monitor. As a result, the fake update will modify the age
at the monitor and the system will behave like a system with no packet losses
(but with a different service rate). This explains why the above result cannot
capture the impact of packet losses.

###  IV-A Tightness of the Upper bound

We now aim to explore if the upper bound on the average AoI is tight for the
systems we have studied in Section III. We consider $\mu=1$ and, when we
analyze the AoI for multiple sources, we fix the arrival rate of the rest of
the sources to $10$, that is, $\sum_{k>1}\lambda_{k}=10$.

We first study in Figures 15-16, a system formed by 2 parallel queues with
equal arrival rate and service rate, i.e., $K=2,$ $N=1$,
$p_{11}=p_{21}=p_{12}=p_{22}=1/2$ and $\mu_{1}=\mu_{2}=\mu$. For this case, we
get from Theorem 20 that the upper bound for source $1$ is

$\frac{1}{2\mu}\left(3+\frac{\sum_{k>1}\lambda_{k}+2\mu}{\lambda_{1}}\right).$

Figure 15: Upper bound and real average AoI comparison when $\lambda$ varies
from $0.1$ to $10^{3}$ for two parallel M/M/1/2* queues with a single source
($\lambda_{2}=0$).  Figure 16: Upper bound and real average AoI comparison
when $\lambda$ varies from $0.1$ to $10^{3}$ for two parallel M/M/1/2* queues
with multiple sources ($\lambda_{2}=10$). $\mu=1$.

As we observe in Figure 15, the upper bound is tight when the arrival rate of
source $1$ is large enough, whereas in Figure 16, we show that it is always
very close to the real age.

We now focus on the influence of $\mu$ on the tightness of the upper bound of
the average AoI in a system with two parallel M/M/1/2* queues with a single
source. First, we consider $\mu=0.1$ in Figure 17 and we show that, when
$\lambda_{1}$ is larger than 2, the upper bound is very tight. Then, we
consider $\mu=10$ in Figure 18 and we show that, when $\lambda_{1}$ is larger
than 10, the upper bound is accurate. Finally, we consider $\mu=100$ in Figure
19 and, as it can be observed in this illustration, the upper bound is very
close to the real age when $\lambda_{1}$ is larger than 1000.

Figure 17: Upper bound and real average AoI comparison when $\lambda$ varies
from $0.1$ to $10^{3}$ for two parallel M/M/1/2* queues with a single source
($\lambda_{2}=0$). $\mu=0.1$.  Figure 18: Upper bound and real average AoI
comparison when $\lambda$ varies from $0.1$ to $10^{3}$ for two parallel
M/M/1/2* queues with a single source ($\lambda_{2}=0$). $\mu=10$.  Figure 19:
Upper bound and real average AoI comparison when $\lambda$ varies from $0.1$
to $10^{3}$ for two parallel M/M/1/2* queues with a single source
($\lambda_{2}=0$). $\mu=100$.

We also investigate a system formed by one M/M/1/3* queue with arrival rate
$\lambda$, service rate $\mu$ and multiple sources in Figure 20. For this
case, we have that $K=1$ and $N=2$ and, therefore, from Theorem 20, the
average AoI of source $1$ is given by

$\frac{1}{\mu}\left(3+\frac{\sum_{k>1}\lambda_{k}+\mu}{\lambda_{1}}\right).$

As we see in Figure 20, we show that it is always very close to the real age
for any value of $\lambda_{1}$. Thus, this plot confirms that the upper bound
on the average AoI we provide in this paper is very tight when the arrival
rate is large.

Figure 20: Upper bound and real average AoI comparison when $\lambda$ varies
from $0.1$ to $10^{3}$ for one M/M/1/3* queue with multiple source
($\lambda_{2}=10$). $\mu=1$.

###  IV-B AoI Comparison with a Single M/M/1/1 queue

We now consider a system with a single source which is formed by $K$
homogeneous queues without buffer, i.e., $\mu_{j}=\mu$ for all $j=1,\dots,K$.
According to the result of Theorem 20, the average AoI is upper bounded by

$\frac{1}{K\mu}\left(1+\mu\sum_{j=1}^{K}\frac{1}{\lambda_{j}}\right).$ (21)

We now aim to compare the above expression with the average AoI of a single
M/M/1/1 queue with preemption of jobs in service, arrival rate $\lambda/K$ and
service rate $\mu$, which according to Theorem 2(a) in [27] is given by

$\frac{K}{\lambda}+\frac{1}{\mu}.$

In the following result, we compare the above expressions.

######  Proposition 5.

Let $K>1$. Then,

$\frac{K}{\lambda}+\frac{1}{\mu}>\frac{1}{K\mu}\left(1+\mu\sum_{j=1}^{K}\frac{1}{\lambda_{j}}\right).$

###### Proof.

First, we note that, when $\lambda_{j}=\lambda/K$ for all $j$, we have that

$\frac{1}{K\mu}\left(1+\mu\sum_{j=1}^{K}\frac{1}{\lambda_{j}}\right)=\frac{1}{K\mu}\left(1+\frac{\mu
K^{2}}{\lambda}\right)=\frac{1}{K\mu}+\frac{K}{\lambda}.$

Therefore, we aim to show that

$\frac{1}{K\mu}+\frac{K}{\lambda}<\frac{1}{\mu}+\frac{K}{\lambda}\iff K>1$

And the desired result follows since the last expression is always true. ∎

An interesting result is derived from the above proposition. Indeed, when we
consider a system formed by K parallel queues and each of them receives the
same arrival rate, since the expression (21) provides an upper bound on the
average AoI, this result implies that the average AoI of a single M/M/1/1
queue with arrival rate $\lambda/K$ is larger than that of the considered
system.

###  IV-C AoI Comparison with [21]

In Theorem 2 in [21], the author provides the following expression of the
average AoI of a system with homogeneous parallel queues where the incoming
jobs are always sent to the server with the oldest job:

$\frac{1}{\mu}\left(\frac{1}{K}\prod_{i=1}^{K-1}\frac{\rho}{i+\rho}+\frac{1}{\rho}+\frac{1}{\rho}\sum_{l=1}^{K-1}\prod_{i=1}^{l}\frac{\rho}{i+\rho}\right).$
(22)

We now notice that, in our model, the knowledge of the queue with the oldest
job is not considered. Therefore, one might expect that the average AoI is
always smaller in the model in [21]. In the following result, we consider the
regime where $\lambda$ tends to infinity and we compare both models.

######  Proposition 6.

When $\lambda\to\infty$, we have that (22) and (20) tend to $\tfrac{1}{K\mu}.$

###### Proof.

The proof is straightforward from (22) and (20). ∎

From this result, we conclude that, when $\lambda\to\infty$, the improvement
of the average AoI caused by the knowledge of the state of the queues is
negligible.

###  IV-D Optimization of the Upper Bound

In this section, we provide a framework to minimize the upper bound of the
average age of each source $i$. The objective is to find the routing
probabilities $p_{ij}$ in such that the age upper bound for each source $i$ is
minimized. The problem can be formulated as a game framework. More formally,
the problem can be formulated as

$\min_{\mathbf{p}_{i}}\frac{1}{\sum_{j=1}^{K}\mu_{j}}\left(1+KN+\sum_{j=1}^{K}\frac{\sum_{l\neq
i}^{n}\lambda_{l}p_{lj}+\mu_{j}}{\lambda_{i}p_{ij}}\right),\ \forall i$ (23)
$\text{s.t.}\ \sum_{j=1}^{K}p_{ij}=1\ \forall i$ (24)

where $\mathbf{p}_{i}=[p_{i1},...,p_{iK}].$ In the sequel, we will
characterizes the Nash Equilibrium (NE) of the above problem and provide a
solution that achieves the NE. Before defining NE, we first introduce the so-
called best response set valued function ($BR_{i}$) for each source (or
player) $i$, which is given as follows

$\text{($BR_{i}$)}:\ \text{Given}\
\mathbf{p}_{-i}\overset{\Delta}{=}(\mathbf{p}_{1},...,\mathbf{p}_{i-1},\mathbf{p}_{i+1},...\mathbf{p}_{N})$
$\mathbf{p}_{i}\in\arg\max_{\mathbf{p}^{\prime}_{i}}\
U_{i}(\mathbf{p}^{\prime}_{i},\mathbf{p}_{-i})$

where $\mathbf{p}^{\prime}_{i}$ satisfies $\sum_{j=1}^{K}p^{\prime}_{ij}=1$
and
$U_{i}(\mathbf{p}^{\prime}_{i},\mathbf{p}_{-i})=\frac{1}{\sum_{j=1}^{K}\mu_{j}}\left(1+KN+\sum_{j=1}^{K}\frac{\sum_{l\neq
i}^{N}\lambda_{l}p_{lj}+\mu_{j}}{\lambda_{i}p_{ij}}\right)$. In order to show
explicitly the dependence of $BR_{i}$ in $\mathbf{p}_{-i}$, we will use the
notation $BR_{i}(\mathbf{p}_{-i})$ to represent the best response set valued
function of source $i$. In other words, the best response consists of
optimizing the utility of each source with respect only to its own action
vector (i.e. routing probability $\mathbf{p}_{i}$).
Furthermore, we can also define the sources’ joint best-response function as

$BR(\mathbf{p})=\left(BR_{1}(\mathbf{p}_{-1}),...BR_{N}(\mathbf{p}_{-N})\right)$

We now provide the definition of NE and the relation with the sources’ best
response.

######  Definition 1.

A strategy profile $\mathbf{p}=(\mathbf{p}_{1},...,\mathbf{p}_{N})$ is a pure
Nash equilibrium iff $\forall i$,

$\forall\ \mathbf{p}^{\prime}_{i},\
U_{i}(\mathbf{p}^{\prime}_{i},\mathbf{p}_{-i})\geq
U_{i}(\mathbf{p}_{i},\mathbf{p}_{-i})$

######  Definition 2.

A strategy profile is a Nash Equilibrium iff

$\mathbf{p}\in BR(\mathbf{p})$

In words, a NE is a fixed point of the BR dynamic. In the sequel, we will
therefore show that a fixed point of BR exists and it is unique.
It is straightforward to see that
$U_{i}(\mathbf{p}^{\prime}_{i},\mathbf{p}_{-i})$ is convex with respect to
$\mathbf{p}_{i}$. The best response problem for each source $i$ can be solved
easily using the standard Lagrangian technique.
The Lagrangian for each source can be written as follows,

$L_{i}(\mathbf{p}_{i},\delta)=U_{i}(\mathbf{p}_{i},\mathbf{p}_{-i})+\delta(\sum_{j=1}^{K}p_{ij}-1)$

The optimal solution of the above optimization problem for each source $i$,
i.e. ($\mathbf{p}^{*}_{i},\delta^{*}$), can then be obtained by KKT and
complementarity conditions, i.e. by using $\frac{\partial L_{i}}{\partial
p_{ij}}=0$ $\forall j$ and $\delta^{*}(\sum_{j=1}^{K}p^{*}_{ij}-1)=0$. After
some algebraic manipulations, this leads to the following expression of
$\mathbf{p}^{*}_{i}$

$p_{ij}^{*}=\frac{\omega_{ij}}{\sum_{j^{\prime}=1}^{K}\omega_{ij^{\prime}}}$

where $\omega_{ij}=\sqrt{\sum_{l\neq i}^{N}\lambda_{l}p_{lj}+\mu_{j}}$.
Consequently, a NE is simply the solution of the following system of equations

$p_{ij}^{*}=\frac{\sqrt{\sum_{l\neq
i}^{N}\lambda_{l}p^{*}_{lj}+\mu_{j}}}{\sum_{j^{\prime}=1}^{K}\sqrt{\sum_{l\neq
i}^{N}\lambda_{l}p^{*}_{lj^{\prime}}+\mu_{j^{\prime}}}},\ \forall i,j$ (25)

We will show numerically later on that the aforementioned system of equations
has a solution. Before that, we will provide an analysis in the case of large
number of sources, by using mean field analysis, and develop a simple
iterative algorithm allowing each source to find its probabilistic routing
vector $\mathbf{p}_{i}=[p_{i1},...,p_{iK}]^{T}$.

####  IV-D1 Mean Field Analysis

We provide here an analysis in the case of large number of sources and provide
a distributed iterative algorithm that converges to the solution of the system
of equations in (25). In order to use mean field tools, indistinguishable
sources should be considered. We therefore consider that
$\lambda_{l}=\bar{\lambda}$. We first define the following mean field term for
each queue $j$:

$m^{n}_{ij}=\frac{1}{n}\sum_{l\neq i}^{n}p_{lj}$

One can notice that $m^{n}_{ij}\rightarrow m_{j}$ when $n\rightarrow\infty$
$\forall$ $i$, where $m_{j}=\lim_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n}p_{ij}$.
The problem in (23) can be written as

$\min_{\mathbf{p}_{i}}\frac{1}{\sum_{j=1}^{K}\frac{\mu_{j}}{n}}\left(\frac{1+KN}{n}+\sum_{j=1}^{K}\frac{\lambda
m^{n}_{ij}+\frac{\mu_{j}}{n}}{\lambda p_{ij}}\right),\ \forall i$

One can notice also that when $n\rightarrow\infty$,
$\mu_{j}/n\rightarrow\bar{\mu}_{j}$, where $\bar{\mu}_{j}$ represents the
scaled asymptotic service rate per source. In fact, when the number of sources
tends to infinity, the service rate of each queue must be high enough to serve
all sources.
When $n$ is very large, the system of equations in (25) tends to:

$p_{ij}^{*}=\frac{\sqrt{m_{j}+\frac{\bar{\mu}_{j}}{\bar{\lambda}}}}{\sum_{j^{\prime}=1}^{K}\sqrt{m_{j^{\prime}}+\frac{\bar{\mu}_{j^{\prime}}}{\bar{\lambda}}}},\
\forall i,j$ (26)

By summing the previous equations over all source indexes $i$ and taking the
limit when $n\rightarrow\infty$, we get

$m_{j}=\frac{\sqrt{m_{j}+\frac{\bar{\mu}_{j}}{\bar{\lambda}}}}{\sum_{j^{\prime}=1}^{K}\sqrt{m_{j^{\prime}}+\frac{\bar{\mu}_{j^{\prime}}}{\bar{\lambda}}}},\
\forall j$ (27)

By making the variable change
$y_{j}=\sqrt{m_{j}+\frac{\bar{\mu}_{j}}{\bar{\lambda}}}$, the previous
equations can be written as:

$y_{j}=\frac{1}{\sum_{j^{\prime}=1}^{K}y_{j^{\prime}}}+\frac{\bar{\mu}_{j}}{\bar{\lambda}y_{j}},\
\forall j$ (28)

######  Proposition 7.

If the system of equations in (28) has a solution, then this solution is
unique.

###### Proof.

Let’s assume that the system of equations has two different solutions
$\mathbf{y}^{1}=[y^{(1)}_{1},...,y^{(1)}_{K}]^{T}$ and
$\mathbf{y}^{2}=[y^{(2)}_{1},...,y^{(2)}_{n}]^{T}$. Without loss of
generality, we can consider that there exists $\nu>1$ such that
$\mathbf{y}^{1}\leq\nu\mathbf{y}^{2}$ and $\exists$ at least one $j$ for which
$y^{(1)}_{j}<\nu y^{(2)}_{j}$ and $\exists$ at least one queue $k$ for which
$y^{(1)}_{k}=\nu y^{(2)}_{k}$. One can see easily that for any possible
vectors $\mathbf{y}^{1}=[y^{(1)}_{1},...,y^{(1)}_{K}]^{T}$ and
$\mathbf{y}^{2}=[y^{(2)}_{1},...,y^{(2)}_{n}]^{T}$ obtaining $\nu$ to satisfy
the above statement is straightforward. Recall that for queue $k$

$y^{(2)}_{k}=\frac{1}{\sum_{j^{\prime}=1}^{K}y^{(2)}_{j^{\prime}}}+\frac{\bar{\mu}_{k}}{\bar{\lambda}y^{(2)}_{k}}$

By using $\mathbf{y}^{1}\leq\nu\mathbf{y}^{2}$ and $y^{(1)}_{j}<\nu
y^{(2)}_{j}$ for at least one queue $j$, we get

$y^{(2)}_{k}<\frac{\nu}{\sum_{j^{\prime}=1}^{K}y^{(1)}_{j^{\prime}}}+\frac{\nu\bar{\mu}_{k}}{\bar{\lambda}y^{(1)}_{k}}=\nu
y^{(1)}_{k}$

Therefore, we obtain $y^{(2)}_{k}<\nu y^{(1)}_{k}$ which contradicts the fact
that $y^{(2)}_{k}=\nu y^{(1)}_{k}$. Consequently, it is not possible to have
two different solutions to the system of equations in (28). ∎

In order to solve the system of nonlinear equations in (28), we provide an
efficient iterative learning algorithm and prove its convergence to the
solution of (28). The proposed algorithm has a reduced complexity and can be
implemented in a distributed manner. We first consider the following iterative
algorithm:

$y_{j}(t+1)=(1-\alpha)y_{j}(t)+\alpha\frac{1}{\sum_{j^{\prime}=1}^{K}y_{j^{\prime}}(t)}+\alpha\frac{\bar{\mu}_{j}}{\bar{\lambda}y_{j}(t)},\
\forall j$ (29)

where $\alpha$ is a sufficiently small step size. This algorithm is a simple
class of Ishikawa algorithm (using Mann-like iteration but with constant step
size $\alpha$; see [29] for details).
By definition of $\mathbf{y}$, one can notice that $\forall j$
$y_{j}\in\mathcal{S}_{j}=\biggl{[}\sqrt{\frac{\bar{\mu}_{j}}{\bar{\lambda}}},\sqrt{1+\frac{\bar{\mu}_{j}}{\bar{\lambda}}}\biggr{]}$.
Let $\mathcal{S}=\mathcal{S}_{1}\times...\times\mathcal{S}_{K}$ be the set of
feasible values of $\mathbf{y}$. In order to ensure that the solution obtained
by the iterative algorithm lies in the feasible set $\mathcal{S}$, we consider
the following projection $\Pi_{\mathcal{S}}$:
$\hat{\mathbf{y}}=\Pi_{\mathcal{S}}(\mathbf{y})$ defined as,
$\forall j$, if $y_{j}<\sqrt{\frac{\bar{\mu}_{j}}{\bar{\lambda}}}$ then
$\hat{y}_{j}=\max{\\{y_{j},\sqrt{\frac{\bar{\mu}_{j}}{\bar{\lambda}}}\\}}$; if
$y_{j}>\sqrt{1+\frac{\bar{\mu}_{j}}{\bar{\lambda}}}$ then
$\hat{y}_{j}=\min{\\{y_{j},\sqrt{1+\frac{\bar{\mu}_{j}}{\bar{\lambda}}}\\}}$;
and $\hat{y}_{j}=y_{j}$ otherwise. Using the aforementioned projection, the
iterative algorithm becomes

$\hat{y}_{j}(t+1)=\Pi_{\mathcal{S}}\biggl{(}(1-\alpha)\hat{y}_{j}(t)+\alpha\frac{1}{\sum_{j^{\prime}=1}^{K}\hat{y}_{j^{\prime}}(t)}+\alpha\frac{\bar{\mu}_{j}}{\bar{\lambda}\hat{y}_{j}(t)}\biggr{)},\
\forall j$ (30)

The following result proves that the algorithm above converges to the solution
of the system of equations in (28), whenever it exists. Notice that once the
algorithm converges, one can obtain $\mathbf{m}$ from $\mathbf{y}$ from the
relation $y_{j}=\sqrt{m_{j}+\frac{\bar{\mu}_{j}}{\bar{\lambda}}}$ $\forall j$.
The routing probabilities for each source $i$, i.e. $\mathbf{p}_{i}$, can then
be obtained from (26). Finally, one can see that since the aforementioned
algorithm depends only on the average arrival and service rates, without
requiring any information about the instantaneous status of the network, it
can be implemented separately by each source.

######  Proposition 8.

The algorithm in (30) converges to the solution of the system of equations in
(28).

###### Proof.

We denote by
$f_{j}(\mathbf{y}(t))=\frac{1}{\sum_{j^{\prime}=1}^{K}y_{j^{\prime}}(t)}+\frac{\bar{\mu}_{j}}{\bar{\lambda}y_{j}(t)}$.
One can see that $f_{j}(\mathbf{y}(t))$ can be written as
$f_{j}(\mathbf{y}(t))=g(\mathbf{y}(t))+g_{j}(y_{j}(t))$ where
$\tilde{g}(\mathbf{y}(t))=\frac{1}{\sum_{j^{\prime}=1}^{K}y_{j^{\prime}}(t)}$
and $g_{j}(y_{j}(t))=\frac{\bar{\mu}_{j}}{\bar{\lambda}y_{j}(t)}$. We also
denote by $\mathbf{y}^{*}=[y^{*}_{1},...,y^{*}_{K}]\in\mathcal{S}$ the
solution of (28). Recall that by using the algorithm in (30), we have
$y_{j}(t+1))=(1-\alpha)\hat{y}_{j}(t)+\alpha\frac{1}{\sum_{j^{\prime}=1}^{K}\hat{y}_{j^{\prime}}(t)}+\alpha\frac{\bar{\mu}_{j}}{\bar{\lambda}\hat{y}_{j}(t)}$
and $\hat{y}_{j}(t+1)=\Pi_{\mathcal{S}}(y_{j}(t+1))$. We can also write the
following

$\displaystyle y_{j}(t+1)-y_{j}^{*}=(1-\alpha)(\hat{y}_{j}(t)-y_{j}^{*})+$
$\displaystyle\alpha(\tilde{g}(\hat{\mathbf{y}}(t))-\tilde{g}(\mathbf{y}^{*}))+\alpha_{t}(g_{j}(\hat{y}_{j}(t))-g_{j}(y_{j}^{*}))$
(31)

By using the mean value theorem, $\exists\
\mathbf{\bar{y}}(t)\in[\hat{\mathbf{y}}(t),\mathbf{y}^{*}]$ such that

$\tilde{g}(\hat{\mathbf{y}}(t))-\tilde{g}(\mathbf{y}^{*})=\nabla\tilde{g}|_{\mathbf{\bar{y}}(t)}^{T}\left(\hat{\mathbf{y}}(t)-\mathbf{y}^{*}\right)$

Similarly, $\exists\
\mathbf{\tilde{y}}(t)=[\tilde{y}_{1}(t),...,\tilde{y}_{K}(t)]^{T}\in[\hat{\mathbf{y}}(t),\mathbf{y}^{*}]$
such that

$g_{j}(\hat{y}_{j}(t))-g_{j}(y_{j}^{*})=\frac{dg_{j}}{dy_{j}}|_{\tilde{y}_{j}(t)}(\hat{y}_{j}(t)-y_{j}^{*}),\
\forall j$

From all the above, we can therefore write $\mathbf{y}(t+1)-\mathbf{y}^{*}$ as
follows:

$\displaystyle\mathbf{y}(t+1)-\mathbf{y}^{*}=(1-\alpha)(\hat{\mathbf{y}}(t)-\mathbf{y}^{*})+\alpha\mathbf{J}_{t}(\hat{\mathbf{y}}(t)-\mathbf{y}^{*})$
(32)

where $\mathbf{J}_{t}$ is a $K\times K$ matrix with diagonal elements
$-\left(\frac{1}{\sum_{l=1}^{K}\bar{y}_{j}(t)}\right)^{2}-\frac{\bar{\mu}_{j}}{\bar{\lambda}\tilde{y}^{2}_{j}(t)}$
and the non-diagonal elements are
$-\left(\frac{1}{\sum_{l=1}^{K}\bar{y}_{j}(t)}\right)^{2}$. We can show easily
then that $\forall\mathbf{a}\in\mathbbm{R}^{K}$
$\mathbf{a}^{T}\mathbf{J}_{t}\mathbf{a}\leq 0$, i.e. $\mathbf{J}_{t}$ is
negative semi definite. Then, by using
$\|\mathbf{y}(t+1)-\mathbf{y}^{*}\|^{2}=(\mathbf{y}(t+1)-\mathbf{y}^{*})^{T}(\mathbf{y}(t+1)-\mathbf{y}^{*})$,
we get

$\displaystyle\|\mathbf{y}(t+1)-\mathbf{y}^{*}\|^{2}=(1-\alpha)^{2}\|\hat{\mathbf{y}}(t)-\mathbf{y}^{*}\|^{2}+$
$\displaystyle\alpha^{2}\|\mathbf{J}_{t}(\hat{\mathbf{y}}(t)-\mathbf{y}^{*})\|^{2}+2(1-\alpha_{t})\alpha(\hat{\mathbf{y}}(t)-\mathbf{y}^{*})^{T}\mathbf{J}_{t}(\hat{\mathbf{y}}(t)-\mathbf{y}^{*})$
(33)

Since $\mathbf{J}$ is negative semi definite and
$\|\mathbf{J}_{t}(\hat{\mathbf{y}}(t)-\mathbf{y}^{*})\|^{2}\leq\|\mathbf{J}_{t}\|^{2}\|\hat{\mathbf{y}}(t)-\mathbf{y}^{*}\|^{2}$,
we get

$\displaystyle\|\mathbf{y}(t+1)-\mathbf{y}^{*}\|^{2}\leq\biggl{(}(1-\alpha)^{2}+\alpha^{2}\|\mathbf{J}_{t}\|^{2}\biggr{)}\|\hat{\mathbf{y}}(t)-\mathbf{y}^{*}\|^{2}$
(34)

Recall that $\hat{\mathbf{y}}\in\mathcal{S}$ and therefore
$\hat{\mathbf{y}}(t)\neq 0$, which implies that $\bar{\mathbf{y}}\neq 0$ and
$\tilde{\mathbf{y}}\neq 0$ and hence $\|\mathbf{J}_{t}\|^{2}$ is bounded.
Therefore, $\exists$ a sufficiently small $\alpha$ such that $\forall t$
$\alpha^{2}\|\mathbf{J}_{t}\|^{2}+\alpha^{2}-2\alpha<-\epsilon<0$, where
$0<\epsilon<1$. Consequently,

$\|\mathbf{y}(t+1)-\mathbf{y}^{*}\|^{2}\leq\biggl{(}1-\epsilon\biggr{)}\|\hat{\mathbf{y}}(t)-\mathbf{y}^{*}\|^{2}$

Then, by using the inequality
$\|\hat{\mathbf{y}}(t+1)-\mathbf{y}^{*}\|^{2}\leq\|\mathbf{y}(t+1)-\mathbf{y}^{*}\|^{2}$
(since $\mathbf{y}^{*}\in\mathcal{S}$), we get

$\displaystyle\|\hat{\mathbf{y}}(t+1)-\mathbf{y}^{*}\|^{2}$
$\displaystyle\leq\biggl{(}1-\epsilon\biggr{)}\|\hat{\mathbf{y}}(t)-\mathbf{y}^{*}\|^{2}$
$\displaystyle\mathmakebox[\widthof{{}\leq{}}][c]{\vdots}$
$\displaystyle\leq\biggl{(}1-\epsilon\biggr{)}^{t+1}\|\hat{\mathbf{y}}(0)-\mathbf{y}^{*}\|^{2}$
(35)

and $\|\hat{\mathbf{y}}(t)-\mathbf{y}^{*}\|^{2}\rightarrow 0$ when
$t\rightarrow\infty$. This concludes the proof. ∎

From the Proof above, one can see that if $\alpha$ is taken such that
$\alpha^{2}\|\mathbf{J}_{t}\|^{2}+\alpha^{2}-2\alpha<0$ then the algorithm in
(30) converges to the solution of (28). Since for each $j$
$\hat{y}_{j}\geq\sqrt{\frac{\bar{\mu}_{j}}{\bar{\lambda}}}$ and
$\|\mathbf{J}_{t}\|^{2}\leq\sum_{i=1}^{K}\sum_{j=1}^{K}J_{t}(i,j)$, which
implies that $\|\mathbf{J}_{t}\|^{2}\leq
K^{2}\left(\frac{1}{\sum_{l=1}^{K}\sqrt{\frac{\bar{\mu}_{j}}{\bar{\lambda}}}}\right)^{2}+K$.
Consequently, it is sufficient to take
$\alpha<\frac{2}{K^{2}\left(\frac{1}{\sum_{l=1}^{K}\sqrt{\frac{\bar{\mu}_{j}}{\bar{\lambda}}}}\right)^{2}+K+1}$.

####  IV-D2 Numerical Results

We provide here an example to show numerically that the system of equations in
(25) has a solution. While in the previous subsection, a detailed mathematical
analysis is provided to optimize the upper bound in the case of large number
of sources, we show here numerical results for finite number of sources. In
order to solve (25), we consider the following iterative algorithm:
$p_{ij}(t+1)=(1-\alpha)p_{ij}(t)+\alpha\frac{\sqrt{\sum_{l\neq
i}^{N}\lambda_{l}p_{lj}(t)+\mu_{j}}}{\sum_{j^{\prime}=1}^{K}\sqrt{\sum_{l\neq
i}^{N}\lambda_{l}p_{lj^{\prime}}(t)+\mu_{j^{\prime}}}},\ \forall i,j$. One can
see that actually this iterative algorithm is similar to (30), or more
precisely (30) can be obtained from the above algorithm by using Mean Field
analysis (using the derivations in the previous subsection).  In Figure 21 we
show an illustrative example where we consider a system with 6 sources et 10
servers with the following parameters for the sources $\lambda_{1}=100$,
$\lambda_{2}=20$, $\lambda_{3}=50$, $\lambda_{4}=\lambda_{5}=10$ and
$\lambda_{6}=1000$ and for the servers $\mu_{1}=1,$ $\mu_{2}=2,$ $\mu_{3}=3$,
$\mu_{4}=5$, $\mu_{5}=10$, $\mu_{6}=20$, $\mu_{7}=50$, $\mu_{8}=100$,
$\mu_{9}=200$ and $\mu_{10}=1000$. The plots show that (25) has a solution,
that $p_{1,2}$, $p_{1,6}$, $p_{1,7}$, $p_{1,8}$, $p_{1,9}$ and $p_{1,10}$
converge to the solution of (25) and also that in this example the algorithm
converges quickly. This shows that the proposed algorithm can converge also
for finite number of sources with different arrival rates. Providing a formal
proof of convergence would be an interesting topic and is left for future
work.

Figure 21: Convergence of the fixed point algorithm defined from (25).

##  V Conclusion

In this paper, we studied the average AoI of a system of multiple sources and
parallel queues using the SHS method. We considered that the queues do not
communicate between each other and that the sources send their updates to the
queues according to a predefined probabilistic routing scheme. First, we
computed the average AoI for the following systems: i) two parallel M/M/1/1
queues, ii) one M/M/1/1 queue with half arrival and loss rates, and iii) one
M/M/1/1 queue with double service rate. Then, we computed the average AoI for
two parallel M/M/1/2* queues, one M/M/1/3* queue with half arrival and loss
rates, and one M/M/1/3* queue with double service time. We conclude that the
average AoI of the system composed of parallel queues is always smaller than
that of one queue with half arrival and loss rates, and can be as small as
that of one queue with double service rate. We also studied the average AoI of
a system with an arbitrary number of heterogeneous M/M/1/(N+1)* queues and we
provided an upper bound of AoI that is tight when there are multiple sources.
We then provided a framework allowing each source to determine its routing
decision, by using Game Theory and best response method. In the contest of
large number of sources, we simplified the game framework by using Mean Field
Games, provided a simple distributed algorithm and proved its convergence to
the desired fixed point.

## References

  * [1] S. Kaul, M. Gruteser, V. Rai, and J. Kenney, “Minimizing age of information in vehicular networks,” in _8th IEEE Annual Conference on Sensor, Mesh and Ad Hoc Communications and Networks_ , 2011, pp. 350–358.
  * [2] S. Kaul, R. Yates, and M. Gruteser, “Real-time status: How often should one update?” in _Proc. of IEEE INFOCOM_ , 2012, pp. 2731–2735.
  * [3] C. Kam, S. Kompella, and A. Ephremides, “Effect of message transmission diversity on status age,” in _IEEE International Symposium on Information Theory_ , 2014, pp. 2411–2415.
  * [4] C. Kam, S. Kompella, G. D. Nguyen, and A. Ephremides, “Effect of message transmission path diversity on status age,” _IEEE Transactions on Information Theory_ , vol. 62, no. 3, pp. 1360–1374, 2015.
  * [5] S. K. Kaul and R. D. Yates, “Timely updates by multiple sources: The m/m/1 queue revisited,” in _IEEE 54th Annual Conference on Information Sciences and Systems (CISS)_ , 2020, pp. 1–6.
  * [6] M. Moltafet, M. Leinonen, and M. Codreanu, “On the age of information in multi-source queueing models,” _IEEE Transactions on Communications, Early Access_ , 2020.
  * [7] A. Kosta, N. Pappas, V. Angelakis _et al._ , “Age of information: A new concept, metric, and tool,” _Foundations and Trends® in Networking_ , vol. 12, no. 3, pp. 162–259, 2017.
  * [8] A. Maatouk, M. Assaad, and A. Ephremides, “The age of updates in a simple relay network,” in _Proc. of IEEE Information Theory Workshop (ITW)_ , 2018\.
  * [9] Y.Sun, E.Uysal-Biyikoglu, and S.Kompella, “Age-optimal updates of multiple information flows,” _arXiv preprint, arXiv:1801.02394_ , 2018.
  * [10] R. Talak, S. Karaman, and E. Modiano, “Can determinacy minimize age of information?” _arXiv preprint arXiv:1810.04371_ , 2018.
  * [11] Y. Sun, E. Uysal-Biyikoglu, R. D. Yates, C. E. Koksal, and N. B. Shroff, “Update or wait: How to keep your data fresh,” _IEEE Transactions on Information Theory_ , vol. 63, no. 11, pp. 7492–7508, 2017.
  * [12] A. Maatouk, S. Kriouile, M. Assaad, and A. Ephremides, “The age of incorrect information: A new performance metric for status updates,” _To appear in IEEE/ACM Transactions on Networking, arXiv preprint arXiv:1907.06604_ , 2020\.
  * [13] Y.-P. Hsu, E. Modiano, and L. Duan, “Scheduling algorithms for minimizing age of information in wireless broadcast networks with random arrivals: The no-buffer case,” _arXiv e-prints, p. arXiv:1712.07419_ , 2017.
  * [14] A. M. Bedewy, Y. Sun, S. Kompella, and N. B. Shroff, “Age-optimal sampling and transmission scheduling in multi-source systems,” _arXiv e-prints, p. arXiv:1812.09463_ , 2018.
  * [15] I. Kadota, A. Sinha, E. Uysal-Biyikoglu, R. Singh, and E. Modiano, “Scheduling policies for minimizing age of information in broadcast wireless networks,” _IEEE/ACM Transactions on Networking_ , vol. 26, no. 6, pp. 2637–2650, 2018\.
  * [16] A. Maatouk, M. Assaad, and A. Ephremides, “Minimizing the age of information: NOMA or OMA?” in _Proc. of Infocom AoI Workshop 2019, arXiv e-prints, p. arXiv:1901.03020_ , 2019.
  * [17] ——, “Minimizing the age of information in a CSMA environment,” in _Proc. of IEEE WiOPT, arXiv e-prints, p. arXiv:1901.00481_ , 2019.
  * [18] ——, “On the age of information in a CSMA environment,” _To appear in IEEE Transactions on Networking_ , 2020.
  * [19] A. M. Bedewy, Y. Sun, and N. B. Shroff, “Optimizing data freshness, throughput, and delay in multi-server information-update systems,” in _IEEE International Symposium on Information Theory (ISIT)_ , 2016, pp. 2569–2574.
  * [20] Y. Sun, E. Uysal-Biyikoglu, and S. Kompella, “Age-optimal updates of multiple information flows,” in _IEEE Conference on Computer Communications (INFOCOM) Workshops_ , April 2018, pp. 2281–2285.
  * [21] R. D. Yates, “Status updates through networks of parallel servers,” in _IEEE International Symposium on Information Theory (ISIT)_ , June 2018, pp. 2281–2285.
  * [22] A. M. Bedewy, Y. Sun, and N. B. Shroff, “The age of information in multihop networks,” _IEEE/ACM Transactions on Networking_ , vol. 27, no. 3, pp. 1248–1257, 2019.
  * [23] ——, “Age-optimal information updates in multihop networks,” in _IEEE International Symposium on Information Theory (ISIT)_ , 2017, pp. 576–580.
  * [24] S. Farazi, A. G. Klein, , and D. R. Brown, “Fundamental bounds on the age of information in general multi-hop interference networks,” in _IEEE Conference on Computer Communications (INFOCOM) Workshops_ , April 2019, pp. 96–101.
  * [25] A. S. B. Buyukates and S. Ulukus, “Age of information in multihop multicast networks,” _Journal of Communications and Networks_ , vol. 21, no. 3, pp. 256–267, 2019.
  * [26] R. D. Yates, Y. Sun, D. R. B. III, S. K. Kaul, E. Modiano, and S. Ulukus, “Age of information: An introduction and survey,” _https://arxiv.org/pdf/2007.08564.pdf_ , 2020.
  * [27] R. D. Yates and S. K. Kaul, “The age of information: Real-time status updating by multiple sources,” _IEEE Transactions on Information Theory_ , vol. 65, no. 3, pp. 1807–1827, 2019.
  * [28] J. P. Hespanha, “Modelling and analysis of stochastic hybrid systems,” _IEE Proceedings-Control Theory and Applications_ , vol. 153, no. 5, pp. 520–535, 2006.
  * [29] S. Ishikawa, “Fixed points by a new iteration method,” in _Proc. Amer. Math. Soc._ , vol. 44, 1974, pp. 147–150.

##  Appendix A Proof of Theorem 20

$l$ | $\lambda^{l}$ | x | $x^{\prime}=xA_{l}$ | $\bar{v}_{q_{l}}A_{l}$
---|---|---|---|---
0 | $\lambda_{1}p_{1j}$ | $[x_{0}\ \dots\ x_{j1}\dots\ x_{jN}\ x_{jM}\ \dots]$ | $[x_{0}\ \dots\ x_{j1}\dots\ x_{jN}\ 0\ \dots]$ | $[v_{0}\ \dots\ v_{j1}\dots\ v_{jN}\ 0\ \dots]$
1 | $\sum_{k>1}\lambda_{k}p_{kj}$ | $[x_{0}\ \dots\ x_{j1}\dots\ x_{jN}\ x_{jM}\ \dots]$ | $[x_{0}\ \dots\ x_{j1}\dots\ x_{jN}\ x_{jN}\ \dots]$ | $[v_{0}\ \dots\ v_{j1}\dots\ v_{jN}\ v_{jN}\ \dots]$
2 | $\mu_{j}$ | $[x_{0}\ \dots\ x_{j1}\dots\ x_{jN}\ x_{jM}\ \dots]$ | $[x_{j1}\ \dots\ x_{j2}\dots\ x_{jM}\ x_{jM}\ \dots]$ | $[v_{0}\ \dots\ v_{j1}\dots\ v_{jM}\ v_{jM}\ \dots]$
TABLE VI: Table of SHS transitions of Figure 22.

We consider a system with $K$ parallel M/M/1/(N+1)* queues. Without loss of
generality, we compute the average AoI of source 1. The result is of course
true for any source $i$. The arrival rate to queue $j$ from source $1$ is
$\lambda_{1}p_{1j}$ and that from the rest of the sources is
$\sum_{k>1}\lambda_{k}p_{kj}$ for all $j=1,\dots,K$. Let $M=N+1.$ The
continuous state is given by a vector of size $1+K*M$

$\mathbf{x}(t)=[x_{0}(t)\ x_{11}(t)\ \dots\ x_{1M}(t)\ x_{21}(t)\ \ x_{K1}(t)\
\dots\ \\\ x_{KM}(t)\ ],$

where $x_{0}$ represents the current age, $x_{j1}$ the age if the update in
service in the queue $j$ is delivered and $x_{jl}$ the age if the update in
the position $l-1$ of the queue $j$ is delivered. The discrete state is a
Markov Chain with a single state. We note that, when an event (arrival or
departure) occurs in queue $j$, the age of the updates in the rest of the
queues is not modified. This allows us to focus on a queue $j$ to illustrate
the Markov Chain and the SHS transitions, which are presented respectively in
Figure 22 and Table VI.

$0$$0$$2$$1$ Figure 22: The SHS Markov Chain for a system formed by $K$
parallel M/M/1/(N+1)* queues with multiple sources.

We now explain each transition $l$:

  *     * $l=0$

There is an update of source $i$ that arrives to queue $j$. For this case, the
incoming update replaces the update in the last position of the queue and the
age of the incoming update is set to zero, i.e., $x_{jM}^{\prime}=0$.

    * $l=1$

There is an update of another source that arrives to queue $j$. For this case,
the incoming update replaces the update in the last position of the queue and
the age of the incoming update is set to the same value as the age of
penultimate update, , i.e., $x_{jM}^{\prime}=x_{jN}$.

######  Remark 3.

If $N=0$, the last update in the queue is an update that is in service.
Therefore, when an update of other sources arrives to the system, the incoming
update replaces the update in service and $x_{j1}^{\prime}=x_{0}$.

    * $l=2$

The update in service in queue $j$ is delivered and therefore the age of the
monitor changes to $x_{j1}$. Besides, all the elements in the queue move a
position ahead, which causes that their ages change respectively from $x_{j1}$
to $x_{j2}$, from $x_{j2}$ to $x_{j3}$, $\dots$ and from $x_{jN}$ to $x_{jM}$.
Finally, in the last position of the queue, we put a fake update whose age
value is set to $x_{jM}$, that is, the age of the penultimate element in the
queue, i.e., $x_{jM}^{\prime}=x_{jM}$.

As we have just mentioned, when an update of queue $i$ is delivered to the
monitor, a fake update is put in the last position of the queue. We now
explain that these additional fake updates lead to a larger sojourn time of
the incoming updates, which implies that the overall average AoI is larger
and, as a result, the result we obtain provides an upper bound of the real
average AoI of the system (i.e. without fake updates). Consider that all the
updates of queue $i$ are delivered before the arrival of a new update to that
queue. According to the above explained system, the queue is full of fake
updates and, therefore, when a new update arrives to the system, it is
enqueued in the last position of the queue. However, in a system without fake
updates, a new update would find the queue empty and would start the service
upon arrival. As a result, we have that the sojourn time of new updates in a
system with fake updates (i.e., the above presented system) is clearly larger
than the sojourn time without fake updates.

Since the Markov Chain is formed by a single state, the stationary
distribution is trivial. We define the vector $\mathbf{v}=[v_{0}\ v_{11}\
v_{12}\ v_{1M}\ v_{21}\ \dots\ v_{K1}\dots v_{KM}]$ and also $\mathbf{b}$ as
the vector of size $1+K*M$ with all ones. From the result of Theorem 4 in [27]
and the above reasoning, we know that an upper bound of the AoI is given by
$v_{0}$, that is, the first coordinate of the vector $\mathbf{v}$.

In the remainder of the proof, we present the system of equations that
$\mathbf{v}$ satisfies and solve it. We first present the equation of the
first coordinate of $\mathbf{v}$:

$\sum_{j=1}^{K}(\lambda_{1}p_{1j}+\sum_{k>1}\lambda_{k}p_{kj}+\mu_{j})v_{0}=\\\
1+\sum_{j=1}^{K}\left((\lambda_{1}p_{1j}+\sum_{k>1}\lambda_{k}p_{kj})v_{0}+\mu_{j}v_{j1}\right),$

which can be alternatively written as

$v_{0}\sum_{j=1}^{K}\mu_{j}=1+\sum_{j=1}^{K}\mu_{j}v_{j1}.$ (36)

Let $l^{\prime}=l+1.$We now present that, for all $j=1,\dots,K$ and all
$l=1,\dots,M$, the following equation is satisfied:

$\displaystyle\sum_{m=1}^{K}(\lambda_{1}p_{1m}+\sum_{k>1}\lambda_{k}p_{km}+\mu_{m})v_{ml}=1+$
$\displaystyle\sum_{m\neq
j}(\lambda_{1}p_{1m}+\sum_{k>1}\lambda_{k}p_{km}+\mu_{m})v_{ml}+$
$\displaystyle(\lambda_{1}p_{1j}+\sum_{k>1}\lambda_{k}p_{kj})v_{jl}+\mu_{j}v_{ml^{\prime}}$
$\displaystyle\iff\mu_{j}v_{jl}=1+\mu_{j}v_{jl^{\prime}}$ (37)

Using recursively the last expression for $l$ equals $1$ to $N$, we get that

$\mu_{j}v_{j1}=N-1+\mu_{j}v_{mN}.$ (38)

We now focus on the last position of queue $j$ and the equation that it must
satisfy is the following:

$\displaystyle\sum_{m=1}^{K}(\lambda_{1}p_{1m}+\sum_{k>1}\lambda_{k}p_{km}+\mu_{m})v_{mM}=1+$
$\displaystyle\sum_{m\neq
j}(\lambda_{1}p_{1m}+\sum_{k>1}\lambda_{k}p_{km}+\mu_{m})v_{mM}+\sum_{k>1}\lambda_{k}p_{kj}v_{jN}+\mu_{j}v_{jM}$
$\displaystyle\iff(\lambda_{1}p_{1j}+\sum_{k>1}\lambda_{k}p_{kj})v_{jM}=1+\sum_{k>1}\lambda_{k}p_{kj}v_{jN}$

Besides, from (37), for l=N, we have that

$\mu_{j}v_{jN}=1+\mu_{j}v_{jM}.$

Using the last two expressions, we get that

$\displaystyle\mu_{j}v_{jN}$ $\displaystyle=1+\mu_{j}v_{jM}$
$\displaystyle=1+\mu_{j}\frac{1+\sum_{k>1}\lambda_{k}p_{kj}v_{jN}}{\lambda_{1}p_{1j}+\sum_{k>1}\lambda_{k}p_{kj}}.$

The last expression is equivalent to the following one:

$\mu_{j}v_{jN}\left(1-\frac{\lambda_{2}p_{2j}}{\lambda_{1}p_{1j}+\sum_{k>1}\lambda_{k}p_{kj}}\right)=1+\frac{\mu_{j}}{\lambda_{1}p_{1j}+\sum_{k>1}\lambda_{k}p_{kj}}\iff$
$\mu_{j}v_{jN}\left(\frac{\lambda_{1p_{1j}}}{\lambda_{1p_{1j}}+\sum_{k>1}\lambda_{k}p_{kj}}\right)=\frac{\lambda_{1}p_{1j}+\sum_{k>1}\lambda_{k}p_{kj}+\mu_{j}}{\lambda_{1}p_{1j}+\sum_{k>1}\lambda_{k}p_{kj}}\iff$
$\mu_{j}v_{jN}\lambda_{1}p_{1j}=\lambda_{1}p_{1j}+\sum_{k>1}\lambda_{k}p_{kj}+\mu_{j}\iff$
$\mu_{j}v_{jN}=1+\frac{\sum_{k>1}\lambda_{k}p_{kj}+\mu_{j}}{\lambda_{1}p_{1j}}.$

Using the last expression with (38) and (36), the desired result follows for
$i=1$.

Figure 23:  Simulation results of the M/M/1/1 queue with different values of
arrival rates; x-axis in logarithmic scale. Figure 24: Simulation results of
the M/M/1/3* queue with different values of arrival rates; x-axis in
logarithmic scale.

##  Appendix B Simulations for M/M/1/1 and M/M/1/3* queues

We provide in this section a comparison between the results obtained by SHS
and those obtained by simulation to show the accuracy of the SHS method. We
consider two different examples with either M/M/1/1 or M/M/1/3* queues. We
consider a system with a single source and without packet losses. We set
$\mu=2$. As expected, the curves presented in Figures 23 and 24 show that, in
all the considered cases, the results obtained by simulation coincide with
those of SHS.

##  Appendix C Table of transitions of Figure 9

$l$ | $q_{l}\to q_{l^{\prime}}$ | $\lambda^{l}$ | $x^{\prime}=xA_{l}$ | $\bar{v}_{q_{l}}A_{l}$
---|---|---|---|---
0 | $00\to 01$ | $\lambda_{1}p_{12}$ | $[x_{0}\ 0\ 0\ 0\ 0]$ | $[v_{00}(0)\ 0\ 0\ 0\ 0]$
|  | $\sum_{k>1}\lambda_{k}p_{k2}$ | $[x_{0}\ 0\ 0\ x_{0}\ 0]$ | $[v_{00}(0)\ 0\ 0\ v_{00}(0)\ 0]$
1 | $01\to 00$ | $\mu_{2}$ | $[x_{21}\ 0\ 0\ 0\ 0]$ | $[v_{01}(3)\ 0\ 0\ 0\ 0]$
|  | $\theta_{2}$ | $[x_{0}\ 0\ 0\ 0\ 0]$ | $[v_{01}(0)\ 0\ 0\ 0\ 0]$
2 | $01\to 02$ | $\lambda_{1}p_{12}$ | $[x_{0}\ 0\ 0\ x_{21}\ 0]$ | $[v_{01}(0)\ 0\ 0\ 0\ 0]$
|  | $\sum_{k>1}\lambda_{k}p_{k2}$ | $[x_{0}\ 0\ 0\ x_{21}\ x_{21}]$ | $[v_{01}(0)\ 0\ 0\ v_{00}(0)\ 0]$
3 | $02\to 01$ | $\mu_{2}$ | $[x_{21}\ 0\ 0\ x_{22}\ 0]$ | $[v_{02}(3)\ 0\ 0\ v_{02}(4)\ 0]$
|  | $\theta_{2}$ | $[x_{0}\ 0\ 0\ x_{22}\ 0]$ | $[v_{02}(0)\ 0\ 0\ v_{02}(4)\ 0]$
4 | $02\to 02$ | $\lambda_{1}p_{12}$ | $[x_{0}\ 0\ 0\ x_{21}\ 0]$ | $[v_{02}(0)\ 0\ 0\ v_{02}(3)\ 0]$
|  | $\sum_{k>1}\lambda_{k}p_{k2}$ | $[x_{0}\ 0\ 0\ x_{21}\ x_{21}]$ | $[v_{02}(0)\ 0\ 0\ v_{02}(3)\ v_{02}(3)]$
5 | $00\to 10$ | $\lambda_{1}p_{11}$ | $[x_{0}\ 0\ 0\ 0\ 0]$ | $[v_{00}(0)\ 0\ 0\ 0\ 0]$
|  | $\sum_{k>1}\lambda_{k}p_{k1}$ | $[x_{0}\ x_{0}\ 0\ 0\ 0]$ | $[v_{00}(0)\ v_{00}(0)\ 0\ 0\ 0\ 0]$
6 | $10\to 00$ | $\mu_{1}$ | $[x_{12}\ 0\ 0\ 0\ 0]$ | $[v_{01}(1)\ 0\ 0\ 0\ 0]$
|  | $\theta_{1}$ | $[x_{0}\ 0\ 0\ 0\ 0]$ | $[v_{01}(0)\ 0\ 0\ 0\ 0]$
7 | $01\to 11$ | $\lambda_{1}p_{11}$ | $[x_{0}\ 0\ 0\ x_{21}\ 0]$ | $[v_{01}(0)\ 0\ 0\ v_{01}(3)\ 0]$
|  | $\sum_{k>1}\lambda_{k}p_{k1}$ | $[x_{0}\ x_{0}\ 0\ x_{21}\ 0]$ | $[v_{01}(0)\ v_{01}(0)\ 0\ 0\ v_{01}(3)\ 0]$
8 | $11\to 01$ | $\mu_{1}$ | $[x_{11}\ 0\ 0\ x_{21}\ 0]$ | $[v_{11}(1)\ 0\ 0\ v_{11}(3)\ 0]$
|  | $\theta_{1}$ | $[x_{0}\ 0\ 0\ x_{21}\ 0]$ | $[v_{11}(0)\ 0\ 0\ v_{11}(3)\ 0]$
9 | $02\to 12$ | $\lambda_{1}p_{11}$ | $[x_{0}\ 0\ 0\ x_{21}\ x_{22}]$ | $[v_{02}(0)\ 0\ 0\ v_{02}(3)\ v_{02}(4)]$
|  | $\sum_{k>1}\lambda_{k}p_{k1}$ | $[x_{0}\ x_{0}\ 0\ x_{21}\ x_{22}]$ | $[v_{02}(0)\ v_{02}(0)\ 0\ 0\ v_{02}(3)\ v_{02}(4)]$
10 | $12\to 02$ | $\mu_{1}$ | $[x_{11}\ 0\ 0\ x_{21}\ x_{22}]$ | $[v_{12}(1)\ 0\ 0\ v_{12}(3)\ v_{12}(4)]$
|  | $\theta_{1}$ | $[x_{0}\ 0\ 0\ x_{21}\ x_{22}]$ | $[v_{12}(0)\ 0\ 0\ v_{12}(3)\ v_{12}(4)]$
11 | $10\to 11$ | $\lambda_{1}p_{12}$ | $[x_{0}\ x_{11}\ 0\ 0\ 0]$ | $[v_{10}(0)\ v_{10}(1)\ 0\ 0\ 0]$
|  | $\sum_{k>1}\lambda_{k}p_{k2}$ | $[x_{0}\ x_{11}\ 0\ x_{0}\ 0]$ | $[v_{10}(0)\ v_{10}(1)\ 0\ v_{10}(0)\ 0]$
12 | $11\to 10$ | $\mu_{2}$ | $[x_{21}\ x_{11}\ 0\ 0\ 0]$ | $[v_{11}(3)\ v_{11}(1)\ 0\ 0\ 0]$
|  | $\theta_{2}$ | $[x_{0}\ x_{11}\ 0\ 0\ 0]$ | $[v_{11}(0)\ v_{11}(1)\ 0\ 0\ 0]$
13 | $11\to 12$ | $\lambda_{1}p_{12}$ | $[x_{0}\ x_{11}\ 0\ x_{21}\ 0]$ | $[v_{11}(0)\ v_{11}(1)\ 0\ v_{11}(3)\ 0]$
|  | $\sum_{k>1}\lambda_{k}p_{k2}$ | $[x_{0}\ x_{11}\ 0\ x_{21}\ x_{21}]$ | $[v_{11}(0)\ v_{11}(1)\ 0\ v_{11}(3)\ v_{11}(3)]$
14 | $12\to 11$ | $\mu_{2}$ | $[x_{21}\ 0\ 0\ x_{22}\ 0]$ | $[v_{12}(3)\ 0\ 0\ v_{12}(4)\ 0]$
|  | $\theta_{2}$ | $[x_{0}\ 0\ 0\ x_{22}\ 0]$ | $[v_{12}(0)\ 0\ 0\ v_{12}(4)\ 0]$
15 | $12\to 12$ | $\lambda_{1}p_{12}$ | $[x_{0}\ x_{11}\ 0\ x_{21}\ 0]$ | $[v_{12}(0)\ v_{12}(1)\ 0\ v_{12}(3)\ 0]$
|  | $\sum_{k>1}\lambda_{k}p_{k2}$ | $[x_{0}\ x_{11}\ 0\ x_{21}\ x_{21}]$ | $[v_{12}(0)\ v_{12}(1)\ 0\ v_{12}(3)\ v_{12}(3)]$
16 | $10\to 20$ | $\lambda_{1}p_{11}$ | $[x_{0}\ x_{11}\ 0\ 0\ 0]$ | $[v_{10}(0)\ x_{10}(1)\ 0\ 0\ 0]$
|  | $\sum_{k>1}\lambda_{k}p_{k1}$ | $[x_{0}\ x_{11}\ x_{11}\ 0\ 0]$ | $[v_{10}(0)\ v_{10}(1)\ v_{10}(1)\ 0\ 0\ 0]$
17 | $20\to 10$ | $\mu_{1}$ | $[x_{12}\ x_{22}\ 0\ 0\ 0]$ | $[v_{20}(1)\ v_{20}(2)\ 0\ 0\ 0]$
|  | $\theta_{1}$ | $[x_{0}\ x_{22}\ 0\ 0\ 0]$ | $[v_{20}(0)\ x_{20}(2)\ 0\ 0\ 0]$
18 | $11\to 21$ | $\lambda_{1}p_{11}$ | $[x_{0}\ x_{11}\ 0\ x_{21}\ 0]$ | $[v_{11}(0)\ v_{11}(1)\ 0\ v_{11}(3)\ 0]$
|  | $\sum_{k>1}\lambda_{k}p_{k1}$ | $[x_{0}\ x_{11}\ x_{11}\ x_{21}\ 0]$ | $[v_{11}(0)\ v_{11}(1)\ v_{11}(1)\ v_{11}(3)\ 0]$
19 | $21\to 11$ | $\mu_{1}$ | $[x_{11}\ x_{12}\ 0\ x_{21}\ 0]$ | $[v_{21}(1)\ v_{21}(2)\ 0\ v_{21}(3)\ 0]$
|  | $\theta_{1}$ | $[x_{0}\ x_{12}\ 0\ x_{21}\ 0]$ | $[v_{21}(0)\ v_{21}(2)\ 0\ v_{11}(3)\ 0]$
20 | $12\to 22$ | $\lambda_{1}p_{11}$ | $[x_{0}\ x_{11}\ 0\ x_{21}\ x_{22}]$ | $[v_{12}(0)\ v_{12}(1)\ 0\ v_{12}(3)\ v_{12}(4)]$
|  | $\sum_{k>1}\lambda_{k}p_{k1}$ | $[x_{0}\ x_{11}\ 0\ x_{21}\ x_{22}]$ | $[v_{12}(0)\ v_{12}(1)\ v_{12}(1)\ v_{12}(3)\ v_{12}(4)]$
21 | $22\to 12$ | $\mu_{1}$ | $[x_{11}\ x_{12}\ 0\ x_{21}\ x_{22}]$ | $[v_{22}(1)\ v_{22}(2)\ 0\ v_{22}(3)\ v_{12}(4)]$
|  | $\theta_{1}$ | $[x_{0}\ 0\ 0\ x_{21}\ x_{22}]$ | $[v_{22}(0)\ v_{22}(2)\ 0\ v_{22}(3)\ v_{22}(4)]$
22 | $20\to 20$ | $\lambda_{1}p_{11}$ | $[x_{0}\ x_{11}\ 0\ 0\ 0]$ | $[v_{20}(0)\ v_{20}(1)\ 0\ 0\ 0]$
|  | $\sum_{k>1}\lambda_{k}p_{k1}$ | $[x_{0}\ x_{11}\ x_{11}\ x_{21}\ 0]$ | $[v_{22}(0)\ v_{20}(1)\ v_{20}(1)\ 0\ 0\ 0]$
23 | $20\to 21$ | $\lambda_{1}p_{12}$ | $[x_{0}\ x_{11}\ x_{12}\ 0\ 0]$ | $[v_{20}(0)\ v_{20}(1)\ v_{20}(2)\ 0\ 0]$
|  | $\sum_{k>1}\lambda_{k}p_{k2}$ | $[x_{0}\ x_{11}\ x_{12}\ x_{0}\ 0]$ | $[v_{20}(0)\ v_{20}(1)\ v_{20}(2)\ v_{20}(0)\ 0]$
24 | $21\to 20$ | $\mu_{2}$ | $[x_{21}\ x_{11}\ x_{12}\ 0\ 0]$ | $[v_{21}(3)\ v_{21}(1)\ v_{21}(2)\ 0\ 0]$
|  | $\theta_{2}$ | $[x_{0}\ x_{11}\ x_{12}\ 0\ 0]$ | $[v_{21}(0)\ v_{21}(1)\ v_{21}(2)\ 0\ 0]$
25 | $21\to 21$ | $\lambda_{1}p_{11}$ | $[x_{0}\ x_{11}\ 0\ x_{21}\ 0]$ | $[v_{21}(0)\ v_{21}(1)\ 0\ v_{21}(3)\ 0]$
|  | $\sum_{k>1}\lambda_{k}p_{k1}$ | $[x_{0}\ x_{11}\ x_{11}\ x_{21}\ 0]$ | $[v_{21}(0)\ v_{21}(1)\ v_{21}(1)\ v_{21}(3)\ 0]$
26 | $21\to 22$ | $\lambda_{1}p_{12}$ | $[x_{0}\ x_{11}\ x_{12}\ x_{21}\ 0]$ | $[v_{21}(0)\ v_{21}(1)\ v_{21}(2)\ x_{21}(3)\ 0]$
|  | $\sum_{k>1}\lambda_{k}p_{k2}$ | $[x_{0}\ x_{11}\ x_{12}\ x_{21}\ x_{21}]$ | $[v_{21}(0)\ v_{21}(1)\ v_{21}(2)\ v_{21}(3)\ v_{21}(3)]$
27 | $22\to 21$ | $\mu_{2}$ | $[x_{21}\ x_{11}\ x_{12}\ x_{21}\ 0]$ | $[v_{22}(3)\ v_{22}(1)\ v_{22}(2)\ x_{22}(4)\ 0]$
|  | $\theta_{2}$ | $[x_{0}\ x_{11}\ x_{12}\ x_{21}\ 0]$ | $[v_{22}(0)\ v_{22}(1)\ v_{22}(2)\ v_{22}(4)\ 0]$
28 | $22\to 22$ | $\lambda_{1}p_{11}$ | $[x_{0}\ x_{11}\ 0\ x_{21}\ x_{22}]$ | $[v_{22}(0)\ v_{22}(1)\ 0\ v_{22}(3)\ v_{22}(4)]$
|  | $\sum_{k>1}\lambda_{k}p_{k1}$ | $[x_{0}\ x_{11}\ x_{11}\ x_{21}\ x_{22}]$ | $[v_{22}(0)\ v_{22}(1)\ v_{22}(1)\ v_{22}(3)\ v_{22}(4)]$
29 | $22\to 22$ | $\lambda_{1}p_{12}$ | $[x_{0}\ x_{11}\ x_{12}\ x_{21}\ 0]$ | $[v_{22}(0)\ v_{22}(1)\ v_{22}(2)\ v_{22}(3)\ 0]$
|  | $\sum_{k>1}\lambda_{k}p_{k2}$ | $[x_{0}\ x_{11}\ x_{12}\ x_{21}\ x_{21}]$ | $[v_{22}(0)\ v_{22}(1)\ v_{22}(2)\ v_{22}(3)\ v_{22}(3)]$
TABLE VII: Table of transitions of Figure 9.