Characterizing the Age of Information with Multiple Coexisting Data Streams
Yoshiaki Inoue Michel Mandjes




op-tical net-works semi-conduc-tor










Characterizing the Age of Information with 
Multiple Coexisting Data Streams

    
Yoshiaki Inoue and Michel Mandjes
This work was supported in part by JSPS KAKENHI Grant Number
21H03399 and the NWO Gravitation project NETWORKS under grant agreement no. 024.002.003.

Y. Inoue is with 
Department of Information and Communications Technology, 
Graduate School of Engineering, Osaka University, Suita 565-0871, Japan 
(e-mail: yoshiaki@comm.eng.osaka-u.ac.jp).

M. Mandjes is with Mathematical Institute, Leiden University, P.O. Box 9512, 2300 RA Leiden, 
The Netherlands (e-mail: m.r.h.mandjes@math.leidenuniv.nl). He is also affiliated with Korteweg-de Vries Institute for Mathematics, University of Amsterdam, Amsterdam, The Netherlands; Eurandom, Eindhoven University of Technology, Eindhoven, The Netherlands; Amsterdam Business School, Faculty of Economics and Business, University of Amsterdam, Amsterdam, The Netherlands.


    
====================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================





In this paper we analyze the distribution of the Age of Information (AoI) of a tagged data stream sharing a processor with a set of other data streams. We do so in the highly general setting in which the interarrival times pertaining to the tagged stream can have any distribution, and also the service times of both the tagged stream and the background stream are generally distributed. The packet arrival times of the background process are assumed to constitute a Poisson process, which is justified by the fact that it typically is a superposition of many relatively homogeneous streams. The first major contribution is that we derive an expression for the Laplace-Stieltjes transform of the AoI in the resulting GI+M/GI+GI/1 model. Second, we use stochastic ordering techniques to identify tight stochastic bounds on the  AoI. In addition, when approximating the tagged stream's inter-generation times through a phase-type distribution (which can be done at any precision), we present a computational algorithm for the mean AoI. 
As illustrated through a sequence of numerical experiments, the analysis enables us to assess the impact of background traffic on the AoI of the tagged stream. 




§ INTRODUCTION

In this paper we consider the situation of an information source that
is equipped with a sensor, feeding into a processor (often referred to
as a server), and a monitor (cf. Figure <ref>). Due to the fact that the (typically packetized) information has to be processed by the server, it arrives with some delay at the monitor. 
The Age of Information (AoI) <cit.> is a performance measure that quantifies the timeliness of the monitor’s knowledge of the information produced by the source. To conveniently analyze the AoI, queueing theory has proven to be particularly useful <cit.>.

In the theoretical analysis of the AoI, the focus has been primarily on a single data stream in isolation, in that the stream under consideration is the sole user of the processor <cit.>. 
A few exceptions are studies on the AoI in multi-source queueing systems with Poisson input, where inter-generation times of information packets are assumed to be exponentially distributed <cit.>.
In a more typical scenario, however, a tagged data stream corresponds to packet-generation intervals with relatively low variability, sharing the resource with a massive number of competing streams. This may entail that the characteristics of these background streams may significantly impact the AoI of the tagged stream. The main objective of this paper is to quantify this effect. In our analysis we rely on the well-known result that, with the number of background streams being typically large and relatively homogeneous, the corresponding aggregate packet arrival process can be accurately approximated by a Poisson process <cit.>. 





We proceed by discussing the three main contributions of this work.  

    
  ∘ The first part of our analysis considers the AoI of the tagged stream in the context of our general model, which we refer to as a queueing system of the GI+M/GI+GI/1 type. We succeed in evaluating the distribution of the AoI in this multi-source context through its Laplace-Stieltjes transform (LST), with the resulting expression being in terms of the (known) transform of the stationary system delay. 
    
  ∘ While techniques are available to accurately invert LSTs to obtain the underlying random variable's cumulative distribution function, one would like obtain more explicit insight into the impact of the background data streams on the AoI of the tagged stream. Relying on stochastic ordering results, we manage to derive closed-form lower and upper bounds on this AoI. 
    
  ∘ Informally, the class of phase-type distributions contains convolutions and mixture of exponential distributions. Owing to the fact that any distribution on the positive half-line can be approximated arbitrarily closely by a phase-type distribution
    <cit.>, we focus on the model in which the inter-sampling times of the tagged stream are of phase-type. We show how to stably compute the mean AoI, in the resulting multi-source model of the type PH+M/GI+GI/1. 


The model is more realistic than existing ones due to the fact that it allows renewal arrivals (of the tagged stream, that is), rather than just Poisson arrivals. In particular, it thus covers the practically relevant case when the tagged stream has constant interarrival times, and in addition the case that the interarrival times are `almost constant' (i.e., constant plus some perturbation). 

In previous work, all papers on multi-source models, involving a tagged stream sharing the processor with background streams, assume Poisson arrivals. The consequence is that one can deal with the heterogeneity of the service times in a straightforward manner: at any arrival one can probabilistically determine whether the arrived packet corresponds to the tagged stream or to the background stream. Put differently, for any multi-source M+M/GI+GI/1 model one can work with an equivalent single source M/GI/1 model. It is remarked that in <cit.> the AoI under renewal arrivals (i.e., the interarrival times stemming from an arbitray distribution) is extensively studied, but in the context of a single-source model.

We proceed by an (non-exhaustive) account of papers on AoI in the multi-source context, primarily focusing on papers in which the underlying queueing discipline is first-come first-served (FCFS). 
The pioneering paper <cit.> quantifies the AoI of a tagged stream interfering with a background stream, but in a setting in which interarrival times as well as service times are exponentially distributed. For the resulting M+M/M+M/1 model, and actually also for the analogous model in which more than two streams share the queueing resource, the authors succeed in providing a closed-form for the mean of the AoI of each of the streams. This analysis is extended in <cit.> to general service times (but still Poisson arrivals). The fact that it is reasonable to drop older packets when new packets arrive at the queue, justifies the interest in models with finite buffer capacity. In the single-stream case, it is natural to consider the case that the buffer can contain precisely one packet. The corresponding multi-source model (with a buffer size of one, that is) is dealt with in <cit.>, but is remarked that in this setting with multiple coexisting sources it would be more appropriate to work with a setup in which all interfering streams are given a dedicated buffer of size one.
In <cit.>, relying on Palm calculus techniques, an expression for the joint transform of the AoI s is obtained, which is valid under fairly general conditions. It should be noted, however, that its application to specific situations requires additional non-trivial analysis, and so far such analysis has succeeded only in special cases (namely, the input is Poisson and the service discipline corresponding to preemptive last-come first-served). We finally mention <cit.>, in which the uncertainty in the arrival and service processes is captured through the use of uncertainty sets. The analysis leads to bounds on the Peak-AoI (defined as the largest value of the AoI, before it drops due to a new update), accurate under both light and heavy loads.  This literature review confirms that, before our work, no results were available for the AoI distribution in systems of the practically relevant GI+M/GI+GI/1 type.


This paper has been organized as follows. Section <ref> introduces the model, convenient notation, and some preliminaries. Then Section <ref> provides an exact characterization of the AoI distribution through its LST. Explicit stochastic bounds on the AoI are presented in Section <ref>, including appealing lower and upper bounds on the mean AoI. In Section <ref> we specialize to the PH+M/GI+GI/1 systems, which, as argued above, is dense in the class of all GI+M/GI+GI/1 systems (in that any instance in the latter class can be approximated arbitrarily closely by an instance in the former class); we develop computational algorithms to obtain the mean AoI. In Section <ref> we systematically assess the performance of the PH+M/GI+GI/1 system, also quantifying the errors that we would face when we would have worked with the simpler M+M/GI+GI/1 system. Section <ref> gives concluding remarks, including directions for potential follow-up research. 



§ MODEL, NOTATION, AND PRELIMINARIES


Throughout this paper, we follow a convention that for a non-negative
random variable Y, its cumulative distribution function (CDF) and
probability density function (if exists) are denoted by F_Y(·)
and f_Y(·). We also define F_Y(·) as
its complementary CDF and f_Y^*(·) as its LST:


    F_Y(x) = 1 - F_Y(x),
      
    x ≥ 0,
        
    f_Y^*(s) 
    = 
    [e^-sY] = ∫_0^∞ e^-sx F_Y(x),
      
    s ≥ 0.



Suppose that a sensor observes (i.e., samples) the time-varying state of an
information source at a rate λ > 0.
At the sampling time instants, the sensor generates an information
packet containing the obtained data and transmits this to a remote monitor. 
Each information packet experiences some delay at an
intermediate communication channel or a processor which performs 
a computation to extract information from the raw data.
More specifically, let (α_n)_n=0,1,… denote the
sequence of generation times of information packets, where 
α_n-1≤α_n (n = 1,2,…).
We refer to the information packet generated at time t = α_n as
the nth packet. By definition, the inter-generation time G_n
between the (n-1)th and nth packets is given by


    G_n = α_n - α_n-1,
      
    n = 1,2,….


The nth packet (more precisely, the information it contains) is
received by the monitor at time β_n, i.e., the system delay
D_n experienced by the nth packet is given by


    D_n = β_n - α_n,
      
    n = 0,1,….


We hereafter assume that the time axis is taken so 
that


    α_0 > 0,


which is convenient in defining the AoI, as we will see later.





In this paper, we model the sequence of system delays
(D_n)_n = 0,1,… via the use of a single-server queue.
This queue is fed by a tagged input, corresponding to 
the sensor-monitor pair we focus on, as well as 
additional exogenous inputs, which represent the
background traffic load generated by sources other than the
tagged input (cf. Figure <ref>). More specifically, we
consider a queueing system in which the arrival process consists of 
two different input streams, in the sequel denoted
by (X(t))_t ≥ 0 and (X_(t))_t ≥ 0,
where X(t) and X_(t) denote the cumulative amount of work
brought into the system by the sensor and the background traffic, respectively.
The workload in system V(t) at time t is then given by


    V(t) = V(0) + X(t) + X_(t) - μ t + μ∫_u ∈ [0,t)𝕀V(u) = 0  u,
      
    t ≥ 0,


where μ> 0 denotes the service rate. 
See Figure <ref> for a graphical explanation. 
Assuming that the input processes (X(t))_t ≥ 0 and
(X_(t))_t ≥ 0 are càdlàg (i.e., right continuous with left limits),
so is the associated workload process (V(t))_t ≥ 0.



Letting H_n (n = 0,1,…) denote the service requirement of the
nth packet generated by the sensor, X(t) can be represented as


    X(t) = ∑_n=1^∞𝕀0 ≤α_n ≤ t H_n.


In this paper, we let the inter-generation times pertaining to the sensor be denoted by
(G_n)_n=1,2,…, and the corresponding service requirements by
(H_n)_n=0,1,…. We assume that these are two sequences of independent and identically distributed
(i.i.d.) random variables, with the respective CDFs F_G(·) and F_H(·), where its is in addition assumed that both sequences are independent of each other.
On the other hand, we assume that the background traffic
(X_(t))_t ≥ 0 is a compound Poisson process: it has
upward jumps following a Poisson process with rate
λ_, where the jump sizes are i.i.d. with CDF F_H_(·).
The background traffic load is then characterized as


    (X_(t) ≤ x) 
    =
    ∑_k=0^∞e^-λ_t(λ_ t)^k/k!·
    F_H_^⋆ k(x),
      
    x ≥ 0,


where F_H_^⋆ k(·) denotes the k-fold convolution of
F_H_(·), being defined recursively via


    F_H_^⋆ 1(x) = F_H_(x),
      
    F_H_^⋆ k(x) 
    = 
    ∫_0^x F_H_^⋆ k-1(x-y)  F_H_(y),
      
    k = 2,3,….


The processes (X(t))_t ≥ 0 and (X_(t))_t ≥ 0 are assumed independent.
The resulting queueing model could be denoted by GI+M/GI+GI/1, where the +
symbol stands for `superposition': a renewal input stream with generally distributed service times shares a single-server queueing resource with a Poisson input stream with generally distributed service times. As mentioned in the introduction, the Poissonian nature of the background arrivals is motivated by the fact that the number of background streams is typically large and fairly homogeneous. 

In terms of the workload process, the
system delay D_n of the nth packet is given by


    D_n = V(α_n-) + H_n/μ 
    = 
    V(α_n)/μ,


where the second equation follows from the right-continuity of
(V(t))_t ≥ 0. The AoI A(t) at time t is
defined as the elapsed time since the generation time of the last
packet received by the monitor (cf. Figure <ref>):


    A(t) = 
    {[                    t - A(0),                 0 ≤ t < β_0,; t - max{α_n; α_n + D_n ≤ t},                      t ≥β_0. ].


Recall that α_n + D_n = β_n represents the reception time
of the nth packet, so that t ≥β_0 in the second case of 
(<ref>) is required for the set {α_n; α_n + D_n ≤ t} to be
non-empty. Also, note that the condition (<ref>) ensures 
β_0 > 0. Let A_,n (n=1,2,…) denote the nth peak AoI, i.e.,
the value of the AoI just before the nth packet is received by the
monitor:


    A_,n 
    = 
    lim_t →β_n- A(t) 
    = 
    β_n - α_n-1   =
    β_n - α_n + α_n - α_n-1
       =
    D_n + G_n,


where the last equation follows from (<ref>) 
and (<ref>).

Let G, H, and H_ respectively denote generic random
variables following the CDFs F_G(·), F_H(·), and
F_H_(·). Also recall that λ := 1/[G] denotes the
sampling rate of the sensor. We define the traffic intensity of
respective streams as


    ρ := λ[H]/μ,
      ρ_ := λ_[H_]/μ.


In the rest of this paper, we assume that the stability condition


    ρ + ρ_ < 1,


is in place. Because the system is regenerative and the regeneration
time is almost surely (a.s.) finite under (<ref>), the time-averaged
distributions of (D_n)_n = 0,1,…,
(A_,n)_n=1,2,…, and (A(t))_t ≥ 0 agree with 
their respective stationary distributions F_D(·), F_A_(·),
and F_A(·), i.e., the following relations hold a.s.:


    lim_N →∞1/N∑_n=0^N-1𝕀D_n ≤ x   =
    F_D(x),
        
    x ≥ 0,
    
    lim_N →∞1/N∑_n=1^N𝕀A_,n≤ x   =
    F_A_(x), 
        
    x ≥ 0,
    
    lim_T →∞1/T∫_0^T 𝕀A(t) ≤ x t
       =
    F_A(x),
        
    x ≥ 0.


We thus define D, A_, and A as 
generic random variables corresponding to the stationary versions of D_n, A_,n, and A(t),
respectively.



§ EXACT FORMULA FOR THE AOI DISTRIBUTION

In this section we derive an exact expression for the LST of the
stationary AoI A, assuming that the tagged stream, corresponding to
the sensor-monitor pair, has a general inter-sampling time distribution. 

We utilize the following characterization of the stationary
distribution of the AoI <cit.>:


    F_A(x)    = 
    1/[G]∫_0^x 
    {F_A_(y) - F_D(y)} y,
      
    x ≥ 0,
    
    
    f_A^*(s)    = f_D^*(s) - f_A_^*(s)/s[G],
      
    s > 0.


Note here that (<ref>) implies [A_] - [D] =
[G], with which we can verify that F_A(·) and 
f_A^*(·) given by (<ref>) and (<ref>)
satisfy lim_x →∞ F_A(x) = 1 and lim_s ↓ 0 f_A^*(s)
= 1. By (<ref>) and (<ref>), the analysis
of the stationary AoI A is reduced to that of the stationary system
delay D and the stationary peak AoI A_. 

Let B_ denote the length of a busy period of
the M/G/1 queue only with the background traffic load; 
it is a standard result from queueing theory that its 
LST f_B_^*(s) satisfies the Kendall functional equation <cit.>


    f_B_^*(s) 
    = 
    f_H_^*( 
    s + λ_ - λ_f_B_^*(s)/μ).


It is easily seen that for each s ≥ 0, 
this equation and the condition |f_B_^*(s)| ≤ 1 
uniquely determine the value of f_B_^*(s); to this end, observe that 
f_B_^*(s) is the solution of the equation y=f^*(y) with a function

    f^*(y)=f_H_^*( 
    s + λ_ - λ_y/μ),


which is a convex function satisfying f^*(0)>0 and f^*(1)= f_H_^*(s/μ)<1.




The LST of the stationary peak AoI f_A_^*(·) is given in
terms of the LST of the stationary system delay f_D^*(·) by


    f_A_^*(s) 
    = 
    f_H^*(s/μ)
    ∫_y=0^∞
    e^-syℒ_ω^-1[
    1/ϕ(s) - ω(
    s/ψ(ω)·
    f_D^*(ψ(ω)) - f_D^*(s)
    )
    ] (y)
      F_G(y),


where ϕ(s) and ψ(ω) are defined as


    ϕ(s)    := s - λ_ + λ_f_H_^*(s/μ),
    
    ψ(ω)    := ω + λ_ - λ_f_B_^*(ω),


and ℒ_ω^-1[f(s,ω)] denotes the inverse
Laplace transform of f(s,ω) with respect to ω:


    f(s,ω)
    =
    ∫_0^∞ℒ^-1[f(s,ω)](t) e^-ω t t.




From (<ref>) and (<ref>), we have 


    A_,n   = V(α_n-) + H_n/μ + G_n
    = 
    V̂(α_n-)+ G_n
    + H_n/μ,


where V̂(t) denotes the virtual waiting time:


    V̂(t) = V(t)/μ.


Because (X(t))_t ≥ 0 has no jumps in t ∈ [α_n-1,
α_n), in this interval (<ref>) can be rewritten as


    V(t) = V(α_n-1) + X_(t) - μ t 
    + μ∫_u ∈ [α_n-1,t)𝕀V(u) = 0  u,
      
    t ∈ [α_n-1, α_n),


i.e., (V(t))_t ∈ [α_n-1, α_n) behaves as if it is the workload process of an ordinary M/GI/1
queue with the arrival rate λ_, the service requirement distribution
F_H_(·), and the service rate μ. 
Let (V_(t))_t ≥ 0 denote the workload process of such an
M/G/1 queue and let V̂_(t) = V_(t)/μ denote the
corresponding virtual waiting time. It is known that
the LST of V̂_(t) satisfies <cit.>


    ∫_t=0^∞[e^-sV̂_(t)]
    e^-ω t t
       =
    1/ϕ(s) - ω(
    s/ψ(ω)·
    f_V̂_(0)^*(ψ(ω))
    - f_V̂_(0)^*(s)
    ).


Therefore, we have from (<ref>),


    [e^-sA_,n]
       =
    f_H^*(s/μ)
    ·[e^-s(V̂(α_n-) + G_n)]
    
       =
    f_H^*(s/μ)
    ∫_x=0^∞[e^-s(V̂(α_n-) + G_n)|V̂(α_n-1) = x]
      F_V̂(α_n-1)(x)
    
       =
    f_H^*(s/μ)
    ∫_x=0^∞[e^-s(V̂_(G_n) + G_n)|V̂_(0) = x]
      F_D_n-1(x)
    
       =
    f_H^*(s/μ)
    ∫_x=0^∞ F_D_n-1(x)
    ∫_y=0^∞
    e^-sy [e^-sV̂_(y)|V̂_(0) = x]
      F_G(y)
    
       =
    f_H^*(s/μ)
    ∫_y=0^∞
    e^-sy{∫_x=0^∞[e^-sV̂_(y)|V̂_(0) = x]
      F_D_n-1(x)
    } F_G(y),


which together with (<ref>) implies (<ref>).





The LST f_A^*(·) of the stationary AoI is given in terms of the
LST f_D^*(·) of the stationary system delay by


    f_A^*(s) 
    = 
    f_H^*(s/μ)
    ∫_y=0^∞1 - e^-sy/[G]s·ℒ_ω^-1[
    1/ϕ(s) - ω(
    s/ψ(ω)·
    f_D^*(ψ(ω)) - f_D^*(s)
    )
    ] (y)
      F_G(y).




Similarly to the proof of Lemma <ref>, 
we observe from D_n = (V(α_n-) + H_n)/μ that its
distribution equals that of the sum of H/μ and
V̂_(G_n) conditioned on V̂_(0) = D_n-1.
We thus have for the stationary system
delay D (cf. (<ref>)),


    f_D^*(s) 
    = 
    f_H^*(s/μ)
    ∫_y=0^∞ℒ_ω^-1[
    1/ϕ(s) - ω(
    s/ψ(ω)·
    f_D^*(ψ(ω)) - f_D^*(s)
    )
    ] (y)
      F_G(y).


Therefore, we obtain (<ref>) by combining 
(<ref>), (<ref>), and
(<ref>).


The LST of the stationary AoI is thus given in terms of that of the
stationary system delay D. The characterization of D is known in
the literature <cit.>, which takes a form of an
intuitively appealing decomposition formula:



Consider an GI/GI/1 queue which has inter-arrival times following the CDF
F_G(·) and service times with LST 
f_H_⋆^*(s) := f_H^*( (s+λ_-λ_f_B_^*(s))/μ ) (s ≥ 0).
Let W_⋆ denote the stationary waiting time in this GI/GI/1 queue:


    f_W_⋆^*(s)
    =
    exp[
    -∑_k=1^∞1/k∫_0^∞
    (1-e^-sx)  F_U_k(x)
    ],
      
    s ≥ 0,


where U_k denotes the sum of k i.i.d. copies of the difference
H_⋆ - G of the inter-arrival and service times, i.e.,


    F_U_1(x) 
       = 
    ∫_0^∞(H_⋆≥ y-x)  F_G(y),
      
    x ∈ (-∞, ∞),
    
    
    F_U_k(x)
       =
    ∫_-∞^∞
    F_U_k-1(x-y)  F_U_1(y),
      
    x ∈ (-∞, ∞),
     
    k = 2,3,….


We then have for the original model with GI/GI and M/GI input
streams, the LST of the system delay D experienced by the GI/GI
stream is given by


    f_D^*(s) 
    = 
    (1-ρ_)s/s-λ_+λ_f_H_^*(s/μ)·
    f_W_⋆^*(s - λ_+λ_f_H_^*(s/μ))
    · f_H^*(s/μ).




The expression (<ref>) for the stationary waiting time
W_⋆ in the GI/GI/1 queue can be found in <cit.>.
Its mean [W_⋆] is also given explicitly by <cit.>:


    [W_⋆] = ∑_k=1^∞1/k∫_0^∞ x  
    F_U_k(x).




Lemma <ref> readily follows from
<cit.>. Note that although these papers show 
the decomposition formula for only the time-depedent
and stationary virtual waiting time V̂(t) and V, their
discussion clearly carries over to the virtual waiting time just before
arrivals of the GI/GI stream.


Theorem <ref> and Lemma <ref> completely
determine the AoI distribution. In particular, we can evaluate the
value of f_A^*(s) and its derivatives numerically, which enables us
to compute the moments and the CDF F_A(·) of the stationary 
AoI; we refer to e.g. <cit.> for Laplace inversion methods to numerically compute the corresponding CDF.
On the other hand, these formulas also enable us to gain further insights
into the way the model parameters affect the AoI performance. 
In the next section, we explore this direction in detail utilizing
stochastic ordering techniques.



§ STOCHASTIC ORDERINGS AND BOUNDS OF THE AOI


In this section, we derive stochastic comparison results and
closed-form bounds for the AoI based on the stochastic ordering
techniques. In Section <ref>, we briefly review the
definitions of several stochastic orders to be used. 
In Section <ref>, we prove several stochastic ordering
properties of the AoI, which leads to its distributional bounds. 
We then derive closed-form bounds on the mean AoI in Section <ref>.



 §.§ Stochastic Orders



Let us start by briefly reviewing the definitions of some stochastic orders (see
<cit.> for more detailed explanations).
Let Y_0 and Y_1 denote two non-negative random valiables.
Y_0 is said to be smaller than Y_1 in the usual stochastic
order (denoted by Y_0 ≤_ Y_1) when


    F_Y_0(x) ≤F_Y_1(x),
      ∀ x ≥ 0.


It is known that Y_0 ≤_ Y_1 if and only if
[h(Y_0)] ≤[h(Y_1)] holds for all non-decreasing functions
h(·) such that the expectations exist.

On the other hand, Y_0 is said to be smaller than Y_1 in the
increasing convex order (denoted by Y_0 ≤_ Y_1) when


    ∫_x^∞F_Y_0(y)  y 
    ≤∫_x^∞F_Y_1(y)  y,
      ∀ x ≥ 0.


As the name suggests, Y_0 ≤_ Y_1 if and only if 
[h(Y_0)] ≤[h(Y_1)] holds for all non-decreasing convex
functions h(·) such that the expectations exist.
Similarly, Y_0 is said to be smaller than Y_1 in the
increasing concave order (denoted by Y_0 ≤_ Y_1) when


    ∫_0^x F_Y_0(y)  y 
    ≤∫_0^x F_Y_1(y)  y,
      ∀ x ≥ 0,


which is equivalent to that [h(Y_0)] ≤[h(Y_1)] holds for all
non-decreasing concave functions h(·) such that the expectations exist.

Furthermore, Y_0 is said to be smaller than Y_1 in the
convex order (denoted by Y_0 ≤_ Y_1) when 


    Y_0 ≤_ Y_1    [Y_0] = [Y_1].


The convex ordering Y_0 ≤_ Y_1 holds if and only if [h(Y_0)] ≤[h(Y_1)]
for all convex functions h(·) such that the expectations exist.
Because [Y_i] = ∫_0^∞F_Y_i(y)  y, it is
readily seen that Y_0 ≤_ Y_1 is also equivalent to
Y_0 ≥_ Y_1 and [Y_0] = [Y_1]; note that the
inequalities have opposite directions here.

While the usual stochastic order ≤_ compares the magnitude of
random variables, the convex order ≤_ compares random
variables in terms of the variability given the same mean.
The increasing convex and increasing concave orders compare both the
magnitude and variability. Specifically, Y_0 ≤_ Y_1 implies
that Y_0 is smaller and less variable than Y_1, while 
Y_0 ≤_ Y_1 implies that Y_0 is smaller but more variable
than Y_1.

Notice however that ≤_ is stronger than ≤_ in the sense that 
Y_0 ≤_ Y_1 ⇒ Y_0 ≤_ Y_1 (also we have
Y_0 ≤_ Y_1 ⇒ Y_0 ≤_ Y_1).
On the other hand, if Y_0 ≤_ Y_1 and Y_0 ≤_ Y_1,
then Y_0 =_ Y_1 <cit.>, where =_
denotes the equality in distribution. Therefore, the usual stochastic
order and the convex order compare random variables from essentially
different viewpoints.



 §.§ Distributional Bounds for the AoI



As mentioned in the proof of Lemma <ref>, 
the virtual waiting time process (V̂(t))_t ∈ [α_n-1,α_n) during
an inter-sampling time behaves in the same way as the ordinary M/GI/1
queue with the arrival rate λ_ and the service time distribution
F_Ĥ_(·), where Ĥ_ =_ H/μ.
Recall that the virtual waiting time and the system delay are related as 
D_n = V̂(α_n).  We then define (Z(t; D_n))_t ≥ 0 as the
virtual waiting time process of the M/GI/1 queue with only the
background traffic, whose initial value equals to the system
delay D_n:


    F_Z(t; D_n)(x)
    =
    ∫_0^∞(V̂_(t) > x |V̂_(0) = x_0)
      F_D_n(x_0),
      
    x ≥ 0.


To make the notation simpler, we define


    F_Z(t, x; D_n)
    =
    F_Z(t; D_n)(x), 
      
    x ≥ 0.


Let W_n denote the waiting time of the nth packet 
(see (<ref>) and (<ref>)): 


    W_n = V(α_n-)/μ = V̂(α_n-).


By definition, we have 


    F_W_n(x) 
       = 
    ∫_0^∞F_Z(t, x; D_n-1)  F_G(t),
      
    n = 1,2,….


Note that (<ref>) along with 


    D_n = W_n + H_n/μ,


defines a recursive relation for D_n (n=0,1,…).

Imposing the following assumption on the initial virtual waiting time
V̂(0) largely simplifies the derivation of stochastic
comparison results. It is by no means restrictive, as such an assumption on V̂(0) does not affect any of the
results on stationary distributions.



The initial virtual waiting time V̂(0) is distributed according
to the stationary virtual waiting time distribution in the M/GI/1
queue with only the background traffic existing, i.e.,


    [e^-sV̂(0)] = (1-ρ_)s/s-λ_+λ_f_H_^*(s/μ),
      
    s > 0.




Assumption <ref> represents the situation that the system
had been running only with the background traffic until the first
generation time α_0 of the packet by the sensor.


Under this assumption, the following monotonicity property is known in the literature <cit.>.




If Assumption <ref> is satisfied, then Z(t; D_n) is
stochastically decreasing in t for each n, i.e.,


    F_Z(t_1,x; D_n) ≥F_Z(t_2,x; D_n),
      
    0 ≤ t_1 ≤ t_2,
      
    x ≥ 0,
      
    n = 1,2,….




Informally, the virtual waiting time process V̂(t) during
the time-interval t ∈ [α_n-1, α_n)
starts from the system delay D_n-1 = V̂(α_n-1) and it tends
toward the stationary distribution of the M/GI/1 queue.
Lemma <ref> shows that this transition is in fact
decreasing in the usual stochastic order.

Recall that the exact peak AoI distribution is characterized by the
rather complicated expression in (<ref>).
The following lemma shows that ignoring the dependence
between D_n and G_n in (<ref>) yields a simple upper
bound of the peak AoI in the convex order.




Under Assumption <ref>, the nth peak AoI A_,n
is bounded above in the convex
order as


    A_,n≤_ D_n + G,


where G denotes a generic random variable for inter-sampling times
and it is independent of the system delay D_n.




We consider the joint distribution of W_n and the inter-sampling
time G_n between (n-1)st and nth packets in the following form: 


    (W_n > x, G_n ≤ y)
        
       =
    ∫_x_0=0^∞(V̂_(G_n) > x, G_n ≤ y |V̂_(0) = x_0)
      F_D_n-1(x_0)
    
       =
    ∫_t=0^∞𝕀t ≤ y∫_x_0=0^∞(V̂_(t) > x |V̂_(0) = x_0)
      F_D_n-1(x_0)
      F_G(t)
    
       =
    ∫_t=0^∞𝕀t ≤ yF_Z(t,x; D_n-1)
      F_G(t),


where F_Z(t,x; D_n-1) is defined as in (<ref>) and (<ref>).
We then have


    (W_n > x, G_n ≤ y)
       =
    [F_Z(G,x; D_n-1)𝕀G ≤ y]
    
       ≥ [F_Z(G,x; D_n-1)][𝕀G ≤ y],


where the last inequality follows from the fact that 
for any random variable Y and non-increasing functions
g_1(·) and g_2(·), we have [g_1(Y)g_2(Y)] ≥[g_1(Y)][g_2(Y)] (see <cit.>).

Because (<ref>) implies 
(W_n > x) = [F_Z(x,G; D_n-1)],
we further obtain


    (W_n > x, G_n ≤ y)
    ≥(W_n > x)(G_n ≤ y),
      
    x ≥ 0,
      
    y ≥ 0.


Therefore, W_n and G_n are 
negatively quadrant dependent <cit.>.
It then follows from <cit.>
that for any non-decreasing functions
h_1(·) and h_2(·),


    [h_1(W_n), h_2(G_n)] ≤ 0.


We then have from <cit.>,


    W_n + G_n ≤_ W_n + G,


where G on the right-hand side is taken to be independent of
W_n. Because the convex order is preserved 
under the addition of an independent random variable <cit.>,
we obtain from (<ref>), (<ref>), and 
(<ref>), 


    A_,n+1   = W_n + G_n + H_n/μ
       ≤_
    W_n + G + H_n/μ
    =
    D_n + G.




As an application of Lemma <ref>, a similar bound for the
stationary version follows immediately, as stated in the next corollary. 




The stationary peak AoI A_ is bounded above in the convex
order as


    A_≤_ D + G,


where the stationary system delay D and the inter-sampling time
G are taken to be independent.


Since Corollary <ref> refers only to the stationary
distribution, Assumption <ref> is irrelevant here.


With Corollary <ref> and an additional analysis for
a lower bound, we obtain surprisingly simple stochastic bounds for the
stationary AoI:



The stationary AoI A is bounded by the sum of independent random variables as


    D^- + G≤_ A ≤_ D + G,


where D denotes the stationary system delay as before, D^- denotes
the system delay in the ordinary M/GI/1 queue with
only the background traffic existing, and G denotes the 
generic random variable for residual inter-sampling times:


    f_D^-^*(s)    = (1-ρ_)s/s-λ_+λ_f_H_^*(s/μ),
    
    
    f_G(x)    = F_G(x)/[G],
      
    f_G^*(s) = 1 - f_G^*(s)/s[G].








We first derive the upper bound.
From (<ref>), (<ref>), and Corollary <ref>, we have


    ∫_0^x F_A_(y)  y 
       =
    [A_] - ∫_x^∞F_A_(y)  y 
    
       =
    [D]+[G] - ∫_x^∞F_A_(y)  y
    
       ≥[D]+[G] - ∫_x^∞F_D+G(y)  y
    =
    ∫_0^x F_D+G(y)  y,
      ∀ x ≥ 0.


It then follows from (<ref>) that


    F_A(x) 
       ≥1/[G]{∫_0^x F_D+G(y)  y
    -
    ∫_0^x F_D(y)  y
    }
       =
    1/[G]{∫_y=0^x 
    (
    F_D(y)
    +
    ∫_u=0^y F_G(y-u)   F_D(u)
    )  y
    -
    ∫_0^x F_D(y)  y
    }
       =
    ∫_y=0^x 
    (
    ∫_u=0^y 
    F_G(y-u)/[G]  F_D(u)
    )
     y
    =
    F_D+G(x),
      ∀ x ≥ 0,


which proves the upper bound A ≤_ D+G.

We then consider the lower bound. Consider a modified workload process
(V^-(t))_t ≥ 0
where the workload brought by the sensor after time 0 is replaced by zero (cf. (<ref>)):


    V^-(t) = V(0) + X_^-(t) - μ t + μ∫_u ∈ [0,t)𝕀V^-(u) = 0 u,
      
    t ≥ 0.


It is obvious that we have V^-(t) ≤ V(t) (for t ≥ 0).
We define D_n^- (for n = 0,1,…) by


    D_n^- = V^-(α_n)+H_n/μ,


which, together with (<ref>), implies D_n^- ≤ D_n for each sample
path. Let us then define A^-(t) as


    A^-(t) := t - max{α_n; α_n + D_n^- ≤ t},


which clearly satisfies


    A^-(t) ≤ A(t),


for each sample path (cf. (<ref>)).
Furthermore, we define the corresponding peak AoI as (cf.(<ref>))


    A_,n^- 
    := 
    lim_t →α_n+D_n^- -A^-(t)
    =
    D_n^- + G_n.



Now suppose that the initial workload V(0) is distributed according
to the stationary distribution of (V^-(t))_t ≥ 0.
Under this condition, D_n^- defined in (<ref>) has the same
distribution as D^- and it is independent of G_n.
Therefore, we have from (<ref>) that


    F_A_,n^-(y)
    = 
    F_D^-(y)
    + ∫_0^y F_G(y-u)   F_D^-(u).


Applying (<ref>) to the modified AoI process
(A^-(t))_t ≥ 0, we obtain its stationary random variable
A^- as


    F_A^-(x) 
    =
    ∫_y=0^x
    (
    ∫_u=0^y F_G(y-u)/[G] F_D^-(u)
    ) y,


i.e., we have A^- =_ D^- + G.
Because (<ref>) still holds under this initial
condition, their time-averaged distributions satisfy 


    1/T∫_0^T 𝕀A^-(t) > x t 
    ≤1/T∫_0^T 𝕀A(t) > x t,
      
    T > 0, x ≥ 0.


Letting T →∞ in this equation (cf. (<ref>)),
we thus obtain the lower bound in (<ref>).




 §.§ Closed-Form Bounds for the Mean AoI



Theorem <ref> enables us to obtain closed-form bounds
for the mean AoI.



The mean AoI [A] is bounded by


    [A] 
       ≥λ_[H_^2]/2(1-ρ_)μ^2 
    + [H]/μ
    +
    [G^2]/2[G],
    
    [A] 
       ≤λ_[H_^2]/2(1-ρ_)μ^2 
    + [H]/μ
    + [G^2]/2[G]
    + (1-ρ_)[W_⋆],


where [W_⋆] is given by (<ref>).


From Lemma <ref> and Theorem <ref>, we can readily verify 
(<ref>) and (<ref>).



Recall that [W_⋆] denotes the mean waiting time in the GI/G/1 queue, where the LST of
the service time H_⋆ is given by f_H_⋆^*(s) =
f_H^*((s+λ_-λ_f_B_^*(s))/μ).
We can then replace [W_⋆] in the upper bound (<ref>)
by an expression only with the first two moments of the model parameters
<cit.>:


    [W_⋆] 
    ≤[H_⋆] +
    ρ_⋆(2-ρ_⋆)[G]/2(1-ρ_⋆)[G],


where


    ρ_⋆ = ρ/1-ρ_,
      [H_⋆] = [H]/μ(1-ρ_),
      [H_⋆] 
    = 
    [H]/μ^2 (1-ρ_)^2
    +
    ρ_[H_^2]/μ^2(1-ρ_)^3.





We can further simplify the upper bound by focusing
on the case where inter-sampling times have relatively small
variability. Note that choosing less variable inter-sampling times is reasonable
because it leads to the reduction in the AoI as known in the literature <cit.>.
Specifically, we derive a simpler bound on the mean AoI assuming
that the inter-sampling time distribution is new better than used in
expectation (NBUE), i.e.,


    [G-y  |  G > y] ≤[G],
      ∀ y ≥ 0.


It is also known that G is NBUE if and only if <cit.>


    G ≤_Exp([G]),


where Exp([G]) denotes an exponentially distributed
random variable with the same mean as G.
The class of NBUE distributions includes gamma distributions with
shape parameter not smaller than 1 and the deterministic distribution,
which are typical choices when modeling inter-sampling times in
monitoring systems.



If the inter-sampling time distribution F_G(·) is NBUE, 
we have


    λ_[H_^2]/2(1-ρ_)μ^2
    + 
    [H]/μ
    +
    [G^2]/2[G]≤[A] ≤λ[H]+λ_[H_^2]/2(1-ρ-ρ_)μ^2
    + 
    [H]/μ
    + 
    [G^2]/2[G].





Owing to Corollary <ref>, it is sufficient to show
that


    λ_[H_^2]/2(1-ρ_)μ^2
    +(1-ρ_)[W_⋆] 
    ≤λ[H^2]+λ_[H_^2]/2(1-ρ-ρ_)μ^2.


Recall that W_⋆ denotes the stationary waiting time in the GI/GI/1 queue as in Lemma <ref>, which has the inter-arrival time distribution F_G(·) and the modified service time distribution F_H_⋆(·) with LST f_H_⋆^*(s) =
f_H^*((s+λ_-λ_f_B_^*(s))/μ). From (<ref>), [W_⋆] is given by


    [W_⋆] = ∑_k=1^∞1/k·[max(0,U_k)],


where U_k is defined as the sum of k i.i.d. copies of the difference
H_⋆ - G of the inter-arrival and service times (see
(<ref>) and (<ref>)).
Because max(0, x) is a convex function of x, (<ref>) implies that if G is NBUE, then 
[W_⋆] is bounded above by the mean stationary waiting time in the M/G/1 queue with the same arrival rate 1/[G] and the same service time distribution, i.e.,


    [W_⋆] 
    ≤λ[H_⋆^2]/2(1-λ[H_⋆])
    =
    1/1-ρ_{λ[H^2] + λ_[H_^2]/2(1-ρ-ρ_)μ^2
    -
    λ_[H_^2]/2(1-ρ_)μ^2},


where the last equality can be verified with a straightforward calculation using


    λ[H_⋆] 
    = ρ/1-ρ_,
      λ[H_⋆^2]
    =
    ρλ_[H_^2]/(1-ρ_)^3μ^2
    +λ[H^2]/(1-ρ_)^2μ^2.


We then obtain (<ref>) from (<ref>).




§ PHASE-TYPE INTER-SAMPLING TIME DISTRIBUTION


In this section, we develop a numerical algorithm to compute the exact mean AoI, which is particularly useful in evaluating the tightness of the closed-form bound derived in the previous section.
To that end, we consider a special case where the inter-sampling times
(G_n)_n=1,2,… follow an M-state (M=1,2,…) phase-type
distribution <cit.> with representation (γ, Γ):


    F_G(x)    = 1 - γexp[Γ x] e,
      
    x ≥ 0,
    
    
    f_G(x)    = γexp[Γ x] (-Γ)e,
      
    x ≥ 0,


where e denotes an M× 1 vector with all elements equal
to one. The mean inter-sampling time is then given by


    [G] = γ (-Γ)^-1e.


The phase-type distribution represents the absorption time of a Markov
chain with the transition rate matrix Γ started from an
initial state distributed according to the probability vector
γ. By definition, we have 
γ≥ 0 and γe = 1. Also, Γ
has negative diagonal elements and non-negative off-diagonal elements.
Because the Markov chain governed by Γ is absorbing,
Γ is a defective transition rate matrix, i.e., its row
sums are non-positive (i.e., Γe≤ 0) with at least
one of them being strictly negative. 

It is noted that one can approximate any distribution on the positive
half-line arbitrarily closely by a distribution from the phase-type
class described above <cit.>. This legitimates considering the queueing model PH+M/GI+GI/1 (in self-evident notation), rather than the formally more general model GI+M/GI+GI/1. 

We impose the common assumption that 
the transition rate matrix Γ +
(-Γe)γ is irreducible, where 
(-Γe)γ corresponds to the event of
absorption and the transition to an initial state for the next inter-sampling time.
It is then ensured that  Γ + (-Γe)γ 
has a unique stationary probability vector π, which  solves


    π{Γ + (-Γe)γ}
    = 0,
      πe = 1.


By multiplying (Γ^-1)e from the right on both sides
of this equation and rearranging terms, we can verify that 


    π(-Γe)
    = 
    1/γ (-Γ)^-1e
    =
    λ,


where the last equality follows from
(<ref>) and λ = 1/[G].

In this model, the total input process 
(X_(t))_t ≥ 0 := (X(t) + X_(t))_t ≥ 0
follows a compound Markovian arrival process (MAP).
To be more specific, let C and D(x) (for x ≥ 0) denote
M × M matrices defined as 


    C = Γ - λ_I,
      D(x) 
    = 
    (-Γe)γ F_H(x) 
    + λ_ F_H_(x)I,


where I denotes an M × M unit matrix.
Let S(t) ∈{1,2,…,M} denote the phase of the system at time
t, i.e., the state of the phase-type distribution of the inter-sampling time.
We define x_(t, x) as a 1 × M vector whose
jth element j=1,2,…,M is given by


    [x_(t)]_j = (X(t) ≤ x, S(t) = j).


The vector (X_(t), S(t))_t ≥ 0 forms a bivariate Markov process and its 
transition law is completely determined by C and D(·): as t ↓ 0+,


    x_(t+ t, x)
    = 
    x_(t, x)  (I+C t)
    +
    ∫_0^x x_(t, x-y)  D(y)   t
    +o( t).


The first and second term on the right-hand side of this
equation correspond respectively the state-transitions without and
with an upward jump in X_(t).


The following results are known for this model <cit.>.
Recall that V(t) denotes the workload in the system, and that
V̂(t)/μ denotes the virtual waiting time. 
Let V̂ and S denote generic random variables following the
(joint) stationary distribution of (V̂(t), S(t)). 
Let v̂(x) (for x ≥ 0) denote a 1 × M vector whose
jth (j = 1,2,…,M) element represents the stationary joint
probability that the virtual waiting time is not greater than x and the phase 
equals j:


    [v̂(x)]_j = (V̂≤ x, S=j).


We define D^*(·) and v̂^*(·) as the LSTs of 
D(·) and v̂(·):


    D^*(s)    = ∫_0^∞ e^-sxD(x),
      
    s ≥ 0,
    
    v̂^*(s)    = ∫_0^∞ e^-sxv̂(x),
      
    s ≥ 0.


We further define and D̂(·) 
and D̂^*(·) as


    D̂(x)    = D(μ x),
      
    x ≥ 0,
    
    D̂^*(s) 
       = 
    ∫_0^∞ e^-sxD̂(x) = D^*(s/μ),
      
    s ≥ 0,


which correspond to D(·) and D^*(·) with
H and H_ replaced with H/μ and H_/μ.




Let (Q^(n))_n=0,1,… denote a sequence of M × M
matrices defined as


    Q^(0) = C,
      Q^(n) 
    = 
    C + ∫_0^∞D̂(x) exp[Q^(n-1)x],
      
    n = 1,2,….


(Q^(n))_n=0,1,… is an elementwise non-decreasing sequence
and its limit Q := lim_n →∞Q^(n) represents a
proper transition rate matrix of an irreducible continuous-time Markov
chain with M states. Let κ denote the stationary
probability vector associated with Q:


    κQ = 0,
      κe = 1.


The vector LST v̂^*(s) of the stationary virtual waiting
time is then given by


    v̂^*(s) =
    (1-ρ-ρ_)κ[I - R^*(s)]^-1,
      
    s ≥ 0,


where R^*(s) is defined as


    R^*(s)
    =
    ∫_x=0^∞e^-sx x ∫_y=x^∞D̂(y) exp[Q(y-x)],


and for s ≥ 0, eigenvalues of R^*(s) have absolute values strictly less
than one, i.e., 
    [I-R^*(s)]^-1 = ∑_m=0^∞
    (R^*(s))^m < ∞
 holds.



For this model, we also have a more straightforward generalized version of the
Pollaczek–Khinchine formula, as follows <cit.>:


    v̂^*(s)[sI+C+D̂^*(s)] 
    =
    (1-ρ-ρ_)s κ.


This equation is readily deduced from (<ref>) and
(<ref>), by noting that (<ref>) implies


    [I-R^*(s)](sI+Q) =
    sI+C+D̂^*(s).


However, the advantage of the expression (<ref>) over
(<ref>) is that [I-R^*(s)]^-1 exists for
all s ≥ 0 as stated in Lemma <ref>.
This is not the case for sI+C+D̂^*(s) because 
(<ref>) implies that 
(sI+C+D̂^*(s)) = 0 when s equals an
eigenvalue of -Q, which all have non-negative real parts owing
to the Perron–Frobenius theorem.


For reasons of numerical stability, the matrix exponentials exp[Q^(n)x]
and exp[Qx] in (<ref>) and (<ref>)

should be computed using the uniformization
technique <cit.>:


    exp[Q^(n) x]    = 
    ∑_k=0^∞(Q^(n) x)^k/k!
    =
    ∑_k=0^∞e^-ζ x (ζ x)^k/k!· 
    (I+ζ^-1Q^(n))^k,
    
    exp[Qx] 
       =
    ∑_k=0^∞e^-ζ x (ζ x)^k/k!· 
    (I+ζ^-1Q)^k,


where ζ denotes the maximum absolute value of the diagonal elements of C:


    ζ = max_i |[C]_i,i| = max_i |[Γ]_i,i|+λ_.


Note here that I+ζ^-1Q^(n) represents a substochastic matrix,
i.e., it has non-negative elements with row sums not exceeding one.
Therefore, we rewrite the recursion (<ref>) as


    Q^(n)   =
    C
    +
    ∑_k=0^∞∫_0^∞e^-ζ x (ζ x)^k/k!D̂(x) 
    · 
    (I+ζ^-1Q^(n-1))^k
    
       =
    C +
    ∑_k=0^∞D̂^(k)·
    (I+ζ^-1Q^(n-1))^k,
      
    n = 1,2,…,


where


    D̂^(m)   := 
    ∫_0^∞e^-ζ x(ζ x)^m /m!D̂(x) 
    
       =
    (-Γe)γ∫_0^∞e^-ζ x/μ(ζ x/μ)^m /m! F_H(x)
    +
    Iλ_∫_0^∞e^-ζ x/μ(ζ x/μ)^m /m! F_H_(x)
    
       =
    (-Γe)γĥ^(m)
    + Iλ_ĥ_^(m),

and

    ĥ^(m) 
    := 
    ∫_0^∞e^-ζ x/μ(ζ x/μ)^m /m! F_H(x),
      ĥ_^(m) 
    := 
    ∫_0^∞e^-ζ x/μ(ζ x/μ)^m /m! F_H_(x).


Here it is observed that the second equality follows from the definitions
(<ref>) and (<ref>) of D(·) and
D̂(·). It is in addition noted that ĥ^(m) (and ĥ_^(m), respectively) can be interpreted as the probability
that a Poisson process of rate ζ generates 
m arrivals in a
random interval of length H/μ (H_/μ, respectively). 
As a consequence, the
computation of ĥ^(m) and ĥ_^(m) is usually straightforward given specific distributions of the random variables H and H_ bg.

Recall that the LST of the stationary AoI is given by
(<ref>), where the LST f_D^*(s) of the stationary system delay D
is obtained from the conditional PASTA property <cit.>: 


    f_D^*(s) 
    = 
    v̂^*(s)(-Γe)f_H^*(s/μ)/λ.


Therefore, it follows from (<ref>) that 
the LST f_A^*(s) of the stationary AoI distribution is
given in terms of the probability density function of the inter-sampling time G by


    f_A^*(s) 
    = 
    f_H^*(s/μ)
    ∫_y=0^∞1-e^-sy/s[G]· p(s,y)  f_G(y)  y,


where p(s,y) := ℒ_ω^-1[ p^*(s,ω)
] (y) denotes the inverse Laplace transform of 


    p^*(s, ω)
    :=
    1/ϕ(s) - ω(
    s/ψ(ω)·
    f_D^*(ψ(ω)) - f_D^*(s)
    ),


i.e., 


    p^*(s,ω) = ∫_0^∞ p(s,t)  e^-ω t  t.


Note here that p^*(s, t) represents the LST of the transient virtual waiting time
of the M/GI/1 queue with only the background traffic, whose initial
value follows the same distribution as the stationary system delay D
(cf. (<ref>)):


    p(s,t)
    =
    [e^-sV̂_(t)|V̂_(0) ].


The LST of the AoI f_A^*(s) is thus determined completely by 
(<ref>), (<ref>), 
(<ref>), (<ref>), (<ref>), and (<ref>).
Similarly, we have the following expression for the mean AoI [A]:


    [A]    = lim_s → 0+ (-1)· f_A^*(s)/ s
       =
    [H]/μ
    +
    1/[G]∫_y=0^∞( y/2 + q(y))y
    · f_G(y)  y,


where q(y) (for y ≥ 0) is defined via


    q(y) := (-1)·lim_s → 0+∂ p(s,y)/∂ s.


From (<ref>), we can verify that q(y) represents the
mean transient virtual waiting time (at time y, that is) of the M/GI/1 queue with only the
background traffic, whose initial value equals to the stationary system delay D:


    q(y) = [V̂_(y) |V̂_(0) ].


Because it follows from (<ref>) that


    ∫_0^∞ q(y) e^-ω y y
       =
    (-1)·∫_0^∞lim_s → 0+∂ p(s,y)/∂ s· e^-ω y y
    
       =
    (-1)·lim_s → 0+∂ p^*(s,ω)/∂ s
       =
    -1-ρ_/ω^2
    +
    1/ω(
    f_D^*(ψ(ω))/ψ(ω)
    +[D]
    ),


q(y) is given by


    q(y) = ℒ_ω^-1[
    -1-ρ_/ω^2
    +
    1/ω(
    f_D^*(ψ(ω))/ψ(ω)
    +[D]
    )
    ].



Therefore, evaluating the expression for [A], as given in (<ref>), requires Laplace inversion
as well as integration over y∈[0,∞).
It is known that the accuracy of Laplace-inversion methods is highly
dependent on the shape of the function under consideration, which
sometimes leads to inconsistent numerical results between different
inversion algorithms <cit.>.
This led us to develop an alternative approach that utilizes our
assumption of phase-type distributed inter-sampling times.

The key idea in this alternative approach lies in expanding f_G(·), as given by (<ref>), as an
Poisson-weighted average of non-negative numbers (cf. the uniformization
technique <cit.>); the favorable numerical properties of the resulting computation scheme are due to the fact that the weights used all have the same sign. Recall that ζ was defined in
(<ref>). We define θ := max_i |[Γ]_i,i
| as the maximum absolute value of diagonal elements of Γ,
i.e., we have


    ζ = θ + λ_.


We then rewrite f_G(x) as


    f_G(x) 
       = 
    γ e^-θ xexp[(θI+Γ) x] (-Γ)e
       =
    γ e^-θ x∑_k=0^∞[(θI+Γ) x]^k/k!·
    (-Γ)e
       =
    ∑_k=0^∞e^-θ x (θ x)^k/k!·γ (I+θ^-1Γ)^k (-Γ)e = ∑_k=0^∞e^-θ x (θ x)^k/k!·
    c_k,


where c_k denotes the non-negative number given by


    c_k = 
    γ (I+θ^-1Γ)^k (-Γ) e,
      
    k = 0,1,….


With this expression, we rewrite (<ref>) as


    [A]
    =
    [H]/μ
    +
    ∑_k=0^∞c_k/θ[G]∫_y=0^∞(
    y/2
    +
    q(y) 
    )
    · y ·e^-θ y(θ y)^k/k!·θ  y,



As an immediate consequence of the definition of the gamma-integral,
we have


    ∫_y=0^∞ y^2 ·e^-θ y(θ y)^k/k!·θ  y
    =
    (k+1)(k+2)/θ^2


Also,


    ∫_0^∞
    q(y) 
    · y ·e^-θ y(θ y)^k/k!·θ  y
       =
    (k+1)∫_0^∞
    q(y) 
    ·e^-θ y(θ y)^k+1/(k+1)! y
    = 
    (k+1) q_k+1,


where q_k (for k=0,1,…) is defined as


    q_k := ∫_0^∞
    q(y) 
    ·e^-θ y(θ y)^k/k!  y.


We then obtain from (<ref>),


    [A]
    =
    [H]/μ
    +
    ∑_k=0^∞c_k/θ[G]{(k+1)(k+2)/2θ^2
    +
    (k+1) q_k+1}.


Therefore, what remains is to determine q_k (k=0,1,…). We
take a transform approach to that end: let q^*(z) (for |z| < 1)
denote the z-transform of q_k: 


    q^*(z) 
    := 
    ∑_k=0^∞ q_k z^k
    =
    ∫_0^∞
    q(y) e^-(θ-θ z)y y.


From (<ref>),


    q^*(z) 
       = 
    (-1)·lim_s → 0+∂/∂ s p^*(s, θ-θ z)
    
       =
    (-1)·lim_s → 0+∂/∂ s[
    1/ϕ(s) - (θ-θ z)(
    s/ψ(θ - θ z)·
    f_D^*(ψ(θ - θ z)) - f_D^*(s)
    )
    ]
    
       =
    -1-ρ_/(θ - θ z)^2
    +
    1/θ-θ z(
    f_D^*(ψ(θ - θ z))/ψ(θ - θ z)
    +[D]
    ).



To determine the coefficients q_k of q^*(z), we first consider
coefficients of ψ(θ - θ z) as a function of z.
By definition (<ref>), ψ(θ - θ z) is written as


    ψ(θ - θ z)    = 
    θ - θ z + λ_ - λ_f_B_^*(θ - θ z)
    
       =
    θ + λ_ 
    - (θ + λ_)
    ·{θ/θ + λ_· z 
    + λ_/θ + λ_· f_B_^*(θ - θ z)
    }
       =
    ζ - ζ b(z),


where we used (<ref>) and


    b(z) := θ/θ + λ_· z 
    + λ_/θ + λ_· f_B_^*(θ - θ z),


which is a proper probability generating function. 
We then define b^(i) (i = 0,1,…) as the ith coefficient of
b(z):


    b(z) = ∑_i=0^∞ b^(i) z^i.


The coefficients b^(i) are determined via Lemma <ref>; 
recall that ζ and θ are related as in (<ref>).
The proof of Lemma  <ref> is provided in
Appendix <ref>.




b^(i) (i = 0,1,…) is given by


    b^(i) 
    = 
    𝕀i=1θ/θ+λ_
    +
    λ_/θ+λ_∑_k=i^∞1/k-i+1ki(θ/θ + λ_)^i
    (λ_/θ + λ_)^k-i
    y_k-i+1^(k),


where y_k^(m) given recursively by


    y_1^(m)   = ĥ_^(m),
      
    m = 0,1,…,
    
    
    y_k^(m)   = ∑_i=0^m y_k-1^(i)ĥ_^(m-i),
      
    k = 1,2,…, 
    m = 0,1,…,


with ĥ_^(m) defined as (<ref>).





We then proceed to computing f_D^*(ψ(θ - θ z))/ψ(θ - θ z)
on the right-hand side of (<ref>).
We define b_ℓ^(k) (for ℓ = 0,1,…, 
k=0,1,…) as


    b_0^(k) = 𝕀k=0,
      
    b_ℓ^(k) = ∑_i=0^k b_ℓ-1^(k-i) b^(i),
      ℓ = 1,2,….


By definition (<ref>), it is obvious that


    {b(z)}^ℓ = ∑_m=0^∞ b_ℓ^(m) z^m.


Let (d^(m))_m = 0,1,… denote the coefficients of
f_D^*(ζ - ζ z) as a function of z:


    f_D^*(ζ - ζ z) = ∑_m=0^∞ d^(m) z^m.


We then have from (<ref>) and (<ref>),


    f_D^*(ψ(θ - θ z))/ψ(θ - θ z)
    =
    f_D^*(ζ - ζ b(z))/ζ - ζ b(z)   =
    1/ζ∑_k=0^∞{b(z)}^k
    ∑_m=0^∞ d^(m){b(z)}^m
    
       =
    1/ζ∑_ℓ=0^∞{b(z)}^ℓ∑_m=0^ℓd^(m)
       =
    ∑_i=0^∞
    z^i
    ·1/ζ∑_ℓ=0^∞ 
    b_ℓ^(i)∑_m=0^ℓd^(m).


where we put ℓ = m+k to obtain the second equality.

The sequence d^(m) (m=0,1,…) is determined by Lemma <ref> and Corollary <ref> below. The proof of Lemma <ref> is provided in
Appendix <ref>. 




Let v̂^(m) (for m=0,1,…) denote the mth
coefficient of v̂^*(ζ-ζ z) as a function of
z: 


    v̂^*(ζ-ζ z) = ∑_m=0^∞v̂^(m)
    z^m.


v̂^(m) (m=0,1,…) is determined by the following recursion:


    v̂^(0)   = (1-ρ-ρ_)κ
    (I-R^(0))^-1,
    
    v̂^(m)   = (∑_k=0^m-1v̂^(k)R^(m-k))
    (I-R^(0))^-1,
      
    m = 1,2,…,


where R^(m) (for m=0,1,…) is given by


    R^(m) 
    = 
    ζ^-1∑_k=0^∞D̂^(m+k+1)
    (I+ζ^-1Q)^k,


with Q given in Lemma <ref> and 
D̂^(m) (m=0,1,…) given by (<ref>).


The existence of (I-R^(0))^-1 is ensured by 
Lemma <ref> because R^(0) = R^*(ζ).







The coefficients d^(m) (m=0,1,…) are given by 


    d^(m) 
    = 
    1/λ∑_k=0^m v̂^(k)
    (-Γe)ĥ^(m-k),
      
    m =0,1,…,


where ĥ^(m) (m = 0,1,…) is given by
(<ref>).
Furthermore, the mean system delay [D] is given by


    [D] = 1/ζ∑_m=0^∞ md^(m).




(<ref>) follows immediately from 
(<ref>) and Lemma <ref>.
Also, (<ref>) is obvious because 


    [D] = (-1) ·lim_s → 0+/ s f_D^*(s)
    = 
    1/ζlim_z → 1/ z f_D^*(ζ - ζ
    z).




We are now in a position to present a (numerically stable) computation scheme for the mean
AoI [A].


The mean AoI is given by (<ref>) with q_k
(k = 0,1,…) determined recursively by 


    q_0 
       = 
    -1-ρ_/θ^2
    +
    [D]/θ + 1/θζ∑_ℓ=0^∞ b_ℓ^(0)∑_m=0^ℓd^(m),
    
    
    q_k
       =
    q_k-1
    -1-ρ_/θ^2
    + 
    1/θζ∑_ℓ=0^∞ 
    b_ℓ^(k)∑_m=0^ℓd^(m),
      
    k = 1,2,…,


where d^(m) (m = 0,1,…) and [D] are given in Corollary
<ref>. 



We have from (<ref>) and (<ref>),


    q^*(z) 
    =
    ∑_k=0^∞
    z^k
    (
    -(k+1)(1-ρ_)/θ^2
    +[D]/θ
    +1/θ∑_i=0^k 
    1/ζ∑_ℓ=0^∞ 
    b_ℓ^(i)∑_m=0^ℓd^(m)),


i.e.,


    q_k 
    =
    -(k+1)(1-ρ_)/θ^2
    +[D]/θ
    +1/θ∑_i=0^k 
    1/ζ∑_ℓ=0^∞ 
    b_ℓ^(i)∑_m=0^ℓd^(m),
      
    k = 0,1,….


It is then readily verified that (q_k)_k=0,1,… satisfies the
recursion (<ref>) and (<ref>).


In Figure <ref>, we summarize the resulting computational
algorithm for [A].
We note that in the computations it is necessary to truncate some infinite series; cf. (<ref>) and (<ref>).
To check the numerical precision of the procedure, the interpretation
(<ref>) of q(y) as the
mean transient workload in the M/GI/1 queue is useful.
More specifically, we can verify from (<ref>) that 


    lim_k →∞ q_k 
       = 
    1/θ·lim_k →∞∫_0^∞
    q(y) 
    ·e^-θ y(θ y)^k θ/k!  y
    
       =
    1/θ·lim_k →∞[
    q(E_k+1)
    ]
    
       = 
    1/θ·lim_y →∞ q(y),


where E_k denotes a random variable following a k-stage Erlang distribution
with rate θ.
Because q(y) represents the transient virtual waiting time in the
M/GI/1 queue with only the background stream (cf. (<ref>)),
it tends to the stationary mean virtual waiting time in the M/GI/1
queue as y →∞. We thus have


    lim_k →∞ q_k 
    =
    1/θ·λ_[H_^2]/2(1-ρ_)μ^2.


Therefore, for large values of k, the difference between the computed q_k and the right-hand side of (<ref>) serves as a
straightforward indication of the loss in numerical precision due to the
truncation.





§ NUMERICAL EXPERIMENTS


In this section we provide a sequence of numerical experiments that systematically assess the impact of the model parameters on the expectation of the AoI. 
We work with the setup of Section <ref> where the inter-sampling times
(G_n)_n=1,2,… follow a phase-type distribution. 
The idea is that we rely on the two-moments fit advocated in <cit.>, which identifies a convenient low-dimensional phase-type distribution with a given mean and variance.
With the coefficient of variation of the random variable G defined as

    s_G = √(𝕍 ar[G]/(𝔼[G])^2),

it distinguishes between the cases s_G≤ 1 and s_G>1. 
In the former case, we fit a mixture of two Erlang distributions: with probability p the fitted distribution is of Erlang(k,ν) type, and with probability 1-p of Erlang(k+1,ν) type. 
The integer-valued parameter k is chosen such that s_G lies in the interval between 1/(k+1) and 1/k, and then the parameters p and ν are chosen such that the mixed Erlang random variable has the given mean and variance. In the latter case, we fit a mixture of two Exponential distributions: with probability p the fitted distribution is exponential with parameter ν_1, and with probability 1-p exponential with parameter ν_2. 
Under the additional constraint of `balanced means', one can uniquely pick the parameters p, ν_1, and ν_2 so that the resulting random variable has the given mean and variance.

In a first experiment we study the effect of the `variability' of the inter-sampling times. 
Through this experiment we can quantify the error we would make by wrongly assuming that these inter-sampling times are exponentially distributed. 
We have picked the following parameters.
The sampling rate λ of the tagged stream, which is the reciprocal of the mean inter-sampling time 𝔼[G], was put equal to 0.05.
The (Poissonian) arrival rate λ_ bg corresponding to the background traffic equals 0.8, 0.85, and 0.9.
The service times of the tagged stream as well as the background stream have been chosen deterministic of size 1, and time is normalized such that the service rate μ equals 1.
Figure <ref> shows the mean AoI as a function of the coefficient of variation s_G pertaining to the inter-sampling times. 
As expected, 𝔼[A] is an increasing function of s_G: more variability typically leads to longer queues.
The more important conclusion, however, is that misspecifying our GI+M/GI+GI/1 model by its more elementary M+M/GI+GI/1 counterpart, that was covered by <cit.>, may lead to significant errors.
Indeed, for the parameters considered, the difference between 𝔼[A] in the M+M/GI+GI/1 case (i.e., the special case of the GI+M/GI+GI/1 system in which s_G=1) and the GI+M/GI+GI/1 can already be a factor 2. 

From Figure <ref>, we also observe that these three curves with
different ρ_ appear to be in parallel with each other, i.e.,
the variability s_G of the inter-sampling times has a common impact
on [A] regardless of the background traffic intensity. Notice that
the mean AoI is approximated by the sum of the
mean forward recurrence time [G^2]/(2[G]) = (1+s_G^2)[G]/2 of
inter-sampling times and the mean system delay; cf. <cit.>, and also
(<ref>)–(<ref>) in the present paper. 
Therefore, the effect of s_G on [A] is roughly evaluated as
(1+s_G^2)[G]/2; this is indeed verified in Figure
<ref>, where the `residual component' [A]- (1+s_G^2)[G]/2 as a
function of s_G is shown to be almost constant.





In the second experiment, we pick the parameters as before, but now
with s_G being fixed at 0.3 and varying λ. The result is
shown in Figure <ref>, where the explicit bounds derived in
Theorem <ref> are also plotted.
We observe a pattern that was found for more elementary AoI models: the mean AoI first decreases in λ, and later increases. 
The explanation for the phenomenon is that small values of λ mean that the monitor is provided with relatively little information, leading to untimely knowledge and hence a large AoI, whereas for large values of λ there will be more queueing delay and hence also a large AoI.
The shape of the graph interestingly means that there is an optimal value of λ, i.e., a value of λ that minimizes 𝔼[A].
In this experiment we also observe that the optimal λ decreases in λ_ bg.
To understand this, recall that the mean AoI is approximated by the
sum of the mean forward recurrence time of inter-sampling times and the mean system delay as mentioned above. Then observe that the former is a decreasing function of the sampling rate λ, whereas the latter is an increasing function of λ.
For smaller values of the rate λ_ bg, the system delay significantly increases from a smaller value of λ on, entailing that the optimal λ becomes smaller.

In Figure <ref>, we observe that the explicit upper bound well
resembles the exact curve of the mean AoI [A] as a function of the
sampling rate λ. The lower bound, on the other hand,
significantly underestimates
[A]; this is in particular the case when λ has a large impact on the system congestion, 
i.e., in a situation where ρ is close to the
`remaining capacity' 1-ρ_ of the system.
To obtain further insights related to this perspective, we set
λ=0.05 and scale up λ_ and μ simultaneously
while keeping ρ_ fixed, where the other parameters are the same as those underlying Figure <ref>. Figure <ref> plots the exact value of [A] as well as the corresponding bounds, as functions of the occupancy ρ_
/(ρ+ρ_) of the background traffic. 
From this figure, we observe that the explicit bounds in Theorem <ref> are tight when the background traffic is
dominant and the single sensor accounts for at most 1% of the traffic; 
in practice, modern broadband communication systems simultaneously serve vast
numbers of data flows, thus resulting in the occupancy
of a single sensor typically dropping far below 1%. In this sense, Theorem <ref> provides a simple yet highly accurate measure of
`information freshness' for monitoring systems with coexisting background
streams.








§ DISCUSSION AND CONCLUDING REMARKS

This paper has described a methodology to analyze the AoI in the context of the highly general GI+M/GI+GI/1 model. In a first main result, Theorem <ref>, we provide a closed-form expression for its Laplace-Stieltjes transform, enabling the numerical evaluation of the AoI moments. The second main result (Theorem <ref>) concerns stochastic lower and upper bounds on the AoI, leading to explicit, insightful lower and upper bounds on the mean AoI (Theorem <ref>). We then approximate the tagged stream's inter-generation times via a phase-type distribution, which can be done at any desired precision level. For the resulting PH+M/GI+GI/1 model we succeed, as a third main result, in developing a stable computational algorithm for the mean AoI; see Theorem <ref> for the resulting expression.

While with our GI+M/GI+GI/1 model we have reached a great level of generality, some further extensions can be thought of. In our setup we have, to model the background packet arrival process, relied on the fact that the superposition of many relatively homogeneous streams behaves essentially Poissonian. This means that we have to develop alternative techniques for settings in which the background stream does not correspond to a large traffic aggregate. A relevant specific case could be the one in which the background stream is modeled as the superposition of relative few streams with deterministic inter-arrival times, each of them with a random `phasing' (cf. the N*D/D/1 queue <cit.>). 

A second challenge could lie in the algorithmic computation of the AoI moments. Concrete objectives could concern optimizing the precision and run time of the algorithm for the mean AoI, as presented in Figure <ref>, and its extension to higher moments. 





§ PROOF OF LEMMA <REF>



Using the integral representation of the M/GI/1 busy period
length <cit.>, we have


    f_B_^*(θ - θ z)
       =
    ∫_0^∞ 
    e^-(θ-θ z)x∑_k=1^∞e^-λ_ x (λ_ x)^k-1/k! F_Y_k(x),


where Y_n denotes the sum of k i.i.d random variables following
the CDF F_H_/μ(·). Therefore, 


    f_B_^*(θ - θ z)
       =
    ∑_m=0^∞
    z^m
    ∑_k=1^∞∫_0^∞e^-θ x(θ x)^m/m!·e^-λ_ x (λ_ x)^k-1/k! F_Y_k(x)
    
       =
    ∑_m=0^∞
    z^m
    ∑_k=1^∞1/k·(m+k-1)!/m!(k-1)!·(θ/θ + λ_)^m
    (λ_/θ + λ_)^k-1
       ·∫_0^∞e^-(θ+λ_) x{(θ + λ_)
    x}^m+k-1/(m+k-1)! F_Y_k(x)
    
       =
    ∑_m=0^∞
    z^m
    ∑_k=1^∞1/km+k-1m(θ/θ + λ_)^m
    (λ_/θ + λ_)^k-1
    y_k^(m+k-1)
       =
    ∑_m=0^∞
    z^m
    ∑_k=m^∞1/k-m+1km(θ/θ + λ_)^m
    (λ_/θ + λ_)^k-m
    y_k-m+1^(k),


where


    y_k^(m) = ∫_0^∞e^-ζ x(ζ x)^m/m! F_Y_k(x),
      
    m = 0,1,…, k= 0,1,….


By definition, it is readily verified that y_k^(m) satisfies
the recursion (<ref>) and (<ref>).
We thus obtain (<ref>) from (<ref>).




§ PROOF OF LEMMA <REF>



From (<ref>), we have


    R^*(ζ - ζ z)
       =
    ∫_x=0^∞ e^-(ζ - ζ z)x x
    ∫_y=x^∞D̂(y) exp[Q(y-x)]
    
       =
    ∫_y=0^∞D̂(y) 
    ∫_x=0^y
    exp[Q(y-x)]
    e^-(ζ - ζ z)x x
    
       =
    ∫_y=0^∞D̂(y) 
    ∫_x=0^y
    ∑_k=0^∞e^-ζ (y-x){ζ (y-x)}^k/k!·
    (I+ζ^-1Q)^k
    e^-(ζ - ζ z)x x
    
       =
    ∑_k=0^∞∫_y=0^∞ 
    e^-ζ yD̂(y) 
    (I+ζ^-1Q)^k
    ∫_x=0^y
    {ζ (y-x)}^k/k!
    e^ζ zx x
    
       =
    ∑_k=0^∞∫_y=0^∞ 
    e^-(ζ - ζ z) yD̂(y) 
    (I+ζ^-1Q)^k
    ∫_x=0^y
    {ζ (y-x)}^k/k!
    e^-ζ z(y-x) x
    
       =
    ∑_k=0^∞
    z^-(k+1)ζ^-1∫_y=0^∞ 
    e^-(ζ - ζ z) yD̂(y) 
    (I+ζ^-1Q)^k
    ∫_x=0^y
    e^-ζ z(y-x){ζ z(y-x)}^k /k!·ζ z
     x
    
       =
    ∑_k=0^∞∑_m=k+1^∞
    z^-(k+1)ζ^-1∫_y=0^∞ 
    e^-(ζ - ζ z) y·e^-ζ zy (ζ zy)^m /m!D̂(y) 
    (I+ζ^-1Q)^k
    
       =
    ∑_k=0^∞∑_m=k+1^∞
    z^m-k-1ζ^-1∫_y=0^∞e^-ζ y(ζ y)^m /m!D̂(y) 
    (I+ζ^-1Q)^k
    
       =
    ∑_k=0^∞∑_m=k+1^∞
    z^m-k-1ζ^-1D̂^(m)
    (I+ζ^-1Q)^k
    
       =
    ∑_k=0^∞∑_m=0^∞
    z^mζ^-1D̂^(m+k+1)
    (I+ζ^-1Q)^k
    
       =
    ∑_m=0^∞
    z^m
    ζ^-1∑_k=0^∞D̂^(m+k+1)
    (I+ζ^-1Q)^k,


Therefore, we have proved


    R^*(ζ-ζ z)
    =
    ∑_m=0^∞R^(m) z^m.



Note that we have from (<ref>), 


    v̂(ζ-ζ z)
       =
    (1-ρ-ρ_)κ
    +
    v̂(ζ-ζ z)
    R^*(ζ-ζ z),


i.e.,


    ∑_m=0^∞v̂^(m) z^m
       =
    (1-ρ-ρ_)κ
    +
    ∑_m=0^∞
    z^m
    ∑_k=0^m
    v̂^(k)R^(m-k).


Comparing the coefficients of both sides of this equation and
re-arranging terms, we have established the claim of Lemma <ref>.



99
Kaul2011
S. Kaul, M. Gruteser, V. Rai, and J. Kenney,
“Minimizing Age of Information in Vehicular Networks,”
in  Proc. IEEE SECON 2011, pp. 350–358, 2011.

Kaul2011-2 S. Kaul, R. Yates, and M. Gruteser, 
“On Piggybacking in Vehicular Networks,” 
in Proc. of IEEE GLOBECOM 2011, 2011.

Kaul2012
S. Kaul, R. Yates, and M. Gruteser, 
“Real-Time Status: How Often Should One Update?” 
in Proc. IEEE INFOCOM 2012, pp. 2731–2735, 2012.

Yates2021
R. Yates, Y. Sun,  D. R. Brown, S. K. Kaul, E. Modiano, and S. Ulukus, 
“Age of Information: An Introduction and Survey,” 
IEEE J. Sel. Areas Commun., vol. 39, no. 5, pp. 1183–1210, 2021. 

Costa2016 M. Costa, M. Codreanu, and A. Ephremides, 
“On the Age of Information in Status Update Systems with Packet Management,” 
IEEE Trans. Inf. Theory, vol. 62, no. 4, pp. 1897–1910, 2016.

Inoue2019 Y. Inoue, H. Masuyama, T. Takine, and T. Tanaka, 
“A General Formula for the Stationary Distribution of the Age of 
Information and Its Application to Single-server Queues,”
IEEE Trans. Inf. Theory, vol. 65, no. 12, pp. 8305–8324, 2019.

Champati2019 J. P. Champati, H. Al-Zubaidy, and J. Gross, 
“On the Distribution of AoI for the GI/GI/1/1 and GI/GI/1/2* Systems:
Exact Expressions and Bounds,” 
in Proc. of IEEE INFOCOM 2019, pp. 37–45, 2019.

Yates2019
R. Yates, S. K. Kaul,
“The Age of Information: Real-Time Status Updating by Multiple Sources”,
IEEE Trans. Inf. Theory, vol. 65, no. 3,  pp.  1807-1827, 2019.

Moltafet2020
M. Moltafet, M. Leinonen, M. Codreanu, “On the Age of Information in Multi-Source Queueing Models," IEEE Trans. Comm., vol. 68, no. 8, pp. 5003-5017, 2020.

Inoue2024
Y. Inoue and T. Takine, “Exact Analysis of the Age of Information in the Multi-Source M/GI/1 Queueing System,”
arXiv:2404.05167, 2024.

Chen2022
Z. Chen, D. Deng, C. She, Y. Jia, L. Liang, S. Fang, M. Wang, Y. Li,
“Age of Information: The Multi-Stream M/G/1/1 Non-Preemptive System," IEEE Trans. Comm., vol. 70, no. 4, pp. 2328-2341 2022.

Jiang2021
Y. Jiang, N. Miyoshi, “Joint Performance Analysis of Ages of Information in a Multi-Source Pushout Server," IEEE Trans. Inf. Theory, vol. 68, no. 2, pp. 965-975,  2021.

Liu2021
Z. Liu, Y. Sang, B. Li, B. Ji, “A Worst-Case Approximate Analysis of Peak Age-of-Information via Robust Queueing Approach," IEEE INFOCOM 2021 - IEEE Conference on Computer Communications, Vancouver, BC, Canada,  pp. 1-10, 2021.

Cinlar1972
E. Çinlar, “Superposition of Point Processes,” In: Stochastic Point Processes, ed. P. A. W. Lewis, Wiley, New York, pp. 546–606, 1972.

Asmussen2003
S. Asmussen,
Applied Probability and Queues, 2nd ed.,
Springer, New York, 2003.

Abate1995 
J. Abate, G. Choudhury, and W. Whitt, “Calculating the M/G/1 Busy-Period Density and LIFO Waiting-Time Distribution by Direct Numerical Transform Inversion,” Oper. Res. Lett., vol. 18, no. 3, pp. 113-–119, 1995.

Takagi1991 H. Takagi,
Queueing Analysis: Vacation and Priority Systems, Part 1,
Elsevier Science Publishers, Amsterdam, 1991.

Hooke1972
J. Hooke,
“A Priority Queue with Low-Priority Arrivals General,”
Oper. Res., vol. 20, no. 2, pp. 373–380, 1972.

Ott1984
T. Ott, 
“On the M/G/1 Queue with Additional Inputs,”
J. Appl. Prob., vol. 21, no. 1, pp. 129–142, 1984.

Cohen1982
J. W. Cohen, 
The Single-Server Queue, North-Holland, Amsterdam, 1982.

Abate1992 
J. Abate and W. Whitt, 
“The Fourier-Series Method for Inverting Transforms of Probability Distributions,” 
Queueing Syst., vol. 10, nos. 1-2, pp. 5-–87, 1992.

Shaked2007
M. Shaked and J. G. Shanthikumar, Stochastic Orders, Springer, New York, 2007.












Gurland1967
J. Gurland, 
“An Inequality Satisfied by the Expectation of the Reciprocal
of a Random Variable,”
Am. Stat., 
vol. 21, no. 2, pp. 24–25, 1967.

Behboodian1994
J. Behboodian,
“Covariance Inequality and Its Applications,”
Int. J. Math. Educ. Sci. Technol.,
vol. 25, no. 5, pp. 643–647, 1994.

Lehmann1966
E. L. Lehmann, 
“Some Concepts of Dependence,” 
Ann. Math. Statist., vol. 37, no. 5, pp. 1137–1153, 1966.

Daley1977
D. J. Daley, “Inequalities for Moments of Tails of Random Variables, 
with a Queueing Application,” Z. 
Wahrscheinlichkeit, vol. 41, no. 2,
pp. 139–143, 1977.

Takine2002
T. Takine, 
“Matrix Product-Form solution for an LCFS-PR Single-Server Queue with
Multiple Arrival Streams Governed by a Markov Chain.
Queueing Syst., vol. 42, no. 2, pp. 131–151, 2002.

Takine1994
T. Takine and T. Hasegawa, 
“The Workload in the MAP/G/1 Queue with State-Dependent Services: Its
Application to a Queue with Preemptive Resume Priority,”
Stoch. Models, vol. 10, no. 1, pp. 183–204, 1994.

Tijms1994
H. C. Tijms, 
Stochastic Models, An Algorithmic Approach, Wiley, Chichester, 1994.

Door1988
E. A. van Doorn and G. J. K. Regterschot, 
“Conditional PASTA,” Oper. Res. Lett., vol. 7, pp. 229–232, 1988.

Davies1979
B. Davies and B. Martin,
“Numerical Inversion of the Laplace Transform: A Survey and Comparison of Methods,”
J. Comput. Phys.  vol. 33, no. 1, pp. 1–32, 1979.

Roberts1996
J. Roberts, U. Mocci, and J. Virtamo (eds.), Broadband Network Teletraffic, Springer, New York, 1996.






