*State vector at next time slot* = ⎧ ⎪⎨ ⎪⎩ (*ni* + 1, *mi* + 1) no packets of higher priorities arrive, the service is not over, (*ni* + 1, 0) one packet with higher priority comes, (*mi* + 1, 0) no packets with higher priorities arrive, the service is over. (8)

> when packet service preemption is allowable. After all the state transfers are described and their transition probabilities are determined, we can obtain the stationary equations, which can be solved completely as in [58] or by using the PGF method as in this paper. Like [58], the service time distribution in this case can be arbitrary.

> Although there are multiple sources, it can be seen that the age analysis of each source is easy when the status updating system has no buffer. Notice that in this case, no queue is formed before the server; thus, there is no chance that the packets from different sources are combined. As a result, the packets from every source are totally divided, and the AoI of each source can be analyzed separately. The situation is much more difficult if the system has a buffer. As an example, we consider the AoI of each source of a two-source system, which is depicted in the second picture of Figure 1.

> We can define a six-dimensional state vector (*n*1, *n*2, *m*1, *m*2, *l*1, *l*2) to describe the AoI of two sources simultaneously, where the state components represent the values of two AoIs at the destination, the age of a packet in the server, and the age of a packet in the system's buffer. In every position of the system, apart from the "age", it is necessary to indicate which source the packet comes from. Therefore, a three-dimensional state vector (*n*, *m*, *l*) that does not include this information is not sufficient. Notice that at any time, at most one of *m*1 and *m*2 are non-zero. This is the same for the parameters *l*1 and *l*2. When there is a buffer in front of the server, apparently a queue is formed if a packet arrives and finds that the server is currently busy. Each one of two packets in the system (one is in the server and the other is in the buffer) may come from source *s*1 or *s*2. Of course, these two packets may belong to different sources. Although the problem becomes complex, theoretically, all the state transfers of every initial six-dimensional state vector can be determined explicitly, since the randomness that causes the state vector transfers are limited to random packet arrival, the service of the packet, and the additional packet preemption. Then, according to the balance of probabilities in the steady state, the stationary equations are established; this solves the first half of the AoI analysis. Details of the latter half—that is, deriving the

average AoI from the group of stationary equations—can be found in the procedures in this article.

We find that in [61], the authors obtained the average continuous AoI for the same twosource status updating system in Figure 1b using the SHS method. They added another assumption that the packet in the server and the packet in the buffer must belong to different sources in their second and the third considered situation and named the policies "source-aware packet management". As we have mentioned above, although the packets from two sources are still combined, after adding this restriction, the complexity of the problem has been greatly reduced.

In fact, the state vector defined for a discrete time system has a very clear physical meaning. For the status updating system with the FCFS queue, the first parameter denotes the AoI and the other state components represent the ages of packets in the server and in the buffer of the system. Thus, a (*c* + 1)-dimensional state vector is needed if the size of the system is equal to *c*. Compared with analyzing the AoI of discrete systems, in the SHS method, the defined state vector is sometimes easier, such as in the system with multi-sources. In [61], in order to characterize the AoI of one source in a two-source system, the authors used a four-dimensional state vector [*x*0(*t*), *x*1(*t*), *x*2(*t*), *x*3(*t*)] that describes the evolutions of AoI when different random events occur. As mentioned before, we use the six-dimensional state vector (*n*1, *n*2, *m*1, *m*2, *l*1, *l*2) describing the random changes of both source 1 and source 2. The parameter *n*1 or *n*2 can also be deleted if only one of two sources are analyzed. In our proposed method, we show the correspondence between the dimension of the state vector and the size of the discrete system; this may not be a unique way to define the discrete state vector. Although considering the AoI of the discrete time system has higher computational complexity, the biggest advantage of discrete AoI analysis is that it can obtain the stationary distribution of the AoI.

Except the simple status updating networks given in Figure 1, we have also obtained the average discrete AoI for a status updating system with two-stage service, where for simplicity, in front of each server, no buffer is equipped. For the system with two parallel servers, the age analysis is more difficult, since some packets may become "ineffective" if one packet is generated later but arrives to the destination earlier. Some policies need to be identified to deal with these packets—for instance, deleting the packet directly once it becomes ineffective. If nothing is done, when an ineffective packet is obtained at the receiver, the value of AoI will not be reduced.

Another direction of generalization we shall talk about is the distribution of packet service time (while the packet arrival process is still Bernoulli). Now, taking the size 2 status updating system as an example, we explain how the service time distribution is relaxed to be an arbitrary distribution. Using a three-dimensional state vector (*n*, *m*, *l*), we can fully describe the random dynamics including the AoI at the receiver and the age of two packets in the system if both the packet interarrival time and the service time have memoryless properties. In each time slot, the changes of the AoI's value and the packet ages depend on random packet arrival, which is memoryless and independent, and whether the packet service is over. When the service time distribution is arbitrary, the probability that the service is completed in one time slot is related to the time this packet has experienced in the server. Let *S* be the random variable of service time, and we represent the general distribution as

$$\Pr\{S=i\} = q\_i \qquad (i \ge 1) \tag{9}$$

We assume that, before the current time slot, the packet has stayed in server for *j* time slots; then, the probabilities that determine the state vector transfers should be the following two conditional probabilities:

$$\Pr\{S = j + 1 | S > j\} \text{ and } \Pr\{S > j + 1 | S > j\} \tag{10}$$

Therefore, if we have knowledge about this passed service time *j*, as before, all the state transfers of state vector (*n*, *m*, *l*) can be completely described and the age analysis becomes feasible. Since no one of the three parameters *n*, *m*, and *l* can provide this information, it is natural to introduce an extra component, say *k*, to denote the service time that the packet has consumed and constitute the four-dimensional state vector (*n*, *m*, *l*, *k*). In this way, the possible state transfers of this four-dimensional state vector can be described and the transition probabilities can also be determined. For example, let the initial state vector be (*n*, *m*, *l*, *k*)—we have the state transfers and transition probabilities as

$$\text{Next state vector} = \begin{cases} (n+1, m+1, l+1, k+1) & \text{the service is not over with prob. } 1 - q\_{k+1} / \sum\_{i=k+1}^{\infty} q\_i \\ (m+1, l+1, 0, 0) & \text{the service completes with prob. } q\_{k+1} / \sum\_{i=k+1}^{\infty} q\_i \end{cases} \tag{11}$$

where we assume the queue discipline is FCFS and there is no packet preemption. We show that the four parameters *n*, *m*, *l*, and *k* satisfy the relationships *n* > *m* > *l* ≥ 0 and *n* > *m* ≥ *k* ≥ 0. The first one holds because *n*, *m*, and *l* are three ages of packets generated in chronological order, and *n* > *m* ≥ *k* is satisfied since the packet system time *m* must be larger than or equal to the service time of the packet, which is denoted by *k*. These relations determine which vectors are qualified state vectors. Although we show that the state transfers can be analyzed and the group of stationary equations can be determined by considering the balance of those stationary probabilities; however, it can be expected that solving these equations is not easy. Since the service time probabilities *qi*s are arbitrary, the expression of the average AoI, as we can determine in later work, will not be closed-formed. It is also important to note that the PGF method cannot be used when the service time distribution is not geometric, because the transition probabilities is no longer the same for different state vectors and thus cannot be the common factor.

Summarizing the above discussions, we have proved that on the basis of original memoryless status updating system, by introducing an extra component to denote the time the packet has consumed in the server, the age analysis becomes feasible for the situation where the packet service time is arbitrarily distributed. Although it may be difficult to obtain the expressions of the system's average AoI, the idea is still applicable when we generalize the size 2 system to a status updating system with arbitrary size *c*. In one of our works, we have shown that for a size *c* discrete time status updating system with Bernoulli packet arrivals and geometrically distributed service time, in order to fully characterize the real time transfers of the system's AoI and all the packet ages, a (*c* + 1) dimensional state vector (*n*, *m*1, ··· , *mc*) should be defined. By adding an extra state component *k* that records the service time the packet has experienced in the server, according to previous discussions, the age analysis can be generalized to a size *c* status updating system whose service time distribution is arbitrary (at least we can establish all the stationary equations).

We have to attribute the above idea to [62], in which the authors considered the timely transmission of the updates over an erasure channel. They assume that each update consists of *k* symbols and the symbol erasure in each time slot is an i.i.d. Bernoulli process. The aim of [62] is to design an optimal online transmission scheme to minimize the time average AoI, and the problem is formulated as a Markov Decision Process (MDP). Although the optimization of AoI is not our interest, the state tuple (*δt*, *dt*, *lt*) defined in Section 2. A is very enlightening, based on which the transmission policy at the next time slot is determined. At the *t*-th time slot, the notation *δt* denotes the value of AoI, *dt* is the age of the next update, i.e., the packet at the head of the queue, and *lt* records the number of symbols that has been obtained successfully up to this time slot—these symbols belong to the update that is transmitted currently. A similar timely source coding problem was also discussed in [63], in which the authors also pointed out that the length of the encoded update is equivalent to the service time of the update, and their considered system behaves as a discrete time *Geo*/*G*/1 queue (we use the notation *Ber*/*G*/1). Therefore, the role of *lt* in [62] can be regarded (or redefined) as the service time that the current update has consumed. By adding this knowledge, the distribution of the source in these papers and the service time distribution in the discrete time status updating system which we study in this part can be arbitrary.

In previous paragraphs, we explain the idea and methods used to study the AoI of discrete time status updating systems. We have shown how the discrete AoI is characterized for the basic system, the system with multiple sources, and the system whose service time distribution is arbitrary. As part of AoI theory, we believe that discrete AoI deserves more attention, and it is meaningful to establish analytical results including the AoI's mean and its distribution for more general systems. In particular, the proposed possible relationship in (4) shows that discussing discrete AoI not only has independent theoretical significance but also helps to determine certain results about continuous AoI. If one problem is difficult in the continuous time model, it is a choice to consider it in discrete time settings.

#### *1.4. The Work in the Current Paper*

We have discussed numerous topics of discrete AoI in the previous subsection, and it is inappropriate to consider all the issues in one article. In this paper, we focus on the stationary AoI of a discrete time system with a *Ber*/*Geo*/1/2 and *Ber*/*Geo*/1/2∗ queue and discuss both in a single model. We assume the packet in the buffer can be probabilistically preempted by the fresher packets from the source and define the queue model in this scenario as *Ber*/*Geo*/1/2∗/*η*, where *η* is the preemption probability. In the literature of AoI, the probabilistic packet preemption (replacement) has been studied in [64]. In [65], the probabilistic preemption was considered in the scenario where a CPU is used frequently to deal with the unpredictable tasks. Then, for the case of *η* = 0, the queue model of the system reduces to *Ber*/*Geo*/1/2, while when *η* is equal to 1, the status updating system with *Ber*/*Geo*/1/2∗ queue is obtained. For the general case, we derive the explicit expression of the system's average AoI. By writing the defined PGF as the power series, for two extreme cases of *η* = 0 and *η* = 1, the distribution expressions of two discrete AoIs are determined as well.

The rest of the paper is organized as follows. In Section 2, we describe the model of a discrete time status updating system with probabilistic packet preemption. The stationary distribution and the mean of the system's AoI are also defined. By analyzing the steady state of a three-dimensional stochastic age process, in Section 3, we obtain the explicit formula of the average AoI under general preemption probability using the probability generation function (PGF) method. In Section 4, let *η* = 0 and *η* = 1, and we determine the average AoIs Δ*Ber*/*Geo*/1/2 and Δ*Ber*/*Geo*/1/2∗ from the general expression derived previously in Section 3. Furthermore, in order to obtain the stationary distribution of two discrete AoIs, we write the PGF as power series. Then, the coefficient before *xn* gives the probability that the AoI takes value *n* for each *n* ≥ 1. Numerical results are placed in Section 5. For the general case, we illustrate the relationships between the average AoI and *η* and the traffic intensity *ρd*, respectively. In addition, the mean and the cumulative probabilities of three discrete AoIs including Δ*Ber*/*Geo*/1/1, Δ*Ber*/*Geo*/1/2, and Δ*Ber*/*Geo*/1/2∗ are depicted. These average discrete AoIs and their corresponding average continuous AoIs are also numerically compared in Section 5. Finally, we conclude the paper in Section 6.

#### **2. System Model and Problem Formulation**

We depict the model of the status updating system which uses the *Ber*/*Geo*/1/2∗/*η* queue in Figure 2, in which the packet in the system's buffer can be preempted by a fresher packet from the source *s* with probability *η*. The packet arrivals to the transmitter are assumed to form a Bernoulli stochastic process—that is, in each time slot, a new packet comes with an identical probability, which we denote by *p*. Packet service time follows the geometric distribution with intensity *γ*. The updated packet generated at *s* is transmitted to the destination *d* through the transmitter, in which a random period of time is consumed. The age of information (AoI) at *d* is defined as the time elapsed since the generation time of the last obtained packet. Within the time when no packet is received, the value of AoI increases by 1 after each time slot ends. Every time a packet passes the transmitter and arrives to *d*, the AoI will be reduced to the system time of the obtained packet, which is actually equal to the instantaneous age of this packet.

**Figure 2.** The model of the discrete time status updating system with probabilistic packet preemption in the system's buffer.

Let *a*(*k*) be the value of AoI in the *k*th time slot. The AoI at the next time slot, *a*(*k* + 1), is determined by

$$a(k+1) = \begin{cases} a(k) + 1 & \text{if no packet is obtained,} \\ a(k) + 1 - \mathcal{Y}\_{\bar{\jmath}} & \text{when } \bar{\jmath}\text{th generated packet arrives to } d. \end{cases} \tag{12}$$

where *Yj* is the interarrival time between the (*j* − 1)th and *j*th arriving packet.

Notice that these (*j* − 1)th and *j*th packets may be generated discontinuously, since between them, some updating packets may be discarded when they arrive and find the system full. Actually, this is exactly the difference between the finite and infinite status updating systems. Based on this observation, in [59], we have determined the average AoI and its stationary distribution for an infinite size status updating system with Bernoulli packet arrivals and geometric service time.

Denote the stationary AoI for the system with probabilistic packet preemption as Δ*Ber*/*Geo*/1/2∗/*η*. We define the time average AoI as follows, which is equal to the mean of the AoI because the age process is assumed to be ergodic. We have

$$\overline{\Delta}\_{\text{Ber}/\text{Geo}/1/2^\*/\eta} = \lim\_{T \to \infty} \frac{1}{T} \sum\_{k=1}^T a(k) \tag{13}$$

$$\hat{\mathbf{x}} = \lim\_{T \to \infty} \frac{1}{T} \sum\_{i=1}^{M\_T} i \cdot |\{1 \le k \le T : a(k) = i\}| \tag{14}$$

$$=\sum\_{i=1}^{\infty} i \cdot \pi\_i \tag{15}$$

where |{1 ≤ *k* ≤ *T* : *a*(*k*) = *i*}| is the times that the AoI takes value *i*, and *MT* = max1≤*k*≤*T a*(*k*) is the maximal discrete AoI in *T* time slots. For each *i* ≥ 1,

$$\pi\_i = \lim\_{T \to \infty} \frac{|\{1 \le k \le T : a(k) = i\}|}{T} \tag{16}$$

is the probability that the stationary AoI takes value *i*. In fact, the probability distribution {*πi*, *i* ≥ 1} forms the stationary distribution of the AoI Δ*Ber*/*Geo*/1/2∗/*η*.

The randomness of both packet arrivals and the service time in the server, along with the probabilistic packet preemption in the system's buffer, together make the AoI at the destination change randomly. After one time slot, the value of AoI may increase by 1 if no packet is obtained or drop to the age of the obtained packet at that time if one such packet is successfully received. In order to fully describe these random dynamics of AoI, we propose to use a three-dimensional state vector to simultaneously record the changes of the AoI, the age of a packet in the server, and the age of the packet waiting in the buffer, and then constitute the three-dimensional stochastic process. Next, the steady-state of this multipledimensional age process is analyzed. To obtain the mean and the distribution of AoI, we define the PGF corresponding to the stationary distribution of the three-dimensional age

process, from which both the AoI's mean and its distribution can be obtained. The detailed analysis of the system's AoI is given in Section 3.

#### **3. AoI Analysis for Status Updating System with Probabilistic Packet Preemption**

Define the three-dimensional state vector (*n*, *m*, *l*), where we use *n* to denote the AoI at destination *d*, and the other two parameters *m* and *l* are the ages of the packets in the system's server and the buffer. In the *k*th time slot, if the server is busy while the buffer is empty, then *nk* and *mk* are greater than 0 but *lk* = 0. When both the server and the buffer are empty, we have *mk* = *lk* = 0. In this case, the entire system is empty.

Consider the following three-dimensional age process

$$AoI\_{PP} = \left\{ (n\_k, m\_k, l\_k) : n\_k > m\_k \ge l\_k \ge 0, k \in \mathbb{N} \right\} \tag{17}$$

where the subscript "PP" in expression (17) is the abbreviation of probabilistic preemption. Notice that when the system is empty, the last two parameters *mk* and *lk* are both equal to 0. When there are two packets in the system, i.e., one is in the server and the other is in the buffer, we show that the state components satisfy *nk* > *mk* > *lk* ≥ 1, since in a path from the source to the receiver, the packet ahead always has a greater age. It is clearly shown later that this relationship facilitates the derivation of probability generation function *HPP*(*x*), which is defined in Equation (20).

Define three random variables *A*, *B*, and *F* to represent whether a packet is generated in a time slot, if the service of the packet is completed, and if the arriving packet replaces the original one in the buffer. For each possible initial state vector, according to different realizations of r.v.s (*A*, *B*, *F*), the state transfers of the three-dimensional state vector (*n*, *m*, *l*) can be described specifically. We list all of them using Table 2. For example, the third row of the table shows that a packet of age *l* is in the buffer and a new packet arrives, since the r.v. *A* takes value 1. However, *F* = 0 means that this new packet will not substitute the original one; meanwhile, *B* = 0 implies that the packet service is not over at this time slot. Summarizing all these events, the beginning state vector (*n*, *m*, *l*) will transfer to (*n* + 1, *m* + 1, *l* + 1) at the next time slot, and the transition probability is determined as *p*(1 − *γ*)(1 − *η*). The other cases in the third column of Table 2 are obtained through similar discussions.

From the state transfers given in Table 2 and the corresponding transition probabilities, we can establish all the stationary equations that characterize the steady-state of age process *AoIPP*. Let *π*(*n*,*m*,*l*), *n* > *m* ≥ *l* ≥ 0 be the probability that the process stays at the state vector (*n*, *m*, *l*) when it reaches the steady-state; we show that these stationary probabilities *π*(*n*,*m*,*l*) satisfy the following equations.

$$\begin{cases} \pi\_{(n,m,l)} = \pi\_{(n-1,m-1,l-1)}[(1-p)(1-\gamma)+p(1-\gamma)(1-\eta)] & (n>m>l\geq 2) \\ \pi\_{(n,m,1)} = \pi\_{(n-1,m-1,0)}p(1-\gamma) + \sum\_{j=1}^{m-2}\pi\_{(n-1,m-1,j)}p(1-\gamma)\eta & (n>m\geq 3) \\ \pi\_{(n,2,1)} = \pi\_{(n-1,0,0)}p(1-\gamma) & (n\geq 3) \\ \pi\_{(n,m,0)} = \pi\_{(n-1,m-1,0)}(1-p)(1-\gamma) & (n>m\geq 2) \\ \qquad + \sum\_{k=n}^{\infty}\pi\_{(k,n-1,m-1)}[(1-p)\gamma+p\gamma(1-\eta)] & (n>m\geq 2) \\ \pi\_{(n,1,0)} = \pi\_{(n-1,0,0)}p\gamma(1-\gamma) & \\ \qquad + \sum\_{k=n}^{\infty}\pi\_{(k,n-1,0)}p\gamma+\sum\_{k=n}^{\infty}\pi\_{(k,n-1,j)}p\gamma\eta & (n\geq 3) \\ \pi\_{(2,1,0)} = \pi\_{(1,0,0)}p(1-\gamma) + \sum\_{k=n}^{\infty}\pi\_{(k,n-1,0)}p\gamma & (n\geq 2) \\ \pi\_{(n,0,0)} = \pi\_{(n-1,0,0)}(1-p) + \sum\_{k=n}^{\infty}\pi\_{(k,n-1,0)}(1-p)\gamma & (n\geq 2) \\ \pi\_{(1,0,0)} = \sum\_{n=1}^{\infty}\pi\_{(n,0,0)}p\gamma & (n\geq 2) \end{cases} \tag{18}$$

**Table 2.** The state vector transfers of age process *AoIPP*.

We explain the stationary equations only for a part of the state vectors and show that the other equations in (18) can be determined in a similar manner. Firstly, for the fifth row of (18), the state vector (*n*, 1, 0) can be obtained from (*n* − 1, 0, 0) assuming that a new packet arrives and enters the server directly, but the service does not end in a single time slot. Next, from the current state vector (*k*, *n* − 1, 0), *k* ≥ *n*, if the service of the age (*n* − 1) packet is completed and a new packet arrives at the same time slot, it is observed that the packet of age (*n* − 1) will be sent to the receiver, which makes the AoI change to *n* at next time slot. The new packet enters the server; thus, the middle parameter of the state vector changes to 1. This gives the expected state (*n*, 1, 0). Since in this case, the buffer is empty, when a new packet comes, it occupies the buffer directly, and no packet preemption occurs. At last, we consider the situation where the age process begins with an arbitrary state (*k*, *n* − 1, *j*) where *k* > *n* − 1 > *j* ≥ 1. As long as the packet service is completed and at the same time a new packet arrives preempting the original one in the buffer, again, we will obtain the state vector (*n*, 1, 0) after one time slot. Combining all of the above cases, the stationary equation corresponding to (*n*, 1, 0) is finally determined. In addition to the fifth row, we also explain the last equation in (18). Observing that in order to obtain the state vector (1, 0, 0), the receiver needs a packet of age 1, the system has then to be emptied. This state can be transferred to only from (*n*, 0, 0), and the service time of the newly arrived packet is restricted to be only one time slot.

To derive the expression of the average AoI Δ*Ber*/*Geo*/1/2∗/*η*, we do not solve Equation (18) although this approach is feasible for the AoI analysis of tje current system. In our work [60], we analyzed the AoI of a status updating system with *Ber*/*Geo*/1/1, *Ber*/*Geo*/1/2, and *Ber*/*Geo*/1/2∗ queues, and the expression of the AoI's stationary distribution was determined for each case. There, we completely solved the stationary equations for each system and obtained the explicit expression for every stationary probability. Notice that this work can be regarded as a discrete correspondence of the packet management of continuous AoI in [8,9]. Assuming all the probabilities *π*(*n*,*m*,*l*) have been determined by solving Equation (18), we have

$$\Pr\left\{\Delta\_{\text{Ber}/\text{Geo}/1/2^{\ast}/\eta} = n\right\} = \begin{cases} \pi\_{(1,0,0)} & (n=1) \\ \pi\_{(n,0,0)} + \sum\_{l=0}^{n-2} \sum\_{m=l+1}^{n-1} \pi\_{(n,m,l)} & (n\ge 2) \end{cases} \tag{19}$$

since the probability that the AoI takes each *n* is equal to the sum of all the stationary probabilities with the identical first component. Equation (19) gives the stationary distribution of the AoI, from which we can calculate the average value of AoI as

$$\Delta\_{\text{Ber}/\text{Geo}/1/2^\*/\eta} = \sum\_{n=1}^{\infty} n \cdot \Pr\left\{ \Delta\_{\text{Ber}/\text{Geo}/1/2^\*/\eta} = n \right\}$$

However, the number of calculations to solve Equation (18) may be large, and apart from this, extra computations are required to determine the AoI's distribution according to Formula (19). Since the AoI is denoted by the first component, to obtain the distribution of AoI, we need to sum all the other state components. Notice that when the dimension of defined state vector is bigger, more calculations are required to determine the AoI's distribution. Therefore, we must determine the mean of AoI and its distribution in another way, i.e., the probability generation function (PGF) method.

For 0 < *x* ≤ 1, define the probability generation function

$$H\_{PP}(\mathbf{x}) = \sum\_{n=1}^{\infty} \mathbf{x}^n \Pr\left\{ \Delta\_{\mathrm{Ber}/\mathrm{Geo}/1/2\*} \,|\, \eta = n \right\} \tag{20}$$

and we write *HPP*(*x*) further as

$$\begin{split} H\_{PP}(\mathbf{x}) &= \mathbf{x} \Pr\{\Delta\_{\mathrm{Ber}/\mathrm{Geo}\prime} \gamma\_{1/2^{\ast}} = 1\} + \sum\_{n=2}^{\infty} \mathbf{x}^{n} \Pr\{\Delta\_{\mathrm{Ber}/\mathrm{Geo}\prime} \gamma\_{1/2^{\ast}} = n\} \\ &= \mathbf{x} \pi\_{(1,0,0)} + \sum\_{n=2}^{\infty} \mathbf{x}^{n} \left\{ \pi\_{(n,0,0)} + \sum\_{l=0}^{n-2} \sum\_{m=l+1}^{n-1} \pi\_{(n,m,l)} \right\} \end{split} \tag{21}$$

$$\mathbf{x}^{n} = \sum\_{m=1}^{\infty} \mathbf{x}^{n} \pi\_{(n,0,0)} + \sum\_{m=2}^{\infty} \mathbf{x}^{n} \sum\_{d=0}^{n-2} \sum\_{m=l+1}^{n-1} \pi\_{(n,m,l)} \tag{22}$$

$$\mathbf{x} = \sum\_{n=1}^{m-1} \mathbf{x}^n \pi\_{(n,0,0)} + \sum\_{l=0}^{m-2} \sum\_{m=l+1}^{m} \sum\_{n=m+1}^{m} \mathbf{x}^n \pi\_{(n,m,l)} \tag{23}$$

$$\mathbf{x} = \sum\_{m=1}^{\infty} \mathbf{x}^{n} \pi\_{(n,0,0)} + \sum\_{m=1}^{\infty} \sum\_{n=m+1}^{\infty} \mathbf{x}^{n} \pi\_{(n,m,0)} + \sum\_{l=1}^{\infty} \sum\_{m=l+1}^{\infty} \sum\_{n=m+1}^{\infty} \mathbf{x}^{n} \pi\_{(n,m,l)} \tag{24}$$

where in (21) we have used the probability expressions (19). Equation (23) is obtained by exchanging the summation order in (22). In Equation (24), we divide the PGF *HPP*(*x*) into three parts. It can be seen in the following paragraphs that the entire function (20) is obtained by determining these parts separately.

According to expression (20), immediately, we have

$$H\_{PP}(1) = 1, \quad \left. \frac{\mathrm{d}H\_{PP}(\mathbf{x})}{\mathrm{d}\mathbf{x}} \right|\_{\mathbf{x}=1} = \overline{\Lambda}\_{\mathrm{Ber}/\mathrm{Geo}/1/2^\*/\eta} \tag{25}$$

That is, the average AoI can be obtained from the PGF's derivative at point *x* = 1, and the probability that the steady state AoI equals *n* is determined by the coefficient before the term *xn* for every *n* ≥ 1.

Now, we determine the PGF *HPP*(*x*). For 0 < *x* ≤ 1, define the functions

$$\begin{aligned} h\_1(\mathbf{x}) &= \sum\_{n=1}^{\infty} \mathbf{x}^n \pi\_{(n,0,0)} \\ h\_2(\mathbf{x}) &= \sum\_{m=1}^{\infty} \sum\_{n=m+1}^{\infty} \mathbf{x}^n \pi\_{(n,m,0)} \\ h\_3(\mathbf{x}) &= \sum\_{l=1}^{\infty} \sum\_{m=l+1}^{\infty} \sum\_{n=m+1}^{\infty} \mathbf{x}^n \pi\_{(n,m,l)} \end{aligned}$$

and

$$h\_2^{(m)}(\mathbf{x}) = \sum\_{m=1}^{\infty} \sum\_{n=m+1}^{\infty} \mathbf{x}^m \pi\_{(n,m,0)}$$

We first give the following lemma, from which the PGF *HPP*(*x*) can be determined completely.

**Lemma 1.** *For the functions hi*(*x*)*,* 1 ≤ *i* ≤ 3 *and h*(*m*) 2 (*x*)*, we have*

$$h\_1(\mathbf{x}) = \frac{p\gamma M\_1 \mathbf{x}}{1 - (1 - p)\mathbf{x}} + \frac{(1 - p)\gamma \mathbf{x}}{1 - (1 - p)\mathbf{x}} h\_2^{(m)}(\mathbf{x}) \tag{26}$$

$$h\_2(\mathbf{x}) = \frac{p(1-\gamma)\mathbf{x}}{1 - (1-p)(1-\gamma)\mathbf{x}}h\_1(\mathbf{x}) + \frac{p\gamma\mathbf{x}}{[1-(1-p)(1-\gamma)\mathbf{x}][1-(1-\gamma)\mathbf{x}]}h\_2^{(m)}(\mathbf{x})\tag{27}$$

$$h\_3(\mathbf{x}) = \frac{p(1-\gamma)\mathbf{x}}{1-(1-\gamma)\mathbf{x}}h\_2(\mathbf{x})\tag{28}$$

*and it is determined that*

$$h\_2^{(m)}(\mathbf{x}) = \frac{\left[\gamma + p^2 (1 - \gamma)\eta - (1 - p)(1 - \gamma)(1 - p\eta)\gamma\mathbf{x}\right]M\_2\mathbf{x}}{\left[1 - (1 - p)(1 - \gamma)\mathbf{x}\right][1 - (1 - \gamma)(1 - p\eta)\mathbf{x}]} \tag{29}$$

*in which the numbers M*1 *and M*2 *are given as*

$$M\_1 = \frac{(1-p)\gamma^2}{(p+\gamma-2p\gamma)\gamma+p^2(1-\gamma)^2} \tag{30}$$

$$M2 = \frac{p\gamma(1-\gamma)}{(p+\gamma-2p\gamma)\gamma+p^2(1-\gamma)^2} \tag{31}$$

**Proof.** Lemma 1 is proved in Appendix A.

Using Lemma 1, we calculate the PGF *HPP*(*x*) as follows. Equation (24) shows

$$\begin{split} H\_{PP}(\mathbf{x}) &= h\_1(\mathbf{x}) + h\_2(\mathbf{x}) + h\_3(\mathbf{x}) \\ &= h\_1(\mathbf{x}) + h\_2(\mathbf{x}) + \frac{p(1-\gamma)\mathbf{x}}{1-(1-\gamma)\mathbf{x}} h\_2(\mathbf{x}) \\ &= h\_1(\mathbf{x}) + \frac{1-(1-p)(1-\gamma)\mathbf{x}}{1-(1-\gamma)\mathbf{x}} \left( \frac{p(1-\gamma)\mathbf{x}}{1-(1-p)(1-\gamma)\mathbf{x}} h\_1(\mathbf{x}) \right. \\ &\quad + \frac{p\gamma\mathbf{x}}{[1-(1-p)(1-\gamma)\mathbf{x}][1-(1-\gamma)\mathbf{x}]} h\_2^{(m)}(\mathbf{x}) \Bigg) \\ &= \frac{1-(1-p)(1-\gamma)\mathbf{x}}{1-(1-\gamma)\mathbf{x}} h\_1(\mathbf{x}) + \frac{p\gamma\mathbf{x}}{[1-(1-\gamma)\mathbf{x}]^2} h\_2^{(m)}(\mathbf{x}) \end{split} \tag{32}$$

where in (32) we have substituted Equation (27).

Using Equation (26) and merging the same terms, eventually, we obtain

$$\begin{split} H\_{PP}(\mathbf{x}) &= \frac{p\gamma M\_1 \mathbf{x} [1 - (1 - p)(1 - \gamma)\mathbf{x}]}{[1 - (1 - p)\mathbf{x}][1 - (1 - \gamma)\mathbf{x}]} \\ &+ \frac{\gamma \mathbf{x} \{1 - (1 - p)[2(1 - \gamma) + p\gamma]\mathbf{x} + (1 - p)^2(1 - \gamma)^2 \mathbf{x}^2\}}{[1 - (1 - p)\mathbf{x}][1 - (1 - \gamma)\mathbf{x}]^2} h\_2^{(m)}(\mathbf{x}) \end{split} \tag{33}$$

in which the function *h* (*m*) 2 (*x*) is given in Equation (29).

According to Formula (25), the average AoI of the system with probabilistic packet preemption is calculated in Theorem 1.

**Theorem 1.** *For the discrete time state updating system with a Ber*/*Geo*/1/2∗/*η queue, assuming the packet waiting in the buffer can be preempted by following fresher packets with probability η, then the average age of information of this system is determined as*

$$\overline{\Delta}\_{\text{Ber}/\text{Cov}/1/2^{\circ}\gamma} = \frac{(p+\gamma-p\gamma)(p+\gamma)-p\gamma}{p\gamma}M\_1 + \frac{(p+\gamma-p\gamma)^2 - p\gamma(1-p)}{p\gamma} \cdot \left. \frac{\text{dh}\_2^{(m)}(x)}{\text{dx}} \right|\_{x=1}$$

$$+\left. \frac{\left\{(p+\gamma-p\gamma)[1-3(1-p)(1-\gamma)] - 2p\gamma(1-p)\right\}p\gamma + \text{Poly}\_1}{p^2\gamma^2}M\_2 \tag{34}$$

*in which we define*

$$\text{Poly}\_1 = [(p + \gamma - p\gamma)^2 - p\gamma(1 - p)][2p(1 - \gamma) + (1 - p)\gamma] \tag{35}$$

*and the derivative of h*(*m*) 2 (*x*) *at point 1 is calculated as*

$$\left. \frac{\mathrm{d}h\_2^{(m)}(\mathbf{x})}{\mathrm{d}x} \right|\_{x=1} = \left( \frac{1}{p + \gamma - p\gamma} + \frac{p(1 - \gamma)(1 - p\eta)}{(p + \gamma - p\gamma)[\gamma + p(1 - \gamma)\eta]} \right) \mathrm{M}\_2 \tag{36}$$

*Let p* = *ρd* · *γ and substitute numbers M*1*, M*2*; the average AoI is also written as*

$$\begin{split} \overline{\Delta}\_{\text{Ber}/\text{Gou}/1/2^\*/\eta} &= \frac{1+\rho\_d(1-2\gamma)+\rho\_d^2(1-\gamma)^2+2\rho\_d^3(1-\gamma)^2}{\rho\_d\gamma[1+\rho\_d(1-2\gamma)+\rho\_d^2(1-\gamma)^2]} \\ &+ \frac{(1-\gamma)+\rho\_d(1-\gamma)(1-2\gamma)+\rho\_d^2(1-\gamma)(1-\gamma+\gamma^2)}{1+\rho\_d(1-2\gamma)+\rho\_d^2(1-\gamma)^2} \left\{ \frac{(2-\gamma)-\rho\_d\gamma(1-\gamma)}{\gamma[1+\rho\_d(1-\gamma)]} \\ &- \frac{(1-\gamma)-\rho\_d(1-\gamma)[1+\gamma-(1-\gamma)\eta]+\rho\_d^2\gamma^2(1-\gamma)\eta}{\gamma[1+\rho\_d(1-\gamma)(1+\eta)+\rho\_d^2(1-\gamma)^2\eta]} \right\} \end{split} \tag{37}$$

*where ρd* = *p*/*γ is defined as the discrete traffic load.*

**Proof.** The average AoI is determined by first computing the derivative of *HPP*(*x*) in (33) and then letting *x* = 1. Replacing parameter *p* with *ρd* · *γ*, expression (37) is obtained eventually. Although a certain amount of calculation is required, all the computations are straight-forward.

Here, we only provide the details from obtaining Equation (36). From (29), we have

$$\begin{split} \left. \frac{d\boldsymbol{h}\_{2}^{(m)}(\boldsymbol{x})}{d\boldsymbol{x}} \right|\_{\boldsymbol{x}=1} &= \left. \frac{\mathbf{d}}{d\mathbf{x}} \left[ \frac{\gamma + p^{2}(1-\gamma)\eta - (1-p)(1-\gamma)(1-p\eta)\gamma\mathbf{x}}{[1-(1-p)(1-\gamma)\mathbf{x}][1-(1-\gamma)(1-p\eta)\mathbf{x}]} \right]\_{\boldsymbol{x}=1} \right|\_{\boldsymbol{x}=1} \\ &= \left. \frac{\mathbf{d}}{d\mathbf{x}} \left( M\_{2} \frac{[\gamma + p^{2}(1-\gamma)\eta]\mathbf{x} - (1-p)(1-\gamma)(1-p\eta)\gamma\mathbf{x}^{2}}{1 - [(1-p)(1-\gamma) + (1-\gamma)(1-p\eta)]\mathbf{x} + (1-p)(1-\gamma)^{2}(1-p\eta)\mathbf{x}^{2}} \right) \right|\_{\boldsymbol{x}=1} \\ &= M\_{2} \frac{\text{Poly}\_{2}\cdot(p+\gamma-p\eta)[\gamma+p(1-\gamma)\eta] + [\gamma+p^{2}(1-\gamma)\eta - (1-p)(1-\gamma)(1-p\eta)\gamma] \cdot \text{Poly}\_{3}}{(p+\gamma-p\eta)^{2}[\gamma+p(1-\gamma)\eta]^{2}} \end{split} \tag{38}$$

where

$$\text{Poly}\_2 = \gamma + p^2 (1 - \gamma)\eta - 2(1 - p)(1 - \gamma)(1 - p\eta)\gamma$$

and

$$\begin{aligned} \text{Poly}\_3 &= (1-p)(1-\gamma) + (1-\gamma)(1-p\eta) - 2(1-p)(1-\gamma)^2(1-p\eta) \\ &= (1-p)(1-\gamma)[\gamma+p(1-\gamma)\eta] + (p+\gamma-p\gamma)(1-\gamma)(1-p\eta) \end{aligned}$$

Notice that

$$\begin{aligned} &\gamma + p^2(1-\gamma)\eta - (1-p)(1-\gamma)(1-p\eta)\gamma \\ &=\gamma + p(1-\gamma)\eta - p(1-p)(1-\gamma)\eta - (1-p)(1-\gamma)(1-p\eta)\gamma \\ &=\gamma + p(1-\gamma)\eta - (1-p)(1-\gamma)[\gamma + p(1-\gamma)\eta] \\ &=(p+\gamma-p\gamma)[\gamma + p(1-\gamma)\eta] \end{aligned} \tag{39}$$

Substituting (39) into (38) results in

$$\begin{split} \left. \frac{\mathrm{d}h\_{2}^{(m)}(x)}{\mathrm{d}x} \right|\_{x=1} &= M\_{2} \left( 1 - \frac{(1-p)(1-\gamma)(1-p\eta)\gamma}{(p+\gamma-p\gamma)[\gamma+p(1-\gamma)\eta]} + \frac{(1-p)(1-\gamma)}{p+\gamma-p\gamma} + \frac{(1-\gamma)(1-p\eta)}{\gamma+p(1-\gamma)\eta} \right) \\ &= M\_{2} \left( \frac{1}{p+\gamma-p\gamma} + \frac{p(1-\gamma)(1-p\eta)}{(p+\gamma-p\gamma)[\gamma+p(1-\gamma)\eta]} \right) \end{split} \tag{40}$$

which is exactly Equation (36).

Notice that in definition (20), for each *n* ≥ 1, the coefficient of *xn* is the probability that the AoI equals *n*. In order to obtain these coefficients, we decompose the PGF *HPP*(*x*) into power series. This shows that

$$\begin{split} H\_{PP}(\mathbf{x}) &= \frac{p\gamma M\_1 \mathbf{x}}{\gamma - p} \left( \frac{(1-p)\gamma}{1 - (1-p)\mathbf{x}} - \frac{p(1-\gamma)}{1 - (1-\gamma)\mathbf{x}} \right) \\ &- \left( \frac{(1-p)^2 \gamma^2 M\_2 \mathbf{x}^2}{(\gamma - p)[1 - (1-p)\mathbf{x}]} - \frac{p\gamma(1-p)(1-\gamma)M\_2 \mathbf{x}^2}{(\gamma - p)[1 - (1-\gamma)\mathbf{x}]} + \frac{p\gamma M\_2 \mathbf{x}^2}{[1 - (1-\gamma)\mathbf{x}]^2} \right) \\ &\times \left( \frac{\eta(1-p)(p+\gamma-p\gamma)}{(1-\eta)[1 - (1-p)(1-\gamma)\mathbf{x}]} - \frac{(1-p\eta)[\gamma + p(1-\gamma)\eta]}{(1-\eta)[1 - (1-\gamma)(1-p\eta)\mathbf{x}]} \right) \end{split} \tag{41}$$

when the preemption probability *η* = 1, while for the case *η* = 1, we have

$$\begin{split} H\_{PP}(\mathbf{x}) &= \frac{p\gamma M\_1 \mathbf{x}}{\gamma - p} \left( \frac{(1 - p)\gamma}{1 - (1 - p)\mathbf{x}} - \frac{p(1 - \gamma)}{1 - (1 - \gamma)\mathbf{x}} \right) \\ &+ \left( \frac{(1 - p)^2 \gamma^2 M\_2 \mathbf{x}^2}{(\gamma - p)[1 - (1 - p)\mathbf{x}]} - \frac{p\gamma (1 - p)(1 - \gamma)M\_2 \mathbf{x}^2}{(\gamma - p)[1 - (1 - \gamma)\mathbf{x}]} + \frac{p\gamma M\_2 \mathbf{x}^2}{[1 - (1 - \gamma)\mathbf{x}]^2} \right) \\ &\times \frac{\gamma + p^2 (1 - \gamma) - (1 - p)^2 (1 - \gamma)\gamma \mathbf{x}}{[1 - (1 - p)(1 - \gamma)\mathbf{x}]^2} \end{split} \tag{42}$$

The details of obtaining Equations (41) and (42) are given in Appendix B. In Section 4, along with the average value of the AoI, we determine the AoI's stationary distribution for two extreme cases: *η* = 0 and *η* = 1.

#### **4. Stationary Age of Information under Two Extreme Cases**

In this section, we determine the average AoI of the status updating system without packet preemption by setting *η* = 0, and when the preemption probability *η* is equal to 1, the mean of the AoI for the *Ber*/*Geo*/1/2∗ queue modeled system is also derived. In addition, using Equations (41) and (42), the stationary distributions of the discrete AoI for two cases are also obtained.

**Theorem 2.** *Assuming the packet arrivals form a Bernoulli process and the service time is geometrically distributed, the average AoIs of the discrete time status updating system with Ber*/*Geo*/1/2 *and Ber*/*Geo*/1/2∗ *queues are calculated as*

$$\overline{\Delta}\_{\text{Ber}/\text{Ga}/1/2} = \frac{1}{\gamma} \left( (1 - \gamma) + \frac{1}{\rho\_d} + \frac{2\rho\_d^2 (1 - \gamma)(1 - \gamma/2)}{1 + \rho\_d (1 - 2\gamma) + \rho\_d^2 (1 - \gamma)^2} \right) \tag{43}$$

*and*

$$\overline{\Delta}\_{\text{Ber}/\text{Geo}/1/2^\*} = \frac{1}{\gamma} \Biggl( (1-\gamma) + \frac{1}{\rho\_d} + \frac{\rho\_d^2 (1-\gamma) \left[ 1 + 3\rho\_d (1-\gamma) + \rho\_d^2 (1-\gamma)(1-2\gamma) \right]}{\left[ 1 + \rho\_d (1-2\gamma) + \rho\_d^2 (1-\gamma)^2 \right] \left[ 1 + \rho\_d (1-\gamma) \right]^2} \right) \tag{44}$$

*For each n* ≥ 1*, the distribution of the AoI* Δ*Ber*/*Geo*/1/2 *is given by*

$$\begin{split} \Pr\{\Delta\_{\text{Ber}} \zeta\_{\text{Gat}} | 1/2 = n\} &= \frac{p\gamma^2 (1-p) M\_1}{(\gamma - p)^2} (1-p)^n - \frac{(\gamma - p^2)(1-p)\gamma^2 M\_2}{(\gamma - p)^2} (1-\gamma)^{n-1} \\ &- \frac{p\gamma^2 (1-p) M\_2}{\gamma - p} (n-1)(1-\gamma)^{n-1} + \frac{p\gamma^2 M\_2}{2} n (n-1)(1-\gamma)^{n-2} \end{split} \tag{45}$$

*while when the system has full packet preemption, we show that*

$$\begin{aligned} &\Pr\{\Lambda\_{Rr/Gn/12^{n}} = n\} \\ &= \frac{p\gamma M\_{1}}{\gamma - p} \Big( \gamma (1-p)^{n} - p(1-\gamma)^{n} \Big) \\ &+ \frac{(1-p)[\gamma^{2} + p(1-\gamma)(p+\gamma)]M\_{2}}{\gamma - p} \Big( (1-p)^{n-1} - [(1-p)(1-\gamma)]^{n-1} \Big) \\ &- \frac{(1-p)\gamma[p+2(1-p)\gamma]M\_{2}}{\gamma - p} \Big( (1-\gamma)^{n-1} - [(1-p)(1-\gamma)]^{n-1} \Big) + p\gamma M\_{2} \Big( A(1-\gamma)^{n-2} \Big) \\ &+ B(n-1)(1-\gamma)^{n-2} + \mathbb{C}[(1-p)(1-\gamma)]^{n-2} + D(n-1)[(1-p)(1-\gamma)]^{n-2} \Big) \end{aligned} \tag{46}$$

*in which the coefficients A, B, C, and D are determined by*

$$A = \frac{2 - p}{p^3} \left( (1 - p)^2 \gamma - \frac{2(1 - p)[\gamma + p^2(1 - \gamma)]}{2 - p} \right) \tag{47}$$

$$C = -\frac{(1-p)(2-p)}{p^3} \left( (1-p)^2 \gamma - \frac{2(1-p)[\gamma + p^2(1-\gamma)]}{2-p} \right) \tag{48}$$

$$D = -\frac{p^2(1-p)\cdot A + (1-p)^2[\gamma + p^2(1-\gamma)]}{p(2-p)}\tag{49}$$

*and*

$$B = \left[\gamma + p^2(1-\gamma)\right] - A - \mathbb{C} - D \tag{50}$$

**Proof.** We first derive two average AoIs in Equations (43) and (44) from the general expression (37). Let *η* be 0; then, no packet preemption will occur in the system's buffer. The system's queue model reduces to *Ber*/*Geo*/1/2, and from (37), we can obtain the average AoI Δ*Ber*/*Geo*/1/2.

In this case, it is easy to show the last two terms within the brace of (37) can be calculated to be 1/*γ*. Thus, we have

$$\overline{\Delta}\_{\text{Ber}/\text{Go}/1/2} = \overline{\Delta}\_{\text{Ber}/\text{Go}/1/2^{\*} \, ^{\prime} \gamma \Big|\_{\eta=0}}$$

$$= \frac{1 + \rho\_d (1 - 2\gamma) + \rho\_d^2 (1 - \gamma)^2 + 2\rho\_d^3 (1 - \gamma)^2}{\rho\_d \gamma [1 + \rho\_d (1 - 2\gamma) + \rho\_d^2 (1 - \gamma)^2]}$$

$$\qquad + \frac{(1 - \gamma) + \rho\_d (1 - \gamma)(1 - 2\gamma) + \rho\_d^2 (1 - \gamma)(1 - \gamma + \gamma^2)}{\gamma [1 + \rho\_d (1 - 2\gamma) + \rho\_d^2 (1 - \gamma)^2]}$$

$$= \frac{1 + \rho\_d (2 - 3\gamma) + \rho\_d^2 (1 - \gamma)(2 - 3\gamma) + \rho\_d^3 (1 - \gamma)(3 - 3\gamma + \gamma^2)}{\rho\_d \gamma [1 + \rho\_d (1 - 2\gamma) + \rho\_d^2 (1 - \gamma)^2]}$$

$$= \frac{1}{\rho\_d \gamma} \Big( 1 + \rho\_d (1 - \gamma) + \frac{2\rho\_d^3 (1 - \gamma)(1 - \gamma/2)}{1 + \rho\_d (1 - 2\gamma) + \rho\_d^2 (1 - \gamma)^2} \Big) \tag{51}$$

$$=\frac{1}{\gamma}\left((1-\gamma)+\frac{1}{\rho\_d}+\frac{2\rho\_d^2(1-\gamma)(1-\gamma/2)}{1+\rho\_d(1-2\gamma)+\rho\_d^2(1-\gamma)^2}\right)\tag{52}$$

where in Equation (51) we use the method of long division.

For the other extreme case of *η* = 1, obviously the general expression (37) gives the average AoI Δ*Ber*/*Geo*/1/2∗ . Similarly, we first determine the value of the last two terms within the brace. We show that the difference of the last two terms equals

$$\frac{1}{\gamma} - \frac{\rho\_d^2 (1 - \gamma)}{\gamma [1 + \rho\_d (1 - \gamma)^2]} \tag{53}$$

thus, the average AoI Δ*Ber*/*Geo*/1/2∗ is calculated as

$$\begin{split} & \left. \overline{\Delta}\_{\text{Ber}/\text{Gav}/1/2^{\circ}} = \left. \overline{\Delta}\_{\text{Ber}/\text{Gav}/1/2^{\circ}} \right|\_{\mathfrak{T}=1} \\ &= \overline{\Delta}\_{\text{Ber}/\text{Gav}/1/2} - \frac{(1-\gamma) + \rho\_d(1-\gamma)(1-2\gamma) + \rho\_d^2(1-\gamma)(1-\gamma+\gamma^2)}{1+\rho\_d(1-2\gamma)+\rho\_d^2(1-\gamma)^2} \cdot \frac{\rho\_d^2(1-\gamma)}{\gamma[1+\rho\_d(1-\gamma)]^2} \\ &= \frac{1}{\gamma} \Big( (1-\gamma) + \frac{1}{1+\gamma} + \frac{\rho\_d^2(1-\gamma)(2-\gamma)}{1+\gamma[1+\gamma](1-\gamma)^2} \Big) \end{split} \tag{54}$$

$$\begin{split} \frac{1}{\gamma} \Big( (1-\gamma) + \frac{1}{\rho\_d} + \frac{\rho\_d^2 (1-\gamma)(2-\gamma)}{1 + \rho\_d (1-2\gamma) + \rho\_d^2 (1-\gamma)^2} \\ -\frac{\rho\_d^2 (1-\gamma)^2 + \rho\_d^3 (1-\gamma)^2 (1-2\gamma) + \rho\_d^4 (1-\gamma)^2 (1-\gamma+\gamma^2)}{\left[1 + \rho\_d (1-2\gamma) + \rho\_d^2 (1-\gamma)^2\right] \left[1 + \rho\_d (1-\gamma)\right]^2} \end{split} \tag{55}$$

$$=\frac{1}{\gamma}\left((1-\gamma)+\frac{1}{\rho\_d}+\frac{\rho\_d^2(1-\gamma)\left[1+3\rho\_d(1-\gamma)+\rho\_d^2(1-\gamma)(1-2\gamma)\right]}{1+\rho\_d(1-2\gamma)+\rho\_d^2(1-\gamma)^2}\right)\tag{56}$$

In Equation (54), since in the case of *η* = 0, the difference of the last two terms is 1/*γ*, the average AoI Δ*Ber*/*Geo*/1/2 is obtained. This equation also gives the exact gap between two average AoIs of the system with and without packet preemption. Notice that the latter term in (54) is always positive; then, the average AoI must become lower when the packet preemption strategy is applied. Equation (52) is substituted in (55), and in Equation (56), the expression of the average AoI Δ*Ber*/*Geo*/1/2∗ is finally determined.

Next, the distribution of the discrete AoI is calculated. Before the expressions (45) and (46) are derived, we first verify that both (45) and (46) are proper probability distributions by providing a specific numerical example.

Numerical results of two AoI distributions. Let *p* = 1/4 and *γ* = 1/2.

Firstly, from Equations (30) and (31), two numbers *M*1 and *M*2 are determined to be

$$M\_1 = \frac{12}{17}, \quad M\_2 = \frac{4}{17}.$$

after some simple calculations, for each *n* ≥ 1, the expression (45) gives

$$\Pr\{\Delta\_{\text{Rrr}/\text{Gro}/1/2} = n\} = \frac{9}{17} \left(\frac{3}{4}\right)^n - \frac{21}{68} \left(\frac{1}{2}\right)^n - \frac{3}{68} (n-1) \left(\frac{1}{2}\right)^{n-1} + \frac{1}{136} n (n-1) \left(\frac{1}{2}\right)^{n-2} \tag{57}$$

To obtain the numerical result of Equation (46), it is necessary to determine the four coefficients *A*, *B*, *C*, and *D* according to expressions (47)–(50). We directly find that

$$A = -\frac{39}{2}, \quad \text{C} = \frac{117}{8}, \quad D = \frac{45}{32}, \quad B = 4.$$

After some extra computations, it is shown that

$$\begin{split} \Pr\{\Delta\_{\text{Ber}/\text{Geo}/1/2^{\circ}} = n\} &= \frac{1}{2} \left(\frac{3}{4}\right)^{n} - \frac{105}{34} \left(\frac{1}{2}\right)^{n} \\ &+ \frac{57}{17} \left(\frac{3}{8}\right)^{n} + \frac{2}{17} (n-1) \left(\frac{1}{2}\right)^{n-2} + \frac{45}{1088} (n-1) \left(\frac{3}{8}\right)^{n-2} \end{split} \tag{58}$$

It can be checked directly that the sum of both (57) and (58) from *n* = 1 to ∞ are equal to 1. Therefore, expressions (45) and (46) indeed form the proper probability distributions.

In the following, by decomposing (41) and (42) further into several simplest rational fractions, we derive the explicit expressions of AoI distributions for the system with and without packet preemption.

First of all, for *η* = 0, it is easy to prove that the last part of (41) is equal to

$$\frac{-\gamma}{1 - (1 - \gamma)x}$$

thus, we have following equations. This shows that

*HPP*(*x*)|*η*=0 = *pγM*1*x γ* − *p* (1 − *p*)*γ* 1 − (1 − *p*)*x* − *p*(1 − *γ*) 1 − (1 − *γ*)*x* − (1 − *p*)2*γ*2*M*2*x*2 (*γ* − *p*)[1 − (1 − *p*)*x*] − *pγ*(1 − *p*)(1 − *γ*)*M*2*x*2 (*γ* − *p*)[1 − (1 − *γ*)*x*] + *pγM*2*x*2 [1 − (1 − *γ*)*x*] 2 · −*γ* 1 − (1 − *γ*)*x* = *pγ*2(1 − *p*)*M*1*x* (*γ* − *p*)[1 − (1 − *p*)*x*] − *p*2*γ*(1 − *γ*)*M*1*x* (*γ* − *p*)[1 − (1 − *γ*)*x*] + (1 − *p*)2*γ*3*M*2*x*2 (*γ* − *p*)[1 − (1 − *p*)*x*][1 − (1 − *γ*)*x*] − *pγ*2(1 − *p*)(1 − *γ*)*M*2*x*2 (*γ* − *p*)[1 − (1 − *γ*)*x*]2 + *pγ*2*M*2*x*2 [1 − (1 − *γ*)*x*]3 = *pγ*2(1 − *p*)*M*1*x* (*γ* − *p*)[1 − (1 − *p*)*x*] − *p*2*γ*(1 − *γ*)*M*1*x* (*γ* − *p*)[1 − (1 − *γ*)*x*] + (1 − *p*)2*γ*3*M*2*x*2 (*γ* − *p*) 1 − *p* (*γ* − *p*)[1 − (1 − *p*)*x*] − 1 − *γ* (*γ* − *p*)[1 − (1 − *γ*)*x*] − *pγ*2(1 − *p*)(1 − *γ*)*M*2*x*2 (*γ* − *p*)[1 − (1 − *γ*)*x*]2 + *pγ*2*M*2*x*2 [1 − (1 − *γ*)*x*]3 = *pγ*2(1 − *p*)*M*1*x γ* − *p* + (1 − *p*)3*γ*3*M*2*x*2 (*γ* − *p*)2 ∑∞ *n*=0[(1 − *p*)*x*] *n* − *p*2*γ*(1 − *γ*)*M*1*x γ* − *p* + (1 − *p*)2(1 − *γ*)*γ*3*M*2*x*2 (*γ* − *p*)2 ∑∞ *n*=0[(1 − *γ*)*x*] *n* − *pγ*2(1 − *p*)(1 − *γ*)*M*2*x*2 *γ* − *p* ∑∞ *n*=1 *n*[(1 − *γ*)*x*] *n*−1 + *pγ*2*M*2*x*2 2 ∑∞ *n*=2 *n*(*n* − 1)[(1 − *γ*)*x*] *n*−2 (59)

Taking the coefficient before *xn*, we find that

$$\begin{split} &\Pr\{\Delta\_{\text{RFr}}/\text{Ga}/1/2 = n\} \\ &= \left(\frac{p\gamma^{2}(1-p)M\_{1}}{\gamma-p}(1-p)^{n-1} + \frac{(1-p)^{3}\gamma^{3}M\_{2}}{(\gamma-p)^{2}}(1-p)^{n-2}\right) \\ & - \left(\frac{p^{2}\gamma(1-\gamma)M\_{1}}{\gamma-p}(1-\gamma)^{n-1} + \frac{(1-p)^{2}(1-\gamma)\gamma^{3}M\_{2}}{(\gamma-p)^{2}}(1-\gamma)^{n-2}\right) \\ & - \frac{p\gamma^{2}(1-p)(1-\gamma)M\_{2}}{\gamma-p}(n-1)(1-\gamma)^{n-2} + \frac{p\gamma^{2}M\_{2}}{2}n(n-1)(1-\gamma)^{n-2} \\ &= \frac{p\gamma^{2}(1-p)M\_{1}}{(\gamma-p)^{2}}(1-p)^{n} - \frac{(\gamma-p^{2})(1-p)\gamma^{2}M\_{2}}{(\gamma-p)^{2}}(1-\gamma)^{n-1} \\ & - \frac{p\gamma^{2}(1-p)M\_{2}}{\gamma-p}(n-1)(1-\gamma)^{n-1} + \frac{p\gamma^{2}M\_{2}}{2}n(n-1)(1-\gamma)^{n-2} \end{split} \tag{60}$$

This gives the stationary distribution (45) for the system without packet preemption.

On the other hand, when *η* is equal to 1, the system has full packet preemption. Factoring the PGF in Equation (42), we can also determine the stationary distribution of the AoI Δ*Ber*/*Geo*/1/2∗ by taking the coefficients of terms *xn* for each *n* ≥ 1. We give the explicit decomposition below, from which the distribution of AoI for the system with packet preemption is obtained explicitly.

From Equation (42), we show that

$$\begin{split} & \quad H\_{PP}(\mathbf{x})|\_{\eta=1} \\ &= \frac{p\gamma M\_{1}\mathbf{x}}{\gamma-p} \left( \frac{(1-p)\gamma}{1-(1-p)\mathbf{x}} - \frac{p(1-\gamma)}{1-(1-\gamma)\mathbf{x}} \right) + \frac{(1-p)^{2}\gamma^{2}M\_{2}\mathbf{x}^{2}}{\gamma-p} \left( \frac{\gamma^{2}+p(1-\gamma)(p+\gamma)}{\gamma^{2}[1-(1-p)\mathbf{x}]} \right) \\ & \quad - \frac{(1-\gamma)[\gamma^{2}+p(1-\gamma)(p+\gamma)]}{\gamma^{2}[1-(1-p)(1-\gamma)\mathbf{x}]} - \frac{p(1-\gamma)(p+\gamma-p\gamma)}{\gamma[1-(1-p)(1-\gamma)\mathbf{x}]^{2}} \\ & \quad - \frac{p\gamma(1-p)(1-\gamma)M\_{2}\mathbf{x}^{2}}{\gamma-p} \left( \frac{p+2(1-p)\gamma}{p[1-(1-\gamma)\mathbf{x}]} - \frac{(1-p)[p+2(1-p)\gamma]}{p[1-(1-p)(1-\gamma)\mathbf{x}]} \right) \\ & \quad - \frac{(1-p)(p+\gamma-p\gamma)}{[1-(1-p)(1-\gamma)\mathbf{x}]^{2}} + p\gamma M\_{2}\mathbf{x}^{2} \left( \frac{A}{1-(1-\gamma)\mathbf{x}} + \frac{B}{[1-(1-\gamma)\mathbf{x}]^{2}} \right) \\ & \quad + \frac{C}{1-(1-p)(1-\gamma)\mathbf{x}} + \frac{D}{[1-(1-p)(1-\gamma)\mathbf{x}]^{2}} \end{split} \tag{61}$$

in which we determine

$$A = \frac{2 - p}{p^3} \left( (1 - p)^2 \gamma - \frac{2(1 - p)[\gamma + p^2(1 - \gamma)]}{2 - p} \right) \tag{62}$$

$$C = -\frac{(1-p)(2-p)}{p^3} \left( (1-p)^2 \gamma - \frac{2(1-p)[\gamma + p^2(1-\gamma)]}{2-p} \right) \tag{63}$$

$$D = -\frac{p^2(1-p)\cdot A + (1-p)^2[\gamma + p^2(1-\gamma)]}{p(2-p)}\tag{64}$$

and

$$B = \left[\gamma + p^2(1 - \gamma)\right] - A - \mathbb{C} - D \tag{65}$$

Obtaining the second and the third row of (61) is not hard, while for the last row, we give some derivation details in Appendix C. Following the same procedures as those used to obtain (60), according to Equation (61), the probability that a stationary AoI equals each *n* is determined by the coefficient of the term *xn*.

*HPP*(*x*)|*η*=1 = *pγM*1 *γ* − *p* " (1 − *p*)*γ*(1 − *p*)*n*−1 − *p*(1 − *γ*)(1 − *γ*)*n*−1 # + (1 − *p*)2*γ*2*M*2 *γ* − *p* × " *γ*2 + *p*(1 − *γ*)(*p* + *γ*) *γ*2 (1 − *p*)*n*−2 − (1 − *γ*)[*γ*2 + *p*(1 − *γ*)(*p* + *γ*)] *γ*2 [(1 − *p*)(1 − *γ*)]*n*−2 − *p*(1 − *γ*)(*p* + *γ* − *pγ*) *γ* (*n* − 1)[(1 − *p*)(1 − *γ*)]*n*−2 # − *pγ*(1 − *p*)(1 − *γ*)*M*2 *γ* − *p* × " *p* + 2(1 − *p*)*γ p* (1 − *γ*)*n*−2 − (1 − *p*)[*p* + 2(1 − *p*)*γ*] *p* [(1 − *p*)(1 − *γ*)]*n*−2 − (1 − *p*)(*p* + *γ* − *pγ*)(*n* − 1)[(1 − *p*)(1 − *γ*)]*n*−2 # + *pγM*2 " *A*(1 − *γ*)*n*−2 + *B*(*n* − 1)(1 − *γ*)*n*−2 + *C*[(1 − *p*)(1 − *γ*)]*n*−2 + *D*(*n* − 1)[(1 − *p*)(1 − *γ*)]*n*−2 # = *pγM*1 *γ* − *p* " *γ*(1 − *p*)*n* − *p*(1 − *γ*)*n* # + (1 − *p*)[*γ*2 + *p*(1 − *γ*)(*p* + *γ*)]*M*2 *γ* − *p* " (1 − *p*)*n*−1 − [(1 − *p*)(1 − *γ*)]*n*−1 # − (1 − *p*)*γ*[*p* + 2(1 − *p*)*γ*]*M*2 *γ* − *p* " (1 − *γ*)*n*−1 − [(1 − *p*)(1 − *γ*)]*n*−1 # + *pγM*2 " *A*(1 − *γ*)*n*−2 + *B*(*n* − 1)(1 − *γ*)*n*−2 + *C*[(1 − *p*)(1 − *γ*)]*n*−2 + *D*(*n* − 1)[(1 − *p*)(1 − *γ*)]*n*−2 # (66)

which determines the distribution of AoI Δ*Ber*/*Geo*/1/2∗ .

So far, in Equations (52), (56), (60) and (66), we have obtained all the results in Theorem 2; thus, the proof is completed.

Actually, we have obtained the explicit expression of AoI's distribution for the system with packet preemption in our early work [60]. Earlier in this paper, we explain that solving the stationary equations is feasible for the easy situations but cannot be generalized when the system structure or queue models become complex. In [60], we focused on the discrete time system with three queues, i.e., the *Ber*/*Geo*/1/1, *Ber*/*Geo*/1/2, and *Ber*/*Geo*/1/2∗, and named them "*discrete packet management strategies*". There, we determined the AoI's stationary distribution for each system, and all the cases are dealt with by solving the stationary equations directly. Although the calculations are long—even tedious—these methods still have great significance, especially when the general status updating system is considered where the packet arrival process or the packet service process is arbitrary. It is with these methods that the analysis of discrete AoI can break through the limitation of the memoryless property that is imposed on the packet arrival and packet service processes in the SHS approach.

In [9], based on graphical arguments of the age process, the authors determined the average continuous AoI for the system with *M*/*M*/1/2 and *M*/*M*/1/2∗ queues as

$$
\Delta\_{M/M/1/2} = \frac{1}{\mu} \left( 1 + \frac{1}{\rho} + \frac{2\rho^2}{1 + \rho + \rho^2} \right) \tag{67}
$$

and

$$
\Delta\_{M/M/1/2^\*} = \frac{1}{\mu} \left( 1 + \frac{1}{\rho} + \frac{\rho^2 (1 + 3\rho + \rho^2)}{(1 + \rho + \rho^2)(1 + \rho)^2} \right) \tag{68}
$$

In addition, in previous work [58], we have proved that the mean of the AoI for a bufferless discrete time status updating system is equal to

$$\overline{\Delta}\_{\text{Ber}/\text{Geo}/1/1} = \frac{1}{\gamma} \left( (1 - \gamma) + \frac{1}{\rho\_d} + \frac{\rho\_d}{1/(1 - \gamma) + \rho\_d} \right) \tag{69}$$

while the corresponding continuous system with an *M*/*M*/1/1 queue has the average AoI

$$\overline{\Delta}\_{M/M/1/1} = \frac{1}{\mu} \left( 1 + \frac{1}{\rho} + \frac{\rho}{1+\rho} \right) \tag{70}$$

which was also given in [9].

We list Equations (43), (44) and (67)–(70) in Table 3—notice that this table has been given previously in Table 1 except for the last row, which gives the average continuous and discrete AoI for an infinite size status updating system. The mean of the discrete AoI Δ*Ber*/*Geo*/1/∞ was obtained recently in our work [59]. It is observed that apart from some additional product factors, the expressions of discrete AoI means for the system with Bernoulli packet arrivals and geometric service times are identical to those of the continuous system's average AoI, which uses the Poisson-exponential assumptions. So far, we have obtained enough evidence to propose the following relationship between the mean of discrete and continuous AoIs:

$$
\mu \cdot \overline{\Lambda}\_{M/M/1/\mathfrak{c}} = \gamma \cdot \overline{\Lambda}\_{\text{Ber}/\text{Ge}\alpha/1/\mathfrak{c}}|\_{\gamma=0} \text{ then replacing } \rho\_d \text{ by } \rho \tag{71}
$$

It is interesting and meaningful to verify the correspondence (71) by calculating the average AoI for more continuous time and discrete time systems. For example, determining the mean of AoI assuming the *M*/*M*/1/*c* queue is used in the continuous time system and for the discrete time systems with general *Ber*/*Geo*/1/*c* queues, where the system's size *c* is larger than 2.

**Table 3.** Some formulas of the average continuous and average discrete age of information.

#### **5. Numerical Simulation**

We provide the numerical results in this section. For general preemption probability, in the first two plots of Figure 3, we illustrate the relationships between the average AoI Δ*Ber*/*Geo*/1/2∗/*η* and the packet preemption probability *η*, and the traffic load *ρd*, respectively. The means of three discrete AoIs including Δ*Ber*/*Geo*/1/1, Δ*Ber*/*Geo*/1/2, and Δ*Ber*/*Geo*/1/2∗ are plotted in Figure 3c. For comparison, we also provide the numerical simulations of corresponding average continuous AoIs. At last, for three discrete AoIs, we depict their distribution curves and the cumulative probabilities in Figure 4. Notice that in our work [58], the distribution of the AoI for the bufferless system was obtained as

$$\Pr\{\Delta\_{\text{Ber}} / \text{Gco} / 1/1 = n\} = \frac{p(1-p)\gamma^3[(1-p)^n - (1-\gamma)^n]}{(p+\gamma-p\gamma)(\gamma-p)^2} - \frac{(p\gamma)^2 n(1-\gamma)^n}{(p+\gamma-p\gamma)(\gamma-p)}\tag{72}$$

For three different traffic loads *ρd*, we first draw the graphs between the average AoI Δ*Ber*/*Geo*/1/2∗/*η* and the preemption probability *η*. It is understandable that replacing the packet in a buffer with a fresher one can decrease the average AoI at the destination, and the numerical results in Figure 3a show that this trend is consistent as the preemption probability becomes large; that is, the mean of the AoI is decreasing monotonically when *η* increases. We mark the values at two extreme points where *η* = 0 and *η* = 1, which gives the average AoI Δ*Ber*/*Geo*/1/2 and Δ*Ber*/*Geo*/1/2∗ . Notice that the closer to *η* = 0, the more similar the behavior of the system becomes to that of a system using *Ber*/*Geo*/1/2 queues, and when *η* gradually gets to 1, a status updating system with *Ber*/*Geo*/1/2∗ queue is finally obtained. The three curves in Figure 3a also show that as the traffic load *ρd* increases from 0.4 to 0.45, the average AoI of the system with probabilistic packet preemption is reduced; thus, the timeliness performance is improved.

In Figure 3b, for three settings of preemption probabilities, i.e., *η* = 0, *η* = 0.5 and *η* = 1, the relationships of the average AoI versus traffic intensity *ρd* are illustrated. The topmost curve gives the average AoI of the system without packet preemption because for the case of *η* = 0, the average AoI reduces to Δ*Ber*/*Geo*/1/2. On the other hand, the curve at the bottom corresponds to the AoI's mean of the system that has full packet preemption. In order to make the differences among these graphs more significant, we draw the results in the range *ρd* ≥ 0.45. Three curves in Figure 3b clearly show that the timeliness of the system with complete packet preemption is the best, since when *η* is set to 1, the system's average AoI is the lowest. Since the results in Figure 3a show that the average AoI is monotonically decreasing when *η* increases, the graphs of the AoI's mean for a system with probabilistic packet preemption is located between the blue and the black lines in Figure 3b, such as the red line, which denotes the average AoI Δ*Ber*/*Geo*/1/2∗/0.5. In addition, the gaps between these curves are not significant for small *ρd*s but become large as *ρd* increases.

**Figure 3.** (**a**) Average AoI versus preemption probability *η* (different traffic load). (**b**) Average AoI versus traffic load *ρd* (different preemption intensity). (**c**) Comparisons of average discrete AoI and average continuous AoI.

**Figure 4.** (**a**) Stationary distributions of discrete AoI for bufferless system and the system with and without packet preemption. (**b**) The cumulative probabilities of three discrete AoIs.

From *ρd* = 0.15 to 0.9, we depict both the average discrete AoIs and the corresponding continuous average AoIs in Figure 3c for a bufferless system and size 2 status updating system with and without preemption. Continuous AoIs are denoted by dashed lines, and we use solid lines to represent the discrete AoIs. First of all, all the curves are decreasing as *ρd* becomes large, and the gaps between them are gradually apparent. For three continuous AoIs, it is observed that the average AoI Δ*M*/*M*/1/2∗ is the lowest in all the range of traffic load *ρ*. For the other two status updating systems, it is found that the system with an *M*/*M*/1/2 queue has a lower average AoI when *ρ* takes small values, while for high *ρ*, the average AoI of the bufferless system is smaller, and thus the timeliness is better. These results are the same for the graphs of discrete AoIs. Notice that when the discrete traffic intensity *ρd* is extremely large (near 0.9), the numerical results show that the average AoI Δ*Ber*/*Geo*/1/1 can be even smaller than Δ*Ber*/*Geo*/1/2∗ .

In Figure 3c, both continuous and discrete average AoIs are monotonically decreasing in the whole range of *ρd*; however, the monotonicity of the curve between the average AoI and *ρd* can only be maintained for small-size status updating systems. It is known that the average AoI of an infinite size system, i.e., Δ*M*/*M*/1/∞, is not monotonic when the traffic load varies from 0 to 1. Thus, for a size *c* status updating system with Bernoulli packet arrivals and geometrically distributed service time, there must be a *critical size c*∗ such that when *c* < *c*∗, the mean of the system's AoI Δ*Ber*/*Geo*/1/*c* is monotonically decreasing as *ρd* tends to 1. In contrast, for those cases where the system size *c* ≥ *c*∗, the curve has a valley, and an optimal *ρd* exists at which the average AoI is minimized. Similarly, for the continuous average AoI Δ*M*/*M*/1/*c* of the system with general size *c*, a *c*∗ also exists so that Δ*M*/*M*/1/*c* is always decreasing when *c* < *c*∗ and the graph of Δ*M*/*M*/1/*c* first falls and then rises for those *c*s where *c* ≥ *c*∗. In addition, from the alternation of Δ*M*/*M*/1/1 and Δ*M*/*M*/1/2, and also of Δ*Ber*/*Geo*/1/1 and Δ*Ber*/*Geo*/1/2, we can infer that, as a function of *c*, the graphs of Δ*M*/*M*/1/*c* and Δ*Ber*/*Geo*/1/*c* are not monotonic.

At last, the distribution curves and cumulative probabilities of three discrete AoIs are depicted in Figure 4, in which we set a relatively large *ρd* to make the difference between them clear. On the whole, these curves are similar. In Figure 4a, from the distributions of AoI Δ*Ber*/*Geo*/1/1 to that of Δ*Ber*/*Geo*/1/2, the peak of the curve decreases and the point at which the peak stationary probability is achieved moves slightly to the right. As the AoI becomes large, the distribution curve of the system with the *Ber*/*Geo*/1/1 queue drops more sharply. The distribution corresponding to Δ*Ber*/*Geo*/1/2∗ has the largest peak value pf all of three discrete AoIs, and the descent speed is the fastest when the value of AoI is large. In addition, it seems that this maximal probability is taken at the same discrete AoI as that of the distribution of Δ*Ber*/*Geo*/1/2. We also provide the cumulative probabilities of three discrete AoIs in Figure 4b.

#### **6. Conclusions**

In this paper, we consider the stationary AoI of a size 2 status updating system where the packet waiting in the buffer can be preempted by fresher packets with the given probability *η*. We show that this phenomenon may occur in the energy-harvest (EH) nodes of wireless sensor networks where the charging process is stochastic. We constitute a three-dimensional age process and derive the general expression of the system's average AoI using the PGF method. Let *η* = 0 and *η* = 1; the mean of two discrete AoIs Δ*Ber*/*Geo*/1/2 and Δ*Ber*/*Geo*/1/2 are determined, and the exact distribution expressions of both AoIs are also obtained by writing the PGF as the power series.

We propose the idea and methods for the analysis of discrete AoIs—that is, constituting multiple-dimensional age processes and applying the PGF method. A detailed introduction is given to exhibit the usage of the idea and methods to more discrete time status updating systems. With this paper, we have shown how the AoI of basic discrete system is characterized, while in further work, we will focus on the age analysis of systems with a more general structure, such as systems with multi-sources and systems with multihop packet transmission. As one part of the AoI theory, we believe that the research into discrete AoI deserves more attention.

**Author Contributions:** Writing—original draft preparation, J.Z.; writing—review and editing, Y.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded in part by the National Natural Science Foundation of China under Grant 61901105 and Grant 62171119, in part by the Natural Science Foundation of Jiangsu Province under Grant BK20190343, in part by the Zhi Shan Young Scholar Program of Southeast University.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Proof of Theorem 1**

In this appendix, we derive all the results in Lemma 1 from the stationary Equations (18).

Define that

$$M\_1 = \sum\_{n=1}^{\infty} \pi\_{(n,0,0)}, \; M\_2 = \sum\_{m=1}^{\infty} \sum\_{n=m+1}^{\infty} \pi\_{(n,m,0)}, \; \text{and} \\ M\_3 = \sum\_{l=1}^{\infty} \sum\_{m=l+1}^{\infty} \sum\_{n=m+1}^{\infty} \pi\_{(n,m,l)}, \; \text{and} \\ M\_4 = \sum\_{l=1}^{\infty} \sum\_{n=l+1}^{\infty} \pi\_{(n,l,n)}$$

These three numbers are determined at the first place. According to the last two rows of (18), we have

$$\begin{aligned} M\_1 &= \pi\_{(1,0,0)} + \sum\_{n=2}^{\infty} \pi\_{(n,0,0)} \\ &= \left(\sum\_{n=1}^{\infty} \pi\_{(n,0,0)}\right) p\gamma + \sum\_{n=2}^{\infty} \left\{\pi\_{(n-1,0,0)}(1-p) + \sum\_{k=n}^{\infty} \pi\_{(k,n-1,0)}(1-p)\gamma\right\} \end{aligned}$$

$$\hat{\rho} = p\gamma M\_1 + (1-p)M\_1 + \sum\_{\mathfrak{m}=1}^{\infty} \sum\_{\mathfrak{n}=\mathfrak{m}+1}^{\infty} \pi\_{(\mathfrak{d},\mathfrak{M},0)}(1-p)\gamma \tag{A1}$$

$$0 = p\gamma M\_1 + (1-p)M\_1 + (1-p)\gamma M\_2 \tag{A2}$$

where in (A1), we have used the substitutions *k* = *n*˜ and *n* − 1 = *m*˜ . From Equation (A2), we obtain the first relation

$$p(1 - \gamma)M\_1 = (1 - p)\gamma M\_2 \tag{A3}$$

Next, we deal with the number *M*3 as follows.

$$\begin{split} M\_3 &= \sum\_{l=1}^{\infty} \sum\_{m=l+1}^{\infty} \sum\_{n=m+1}^{\infty} \pi\_{(n,m,l)} \\ &= \sum\_{m=2}^{\infty} \sum\_{n=m+1}^{\infty} \pi\_{(n,m,1)} + \sum\_{l=2}^{\infty} \sum\_{m=l+1}^{\infty} \sum\_{n=m+1}^{\infty} \pi\_{(n,m,l)} \end{split} \tag{A4}$$

Using the first row of (18), the latter sum in Equation (A4) equals

$$\begin{split} & \sum\_{l=2}^{\infty} \sum\_{m=l+1}^{\infty} \sum\_{n=m+1}^{\infty} \pi\_{(n-1,m-1,l-1)} (1-\gamma)(1-p\eta) \\ &= \sum\_{l=1}^{\infty} \sum\_{m=l+1}^{\infty} \sum\_{n=m+1}^{\infty} \pi\_{(n,m,l)} (1-\gamma)(1-p\eta) \\ &= (1-\gamma)(1-p\eta)M\_{3} \end{split} \tag{A5}$$

and from the second and the third row of (18), the first part of (A4) is calculated as

$$\begin{split} &\sum\_{n=3}^{\infty} \pi\_{(n,2,1)} + \sum\_{m=3}^{\infty} \sum\_{n=m+1}^{\infty} \pi\_{(n,m,1)} \\ &= \sum\_{n=3}^{\infty} \pi\_{(n-1,0,0)} p(1-\gamma) \\ &\quad + \sum\_{m=3}^{\infty} \sum\_{n=m+1}^{\infty} \left\{ \pi\_{(n-1,m-1,0)} p(1-\gamma) + \sum\_{j=1}^{m-2} \pi\_{(n-1,m-1,j)} p(1-\gamma) \eta \right\} \\ &= \sum\_{n=2}^{\infty} \pi\_{(n,1,0)} p(1-\gamma) \\ &\quad + \sum\_{m=2}^{\infty} \sum\_{n=m+1}^{\infty} \pi\_{(n,w,0)} p(1-\gamma) + \sum\_{m=3}^{\infty} \sum\_{n=m+1}^{\infty} \sum\_{j=1}^{m-2} \pi\_{(n-1,m-1,j)} p(1-\gamma) \eta \\ &= \sum\_{m=1}^{\infty} \sum\_{n=m+1}^{\infty} \pi\_{(n,w,0)} p(1-\gamma) + \sum\_{m=2}^{\infty} \sum\_{n=0, n\neq 1}^{\infty} \sum\_{l=1}^{n-1} \pi\_{(0,0,l)} p(1-\gamma) \eta \end{split} \tag{A6}$$

where in (A6), we let *n* − 1 = *n*˜, *m* − 1 = *m*˜ , and *j* = ˜ *l*. Equations (A4), (A5) and (A7) together give

$$p(1 - \gamma)M\_2 = \gamma M\_3 \tag{A8}$$

Since the sum of all the stationary probabilities equals 1, thus we have

$$M\_1 + M\_2 + M\_3 = 1\tag{A9}$$

Combining Equations (A3), (A8) and (A9), we can solve that

$$M\_1 = \frac{(1-p)\gamma^2}{(p+\gamma-2p\gamma)\gamma+p^2(1-\gamma)^2} \tag{A10}$$

$$M\_2 = \frac{p\gamma(1-\gamma)}{(p+\gamma-2p\gamma)\gamma+p^2(1-\gamma)^2} \tag{A11}$$

$$M\_3 = \frac{p^2(1-\gamma)^2}{(p+\gamma-2p\gamma)\gamma+p^2(1-\gamma)^2} \tag{A12}$$

We mention that from the fourth, the fifth, and the sixth equation in (18), another relation can also be obtained for the second number *M*2, which is given directly as

$$M\_2 = p(1 - \gamma)M\_1 + p\gamma M\_2 + p\gamma \eta M\_3 + (1 - p)(1 - \gamma)M\_2 + (1 - p\eta)\gamma M\_3 \tag{A13}$$

which is reduced to (A8) when eliminating *M*1 using the relation (A3).

Then, the relationships between functions *hi*(*x*), 1 ≤ *i* ≤ 3 and *h* (*m*) 2 (*x*) are determined through similar procedures. First of all, we see that

$$\begin{split} h\_{1}(\mathbf{x}) &= \sum\_{n=1}^{\infty} \mathbf{x}^{\mu} \pi\_{(\mathbf{u},0,0)} \\ &= \mathbf{x} \pi\_{(1,0,0)} + \sum\_{n=2}^{\infty} \mathbf{x}^{\mu} \left\{ \pi\_{(n-1,0,0)} (1-p) + \sum\_{k=n}^{\infty} \pi\_{(k,n-1,0)} (1-p) \gamma \right\} \\ &= p \gamma M\_{1} \mathbf{x} + \sum\_{n=1}^{\infty} \mathbf{x}^{\mu+1} \pi\_{(n,0,0)} (1-p) + \sum\_{n=2}^{\infty} \mathbf{x}^{\mu} \sum\_{k=n}^{\infty} \pi\_{(k,n-1,0)} (1-p) \gamma \\ &= p \gamma M\_{1} \mathbf{x} + (1-p) x h\_{1}(\mathbf{x}) + \sum\_{n=1}^{\infty} \mathbf{x}^{\mu+1} \sum\_{k=n+1}^{\infty} \pi\_{(k,n,0)} (1-p) \gamma \end{split} \tag{A14}$$

$$=p\gamma\mathcal{M}\_1\mathbf{x} + (1-p)\mathbf{x}h\_1(\mathbf{x}) + (1-p)\gamma\mathbf{x}h\_2^{(m)}(\mathbf{x})\tag{A15}$$

in which we denote *k* = *n*˜ and *n* − 1 = *m*˜ .

From (A15), we obtain

$$h\_1(\mathbf{x}) = \frac{p\gamma M\_1 \mathbf{x}}{1 - (1 - p)\mathbf{x}} + \frac{(1 - p)\gamma \mathbf{x}}{1 - (1 - p)\mathbf{x}} h\_2^{(m)}(\mathbf{x}) \tag{A16}$$

Using stationary Equation (18), we determine function *h*2(*x*) in the following.

*h*2(*x*) = ∑∞ *m*=1 ∑∞ *n*=*m*+1 *xnπ*(*n*,*m*,0) = ∑∞ *n*=2 *xnπ*(*n*,1,0) + ∑∞ *m*=2 ∑∞ *n*=*m*+1 *xn* × ) *π*(*n*−1,*m*−1,0)(1 − *p*)(1 − *γ*) + ∑∞ *k*=*n π*(*k*,*n*−1,*m*−1)(1 − *pη*)*γ* \* = *x*2*π*(2,1,0) + ∑∞ *n*=3 *xn* ) *π*(*n*−1,0,0) *p*(1 − *γ*) + ∑∞ *k*=*n π*(*k*,*n*−1,0) *pγ* + ∑∞ *k*=*n* ∑*n*−2 *j*=1 *π*(*k*,*n*−1,*j*) *pγη*\* + ∑∞ *m*=1 ∑∞ *n*=*m*+1 *xn*+1*π*(*n*,*m*,0)(1 − *p*)(1 − *γ*) + ∑∞ *m*=2 ∑∞ *n*=*m*+1 *xn* ∑∞ *k*=*n π*(*k*,*n*−1,*m*−1)(1 − *pη*)*γ* = *x*2 ) *π*(1,0,0) *p*(1 − *γ*) + ∑∞ *k*=2 *π*(*k*,1,0) *pγ* \* + ∑∞ *n*=2 *xn*+1*π*(*n*,0,0) *p*(1 − *γ*) + ∑∞ *n*=3 *xn* ∑∞ *k*=*n π*(*k*,*n*−1,0) *pγ* + ∑∞ *n*=3 *xn* ∑∞ *k*=*n* ∑*n*−2 *j*=1 *π*(*k*,*n*−1,*j*) *pγη* + (1 − *p*)(1 − *γ*)*xh*2(*x*)+(1 − *pη*)*γxh*(*m*) 3 (*x*) (A17)

Let *k* = *n*˜, *n* − 1 = *m*˜ , and *m* − 1 = ˜ *l*, we have

$$\begin{aligned} &\sum\_{m=2}^{\infty} \sum\_{m=m+1}^{\infty} \mathbf{x}^{n} \sum\_{k=n}^{\infty} \pi\_{(k,n-1,m-1)} (1 - p\eta) \gamma \\ &= \sum\_{l=1}^{\infty} \sum\_{m=l+1}^{\infty} \mathbf{x}^{m+1} \sum\_{\tilde{n}=\tilde{m}+1}^{\infty} \pi\_{(0,n,l)} (1 - p\eta) \gamma \\ &= (1 - p\eta) \gamma x h\_{3}^{(m)}(\mathbf{x}) \end{aligned} \tag{A18}$$

where we define

$$h\_3^{(m)}(\mathbf{x}) = \sum\_{l=1}^{\infty} \sum\_{m=l+1}^{\infty} \sum\_{n=m+1}^{\infty} \mathbf{x}^{\mathbf{m}} \pi\_{(n,m,l)} \tag{A19}$$

and the last term in Equation (A17) is obtained.

Continuing the calculation of (A17), we obtain that

$$\begin{split} h\_{2}(\mathbf{x}) &= p(1-\gamma)\mathbf{x}h\_{1}(\mathbf{x}) + p\gamma \mathbf{x}h\_{2}^{(m)}(\mathbf{x}) + \sum\_{\theta \neq -2}^{\infty} \mathbf{x}^{\theta \mathbf{i} + 1} \sum\_{\theta = \theta + 1}^{\infty} \sum\_{l=1}^{\tilde{m}-1} \pi\_{(\theta, \theta, l)} p\gamma \eta \\ &+ (1-p)(1-\gamma)\mathbf{x}h\_{2}(\mathbf{x}) + (1-p\eta)\gamma \mathbf{x}h\_{3}^{(m)}(\mathbf{x}) \end{split} \tag{A20}$$

In (A20), we have used the substitutions *k* = *n*˜, *n* − 1 = *m*˜ , and *j* = ˜ *l*. Except for the factor *pγηx*, the sum in (A20) is equal to

$$\begin{split} & \quad \mathbf{x}^{2} \sum\_{n=3}^{\infty} \pi\_{(n,2,1)} + \mathbf{x}^{3} \sum\_{n=4}^{\infty} \left[ \pi\_{(n,3,1)} + \pi\_{(n,3,2)} \right] + \mathbf{x}^{4} \sum\_{n=5}^{\infty} \left[ \pi\_{(n,4,1)} + \pi\_{(n,4,2)} + \pi\_{(n,4,3)} \right] + \cdots \\ & = \sum\_{m=2}^{\infty} \mathbf{x}^{m} \sum\_{n=m+1}^{\infty} \pi\_{(n,m,1)} + \sum\_{m=3}^{\infty} \mathbf{x}^{m} \sum\_{n=m+1}^{\infty} \pi\_{(n,m,2)} + \sum\_{m=4}^{\infty} \mathbf{x}^{m} \sum\_{n=m+1}^{\infty} \pi\_{(n,m,3)} + \cdots \\ & = \sum\_{l=1}^{\infty} \sum\_{m=l+1}^{\infty} \pi\_{m=l+1} \mathbf{x}^{m} \sum\_{n=m+1}^{\infty} \pi\_{(n,m,l)} \\ & = h\_{3}^{(m)}(\mathbf{x}) \end{split} \tag{A21}$$

Substituting the result (A21) into Equation (A20) and merging the same terms gives

$$h\_2(\mathbf{x})[1 - (1 - p)(1 - \gamma)\mathbf{x}] = p(1 - \gamma)\mathbf{x}h\_1(\mathbf{x}) + p\gamma \mathbf{x}h\_2^{(\mathbf{m})}(\mathbf{x}) + \gamma \mathbf{x}h\_3^{(\mathbf{m})}(\mathbf{x})\tag{A.22}$$

We compute *h* (*m*) 3 (*x*) while determining the other function *h* (*m*) 2 (*x*) in the end. As before, from the equations in (18), we find that

*h* (*m*) 3 (*x*) = ∑∞ *l*=1 ∑∞ *m*=*l*+1 ∑∞ *n*=*m*+1 *xmπ*(*n*,*m*,*l*) = ∑∞ *m*=2 ∑∞ *n*=*m*+1 *xmπ*(*n*,*m*,1) + ∑∞ *l*=2 ∑∞ *m*=*l*+1 ∑∞ *n*=*m*+1 *xmπ*(*n*−1,*m*−1,*l*−1)(1 − *γ*)(1 − *pη*) = ∑∞ *n*=3 *x*2*π*(*n*,2,1) + ∑∞ *m*=3 ∑∞ *n*=*m*+1 *xm* ) *π*(*n*−1,*m*−1,0) *p*(1 − *γ*) + ∑*m*−2 *j*=1 *π*(*n*−1,*m*−1,*j*) *p*(1 − *γ*)*η* \* + (1 − *γ*)(1 − *pη*)*xh*(*m*) 3 (*x*) = ∑∞ *n*=3 *x*2*π*(*n*−1,1,0) *p*(1 − *γ*) + ∑∞ *m*=2 ∑∞ *n*=*m*+1 *xm*+1*π*(*n*,*m*,0) *p*(1 − *γ*) + ∑∞ *m*=3 ∑∞ *n*=*m*+1 *xm* ∑*m*−2 *j*=1 *π*(*n*−1,*m*−1,*j*) *p*(1 − *γ*)*η* + (1 − *γ*)(1 − *pη*)*xh*(*m*) 3 (*x*) = *p*(1 − *γ*)*xh*(*m*) 2 (*x*) + *p*(1 − *γ*)*ηxh*(*m*) 3 (*x*)+(1 − *γ*)(1 − *pη*)*xh*(*m*) 3 (*x*) = *p*(1 − *γ*)*xh*(*m*) 2 (*x*)+(1 − *γ*)*xh*(*m*) 3 (*x*) (A23)

which shows that

$$h\_3^{(m)}(\mathbf{x}) = \frac{p(1-\gamma)\mathbf{x}}{1-(1-\gamma)\mathbf{x}}h\_2^{(m)}(\mathbf{x})\tag{A24}$$

Combining Equations (A22) and (A24) yields the relation

$$h\_2(\mathbf{x}) = \frac{p(1-\gamma)\mathbf{x}}{1 - (1-p)(1-\gamma)\mathbf{x}}h\_1(\mathbf{x}) + \frac{p\gamma\mathbf{x}}{[1-(1-p)(1-\gamma)\mathbf{x}][1-(1-\gamma)\mathbf{x}]}h\_2^{(m)}(\mathbf{x}) \tag{A25}$$

Now, we deal with function *h*3(*x*). Similarly to the process of obtaining (A23), we have

$$\begin{split} h\_3(\mathbf{x}) &= \sum\_{m=-2}^{\infty} \sum\_{n=-m+1}^{n} \mathbf{x}^n \pi\_{(n,m,1)} + \sum\_{l=-2}^{\infty} \sum\_{m=-l+1}^{\infty} \sum\_{n=-m+1}^{\infty} \mathbf{x}^l \pi\_{(n-1,m-1,l-1)} (1-\gamma)(1-p\eta) \\ &= \sum\_{m=-3}^{\infty} \mathbf{x}^n \pi\_{(n,2,1)} + \sum\_{m=-3}^{\infty} \sum\_{n=-m+1}^{\infty} \mathbf{x}^n \left\{ \pi\_{(n-1,m-1,0)} p(1-\gamma) \\ &\quad + \sum\_{j=1}^{m-2} \pi\_{(n-1,m-1,j)} p(1-\gamma) \eta \right\} + (1-\gamma)(1-p\eta) \mathbf{x} h\_3(\mathbf{x}) \\ &= \sum\_{m=-3}^{\infty} \mathbf{x}^n \pi\_{(n-1,0,0)} p(1-\gamma) + \sum\_{m=-2}^{\infty} \sum\_{n=-m+1}^{\infty} \mathbf{x}^{n+1} \pi\_{(n,m,0)} p(1-\gamma) \\ &\quad + \sum\_{m=-3}^{\infty} \sum\_{n=-m+1}^{\infty} \mathbf{x}^n \sum\_{j=1}^{m-2} \pi\_{(n-1,m-1,j)} p(1-\gamma) \eta + (1-\gamma)(1-p\eta) \mathbf{x} h\_3^{\mathbf{x}}(\mathbf{x}) \\ &= p(1-\gamma) \mathbf{x} h\_2(\mathbf{x}) + p(1-\gamma) \eta \mathbf{x} h\_3^{\mathbf{x}}(\mathbf{x}) + (1-\gamma)(1-p\eta) \mathbf{x} h\_3^{\mathbf{x}}(\mathbf{x}) \\ &= p(1-\gamma) \mathbf{x} h\_2(\mathbf{x}) + (1-\gamma) \mathbf{x} h\_3^{\mathbf{$$

from which we derive

$$h\_3(\mathbf{x}) = \frac{p(1-\gamma)\mathbf{x}}{1-(1-\gamma)\mathbf{x}}h\_2(\mathbf{x})\tag{A27}$$

So far, we have obtained the relations (26)–(28) in Equations (A16), (A25) and (A27). To complete the proof of Lemma 1, the remaining part is the determination of the last function *h* (*m*) 2 (*x*). It is shown that

*h* (*m*) 2 (*x*) = ∑∞ *m*=1 ∑∞ *n*=*m*+1 *xmπ*(*n*,*m*,0) = ∑∞ *n*=2 *xπ*(*n*,1,0) + ∑∞ *m*=2 ∑∞ *n*=*m*+1 *xm* × ) *π*(*n*−1,*m*−1,0)(1 − *p*)(1 − *γ*) + ∑∞ *k*=*n π*(*k*,*n*−1,*m*−1)(1 − *pη*)*γ* \* = *xπ*(2,1,0) + ∑∞ *n*=3 *x* ) *π*(*n*−1,0,0) *p*(1 − *γ*) + ∑∞ *k*=*n π*(*k*,*n*−1,0) *pγ* + ∑∞ *k*=*n* ∑*n*−2 *j*=1 *π*(*k*,*n*−1,*j*) *pγη*\* + ∑∞ *m*=1 ∑∞ *n*=*m*+1 *xm*+1*π*(*n*,*m*,0)(1 − *p*)(1 − *γ*) + ∑∞ *m*=2 ∑∞ *n*=*m*+1 *xm* ∑∞ *k*=*n π*(*k*,*n*−1,*m*−1)(1 − *pη*)*γ* = *x* ) *π*(1,0,0) *p*(1 − *γ*) + ∑∞ *k*=2 *π*(*k*,1,0) *pγ* \* + ∑∞ *n*=2 *xπ*(*n*,0,0) *p*(1 − *γ*) + ∑∞ *n*=3 *x*∑∞ *k*=*n π*(*k*,*n*−1,0) *pγ* + ∑∞ *n*=3 *x*∑∞ *k*=*n* ∑*n*−2 *j*=1 *π*(*k*,*n*−1,*j*) *pγη* + (1 − *p*)(1 − *γ*)*xh*(*m*) 2 (*x*)+(1 − *pη*)*γxh*(*l*) 3 (*x*) (A28)

Similarly, we use the substitutions *k* = *n*˜, *n* − 1 = *m*˜ , and *m* − 1 = ˜ *l*, and we write

$$\begin{split} & \sum\_{m=2}^{\infty} \sum\_{n=m+1}^{\infty} \mathbf{x}^{m} \sum\_{k=n}^{\infty} \pi\_{(k,n-1,m-1)} (1 - p\eta) \gamma \\ &= \sum\_{l=1}^{\infty} \sum\_{m=l+1}^{\infty} \mathbf{x}^{l+1} \sum\_{n=m+1}^{\infty} \pi\_{(\tilde{n}, \tilde{m}, \tilde{l})} (1 - p\eta) \gamma \\ &= (1 - p\eta) \gamma \mathbf{x} h\_{3}^{(l)}(\mathbf{x}) \end{split} \tag{A29}$$

Thus, the last term in Equation (A28) is obtained, where *h* (*l*) 3 (*x*) is defined and determined as

$$\begin{split} h\_{3}^{(1)}(x) &= \sum\_{l=n+1}^{\infty} \sum\_{m=l+1}^{\infty} x^{l} \sum\_{n=m+1}^{\infty} \pi\_{(n,m)} \\ &= \sum\_{n=n+2}^{\infty} x \sum\_{m=n+1}^{\infty} \pi\_{(n,m,1)} + \sum\_{l=2}^{\infty} \sum\_{m=l+1}^{\infty} x^{l} \sum\_{n=m+1}^{\infty} \pi\_{(n-1,m-1,l-1)} (1-\gamma)(1-p\eta) \\ &= x \sum\_{n=0}^{\infty} \pi\_{(n,2)} + \sum\_{m=3}^{\infty} x \sum\_{n=m+1}^{\infty} \left\{ \pi\_{(n-1,m-1,0)} p(1-\gamma) \\ &\quad + \sum\_{j=1}^{m-2} \pi\_{(n-1,m-1,j)} p(1-\gamma) \eta \right\} + (1-\gamma)(1-p\eta) x h\_{3}^{(l)}(x) \\ &= x \sum\_{n=n}^{\infty} \pi\_{(n-1,0,0)} p(1-\gamma) + \sum\_{m=2}^{\infty} x \sum\_{n=m+1}^{\infty} \pi\_{(n,m,0)} p(1-\gamma) \\ &\quad + \sum\_{n=n+3}^{\infty} x \sum\_{m=n+1}^{\infty} \sum\_{j=1}^{m-2} \pi\_{(n-1,m-1,j)} p(1-\gamma) \eta + (1-\gamma)(1-p\eta) x h\_{3}^{(l)}(x) \\ &= p(1-\gamma)M\_{3}^{x} x + p(1-\gamma)\eta M\_{3}^{x} x + (1-\gamma)(1-p\eta) x h\_{3}^{(l)}(x) \\ &= [\gamma + p(1-\gamma)] \eta M\_{3}^{x} x + (1-\gamma)(1-p\eta) x h\_{3}^{(l)}(x) \end{split}$$

In Equation (A30), we have used the relation (A8), which says *p*(1 − *γ*)*M*2 = *γM*3. From (A30), we obtain the result

$$h\_3^{(l)}(\mathbf{x}) = \frac{[\gamma + p(1 - \gamma)\eta]M\_3\mathbf{x}}{1 - (1 - \gamma)(1 - p\eta)\mathbf{x}} \tag{A31}$$

Coming back to the calculation of function *h* (*m*) 2 (*x*), Equation (A28) shows

$$\begin{split} \frac{d}{2} \langle \mathbf{x} \rangle [1 - (1 - p)(1 - \gamma)\mathbf{x}] &= p(1 - \gamma)M\_1 \mathbf{x} + p\gamma M\_2 \mathbf{x} + p\gamma \eta M\_3 \mathbf{x} + (1 - p\eta)\gamma x h\_3^{(l)}(\mathbf{x}) \\ &= \left[ \gamma + p^2 (1 - \gamma)\eta \right] M\_2 \mathbf{x} + (1 - p\eta)\gamma x h\_3^{(l)}(\mathbf{x}), \end{split} \tag{A32}$$

in which we have used Equations (A3) and (A8) to replace numbers *M*1 and *M*3. Substituting expression (A31), after some extra operations, we determine that

$$h\_2^{(m)}(x) = \frac{\left[\gamma + p^2(1-\gamma)\eta - (1-p)(1-\gamma)(1-p\eta)\gamma x\right]M\_2x}{\left[1 - (1-p)(1-\gamma)x\right]\left[1 - (1-\gamma)(1-p\eta)x\right]}\tag{A33}$$

This eventually completes the proof of Lemma 1.

#### **Appendix B. Proof of Equations (41) and (42)**

In Equation (33), we show that

*h*

$$\begin{split} H\_{PP}(\mathbf{x}) &= \frac{p\gamma M\_1 \mathbf{x} [1 - (1 - p)(1 - \gamma)\mathbf{x}]}{[1 - (1 - p)\mathbf{x}][1 - (1 - \gamma)\mathbf{x}]} \\ &+ \frac{\gamma \mathbf{x} \left\{ 1 - (1 - p)[2(1 - \gamma) + p\gamma]\mathbf{x} + (1 - p)^2(1 - \gamma)^2 \mathbf{x}^2 \right\}}{[1 - (1 - p)\mathbf{x}][1 - (1 - \gamma)\mathbf{x}]^2} h\_2^{(m)}(\mathbf{x}) \end{split} \tag{A34}$$

Equations (41) and (42) are obtained by decomposing each part of (A34). For the first part, we assume

$$\begin{split} \frac{1 - (1 - p)(1 - \gamma)\mathbf{x}}{[1 - (1 - p)\mathbf{x}][1 - (1 - \gamma)\mathbf{x}]} &= \frac{A}{1 - (1 - p)\mathbf{x}} + \frac{B}{1 - (1 - \gamma)\mathbf{x}}\\ &= \frac{(A + B) - [A(1 - \gamma) + B(1 - p)]\mathbf{x}}{[1 - (1 - p)\mathbf{x}][1 - (1 - \gamma)\mathbf{x}]} \end{split} \tag{A35}$$

and according to the coefficients of corresponding terms, we obtain

$$A + B = 1, \quad A(1 - \gamma) + B(1 - p) = (1 - p)(1 - \gamma) \tag{A36}$$

which determine *A* and *B* as

$$A = \frac{(1-p)\gamma}{\gamma - p}, \text{ and } B = \frac{p(1-\gamma)}{\gamma - p} \tag{A37}$$

Therefore,

$$\frac{p\gamma M\_1 \mathbf{x} [1 - (1 - p)(1 - \gamma)\mathbf{x}]}{[1 - (1 - p)\mathbf{x}][1 - (1 - \gamma)\mathbf{x}]} = \frac{p\gamma M\_1 \mathbf{x}}{\gamma - p} \left( \frac{(1 - p)\gamma}{1 - (1 - p)\mathbf{x}} - \frac{p(1 - \gamma)}{1 - (1 - \gamma)\mathbf{x}} \right) \tag{A38}$$

For the second part of expression (A34), let

$$\begin{split} &\frac{1-(1-p)[2(1-\gamma)+p\gamma]\mathbf{x}+(1-p)^{2}(1-\gamma)^{2}\mathbf{x}^{2}}{[1-(1-p)\mathbf{x}][1-(1-\gamma)\mathbf{x}]^{2}} \\ &=\frac{A}{1-(1-p)\mathbf{x}}+\frac{B}{1-(1-\gamma)\mathbf{x}}+\frac{\mathbf{C}}{[1-(1-\gamma)\mathbf{x}]^{2}} \\ &=\frac{(A+B+\mathbf{C})-[(A+B)(1-\gamma)+A(1-\gamma)+B(1-p)+\mathbf{C}(1-p)]\mathbf{x}+c\_{2}\mathbf{x}^{2}}{[1-(1-p)\mathbf{x}][1-(1-\gamma)\mathbf{x}]^{2}} \end{split}$$

in which

$$c\_2 = [A(1 - \gamma) + B(1 - p)](1 - \gamma)$$

Thus, we have

$$\begin{cases} 1 = A + B + \mathbb{C} \\ (1 - p)[2(1 - \gamma) + p\gamma] = (A + B)(1 - \gamma) + A(1 - \gamma) + B(1 - p) + \mathbb{C}(1 - p) \\ (1 - p)^2(1 - \gamma) = A(1 - \gamma) + B(1 - p) \end{cases} \tag{A39}$$

Substituting the third relation and *A* + *B* = 1 − *C* into the second equation, we obtain

$$[(1-p)[2(1-\gamma)+p\gamma] = (1-\mathbb{C})(1-\gamma)+(1-p)^2(1-\gamma)+\mathbb{C}(1-p)\tag{A40}$$

from which we can solve that *C* = *p*.

Then, according to the equations

$$A + B = 1 - p \quad \text{and} \quad A(1 - \gamma) + B(1 - p) = (1 - p)^2(1 - \gamma)$$

the other two numbers are obtained to be

$$A = \frac{(1-p)^2 \gamma}{\gamma - p}, \quad B = -\frac{p(1-p)(1-\gamma)}{\gamma - p} \tag{A41}$$

Thus, the factorization of the second part is obtained.

The last part—that is, the function *h* (*m*) 2 (*x*)—is dealt with similarly. Omitting the straight-forward calculations, we directly find that

$$h\_2^{(m)}(\mathbf{x}) = -\frac{\eta(1-p)(p+\gamma-p\gamma)}{(1-\eta)[1-(1-p)(1-\gamma)\mathbf{x}]} + \frac{(1-p\eta)[\gamma+p(1-\gamma)\eta]}{(1-\eta)[1-(1-\gamma)(1-p\eta)\mathbf{x}]} \tag{A42}$$

Notice that (1 − *η*) is contained in the denominator of fractions in Equation (A42); thus, *η* = 1. When *η* = 1, Equation (29) shows that

$$h\_2^{(m)}(\mathbf{x})\Big|\_{\eta=1} = \frac{\gamma + p^2(1-\gamma) - (1-p)^2(1-\gamma)\gamma\mathbf{x}}{[1 - (1-p)(1-\gamma)\mathbf{x}]^2} \tag{A43}$$

Summarizing the above results, both the Equations (41) and (42) are determined.

#### **Appendix C. Factorization of Last Part of Equation (42)**

We write

$$\begin{split} &\frac{p\gamma M\_{2}\mathbf{x}^{2}}{[1-(1-\gamma)\mathbf{x}]^{2}} \cdot \frac{\gamma+p^{2}(1-\gamma)-(1-p)^{2}(1-\gamma)\gamma\mathbf{x}}{[1-(1-p)(1-\gamma)\mathbf{x}]^{2}} \\ &=p\gamma M\_{2}\mathbf{x}^{2}\bigg(\frac{A}{1-(1-\gamma)\mathbf{x}} + \frac{B}{[1-(1-\gamma)\mathbf{x}]^{2}} + \frac{\mathbb{C}}{1-(1-p)(1-\gamma)\mathbf{x}} + \frac{D}{[1-(1-p)(1-\gamma)\mathbf{x}]^{2}}\bigg) \\ &=p\gamma M\_{2}\mathbf{x}^{2}\bigg(\frac{(A+B)-A(1-\gamma)\mathbf{x}}{[1-(1-\gamma)\mathbf{x}]^{2}} + \frac{(\mathbb{C}+D)-\mathbb{C}(1-p)(1-\gamma)\mathbf{x}}{[1-(1-p)(1-\gamma)\mathbf{x}]^{2}}\bigg) \\ &\tag{A44} \end{split} \tag{A44}$$

Merging the terms in the bracket of (A44), according to corresponding coefficients, it is shown that

$$\begin{cases} A + B + \mathbb{C} + D = \gamma + p^2(1 - \gamma) \\ 2(A + B)(1 - p) + 2(\mathbb{C} + D) + A + \mathbb{C}(1 - p) = (1 - p)^2 \gamma \\ (A + B)(1 - p)^2 + 2A(1 - p) + (\mathbb{C} + D) + 2\mathbb{C}(1 - p) = 0 \\ A(1 - p) + \mathbb{C} = 0 \end{cases} \tag{A45}$$

The last row of (A45) shows that *C* = −*A*(1 − *p*), and using the first relationship, the second and the third row of Equation (A45) are equivalent to

$$\begin{cases} 2(1-p)[\gamma + p^2(1-\gamma) - (\mathbb{C}+D)] + 2(\mathbb{C}+D) + A - A(1-p)^2 = (1-p)^2\gamma\\ (1-p)^2[\gamma + p^2(1-\gamma) - (\mathbb{C}+D)] + 2A(1-p) + (\mathbb{C}+D) - 2A(1-p)^2 = 0 \end{cases} \quad \text{(A46)}$$

which gives

$$p(2-p)(\mathbb{C}+D) = -2p(1-p)A - (1-p)^2[\gamma + p^2(1-\gamma)]\tag{A47}$$

and

$$2p(\mathbb{C} + D) = -p(2 - p)A + (1 - p)^2 \gamma - 2(1 - p)[\gamma + p^2(1 - \gamma)] \tag{A48}$$

Combining Equations (A47) and (A48), the coefficient *A* is solved as

$$A = \frac{2 - p}{p^3} \left( (1 - p)^2 \gamma - \frac{2(1 - p)[\gamma + p^2(1 - \gamma)]}{2 - p} \right) \tag{A49}$$

and immediately

$$\mathcal{C} = -A(1-p) = -\frac{(1-p)(2-p)}{p^3} \left( (1-p)^2 \gamma - \frac{2(1-p)[\gamma + p^2(1-\gamma)]}{2-p} \right) \tag{A50}$$

Using Equation (A47), we have

$$\begin{split} D &= \frac{-2p(1-p)A - (1-p)^2[\gamma + p^2(1-\gamma)]}{p(2-p)} - C \\ &= \frac{-2p(1-p)A - (1-p)^2[\gamma + p^2(1-\gamma)]}{p(2-p)} + A(1-p) \\ &= \frac{p^2(1-p)A - (1-p)^2[\gamma + p^2(1-\gamma)]}{p(2-p)} \\ \end{split} \tag{A51}$$

and in the end, the last number *B* is determined by

$$B = \left[\gamma + p^2(1 - \gamma)\right] - A - \mathcal{C} - D \tag{A52}$$

So far, all the coefficients are obtained and the decomposition is totally determined.

#### **References**

## *Article* **Using Timeliness in Tracking Infections †**

**Melih Bastopcu 1 and Sennur Ulukus 2,\***

**Abstract:** We consider real-time timely tracking of infection status (e.g., COVID-19) of individuals in a population. In this work, a health care provider wants to detect both infected people and people who have recovered from the disease as quickly as possible. In order to measure the timeliness of the tracking process, we use the long-term average difference between the actual infection status of the people and their real-time estimate by the health care provider based on the most recent test results. We first find an analytical expression for this average difference for given test rates, infection rates and recovery rates of people. Next, we propose an alternating minimization-based algorithm to find the test rates that minimize the average difference. We observe that if the total test rate is limited, instead of testing all members of the population equally, only a portion of the population may be tested in unequal rates calculated based on their infection and recovery rates. Next, we characterize the average difference when the test measurements are erroneous (i.e., noisy). Further, we consider the case where the infection status of individuals may be dependent, which occurs when an infected person spreads the disease to another person if they are not detected and isolated by the health care provider. In addition, we consider an age of incorrect information-based error metric where the staleness metric increases linearly over time as long as the health care provider does not detect the changes in the infection status of the people. Through extensive numerical results, we observe that increasing the total test rate helps track the infection status better. In addition, an increased population size increases diversity of people with different infection and recovery rates, which may be exploited to spend testing capacity more efficiently, thereby improving the system performance. Depending on the health care provider's preferences, test rate allocation can be adjusted to detect either the infected people or the recovered people more quickly. In order to combat any errors in the test, it may be more advantageous for the health care provider to not test everyone, and instead, apply additional tests to a selected portion of the population. In the case of people with dependent infection status, as we increase the total test rate, the health care provider detects the infected people more quickly, and thus, the average time that a person stays infected decreases. Finally, the error metric needs to be chosen carefully to meet the priorities of the health care provider, as the error metric used greatly influences who will be tested and at what test rate.

**Keywords:** timely infection tracking; age of information; timely tracking of multiple processes; Markovian infection spread model

#### **1. Introduction**

We consider the problem of timely tracking of an infectious disease, e.g., COVID-19, in a population of *n* people. In this problem, a health care provider wants to detect infected people as quickly as possible in order to take precautions such as isolating them from the rest of the population. The health care provider also wants to detect people who have recovered from the disease as soon as possible since these people need to return to work which is especially critical in sectors such as education, food retail, public transportation,

**Citation:** Bastopcu, M.; Ulukus, S. Using Timeliness in Tracking Infections. *Entropy* **2022**, *24*, 779. https://doi.org/10.3390/e24060779

Academic Editor: Luca Faes

Received: 8 March 2022 Accepted: 27 May 2022 Published: 31 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

etc. Ideally, the health care provider should test all people all the time. However, as the total test rate is limited, the question is how frequently the health care provider should apply tests on these people when their infection and recovery rates are known. In a broader sense, this problem is related to timely tracking of multiple processes in a resource-constrained setting where each process takes binary values of 0 and 1 with different change rates.

Recent studies have shown that people who have recovered from infectious diseases such as COVID-19 can be reinfected. Furthermore, the recovery times of individuals may vary significantly. For these reasons, in this problem, the *i*th person becomes infected with rate *λi* which is independent of the others. Similarly, the *i*th person recovers from the disease with rate *μi*. We note that the index *i* may represent a specific individual or a group of individuals that share common features such as age, gender, and profession. Depending on the demographics, coefficients *λi* and *μi* may be statistically known by the health care provider. We denote the infection status of the *i*th person as *xi*(*t*) (shown with the black curves on the left in Figure 1) which takes the value 1 when the person is infected and the value 0 when the person is healthy. The health care provider applies tests to people marked as healthy with rate *si* and to people marked as infected with rate *ci*. Based on the test results, the health care provider forms an estimate for the infection status of the *i*th person denoted by *x*ˆ*i*(*t*) (shown with the blue curves on the right in Figure 1) which takes the value 1 when the most recent test result is positive and the value 0 when it is negative.

**Figure 1.** System model. There are *n* people whose infection status are given by *xi*(*t*). The health care provider applies tests on these people. Based on the test results, estimations for the infection status *x*ˆ*i*(*t*) are generated. Infected people are shown in red and healthy people are shown in green.

We measure the timeliness of the tracking process by the difference between the actual infection status of people and the real-time estimate of the health care provider which is based on the most recent test results. The difference can occur in two different cases: (i) when the person is sick (*xi*(*t*) = 1) and the health care provider maps this person as healthy (*x*ˆ*i*(*t*) = 0), and (ii) when the person recovers from the disease (*xi*(*t*) = 0) but the health care provider still considers this person as infected (*x*ˆ*i*(*t*) = 1). The former case represents the error due to late detection of infected people, while the latter case represents the error due to late detection of healed people. Depending on the health care provider's preferences, detecting infected people may be more important than detecting recovered people (controlling infection), or the other way around (returning people to workforce).

The age of information was proposed to measure timeliness of information in communication systems, and has been studied in the context of queueing systems [1–8], multi-hop and multi-cast networks [9–17], social networks [18], timely remote estimation of random processes [19–25], energy harvesting systems [26–40], wireless fading channels [41,42], scheduling in networks [43–55], lossless and lossy source and channel coding [56–66], vehicular, IoT and UAV systems [67–70], caching systems [71–82], computation-intensive systems [83–90], learning systems [91–93], gossip networks [94–97] and so forth. A more

detailed review of the age of information literature can be found in references [98–100]. Most relevant to our work, the real-time timely estimation of single and multiple counting processes [19,25], a Wiener process [20], a random walk process [101], and binary and multiple states Markov sources [23,51,102] have been studied. The study that is closest to our work is reference [23], where the remote estimation of a symmetric binary Markov source is studied in a time-slotted system by finding the optimal sampling policies via formulating a Markov Decision Process (MDP) for real-time error, AoI and AoII metrics. Different from [23], in our work, we consider real-time timely estimation of multiple nonsymmetric binary sources for a continuous time system. In our work, the sampler (health care provider) does not know the states of the sources (infection status of people), and thus, takes the samples (applies medical tests) randomly (exponential random variables) with fixed rates. Thus, we optimize the test rates of people to minimize the real-time estimation error.

In this paper, we consider the real-time timely tracking of infection status of *n* people. We first find an analytical expression for the long-term average difference between the actual infection status of people and the estimate of the health care provider based on test results. Then, we propose an alternating minimization-based algorithm to identify the test rates *si* and *ci* for all people. We observe that if the total test rate is limited, we may not apply tests on all people equally. Next, we provide an alternative method to characterize the average difference, by finding the steady state of a Markov chain defined by (*xi*(*t*), *x*ˆ*i*(*t*)). By using this alternative method, we determine the average estimation error when there are errors in the test measurements expressed by a false positive rate *p* and a false negative rate *q*. Next, we consider the infection status of two people where an infected person may spread the disease to another person if the infection has not been detected by the health care provider to consequently isolate the infected person. Finally, we consider an age of incorrect information-based error metric where the estimation error increases linearly over time when the health care provider has not detected the changes in the infection status of the people.

Through extensive numerical results, we observe that increasing the total test rate helps track the infection status of people better, and increasing the size of the population increases diversity which may be exploited to improve the performance. Depending on the health care provider's priorities, we can allocate additional tests to people marked as healthy to detect the infections faster or to people marked as infected to detect the recoveries more quickly. In order to combat the test errors, the health care provider may prefer to apply tests to only a selected portion of the population with higher test rates. When the infection status of a person depends on that of another person, the average time that a person remains infected can be reduced by increasing the total test rate as it helps to detect the infected people more quickly. Finally, we observe that depending on the error metric used, the test rate distribution among the population differs greatly, and thus, we should choose an error metric that aligns with the priorities of the health care provider.

#### **2. System Model**

We consider a population of *n* people. We denote the infection status of the *i*th person at time *t* as *xi*(*t*) (black curve in Figure 2a) which takes binary values 0 or 1 as follows,

$$x\_i(t) = \begin{cases} 1, & \text{if the } i \text{th person is infected at time } t, \\ 0, & \text{otherwise.} \end{cases} \tag{1}$$

In this paper, we consider a model where each person can be infected multiple times after recovering from the disease. We denote the time interval that the *i*th person stays healthy for the *j*th time as *Wi*(*j*) which is exponentially distributed with rate *λi*. We denote the recovery time for the *i*th person after being infected with the virus for the *j*th time as *Ri*(*j*) which is exponentially distributed with rate *μi*.

A health care provider wants to track the infection status of each person. Based on the test results at times *ti*,, the health care provider generates an estimate for the status of the *i*th person denoted as *x*ˆ*i*(*t*) (blue curve in Figure 2a) by

$$
\hat{x}\_i(t) = x\_i(t\_{i,\ell}), \quad t\_{i,\ell} \le t < t\_{i,\ell+1}. \tag{2}
$$

When *x*ˆ*i*(*t*) is 1, the health care provider applies the next test to the *i*th person after an exponentially distributed time with rate *ci*. When *x*ˆ*i*(*t*) is 0, the next test is applied to the *i*th person after an exponentially distributed time with rate *si*.

**Figure 2.** (**a**) A sample evolution of *xi*(*t*) and *x*ˆ*i*(*t*), and (**b**) the corresponding Δ*i*(*t*) in (5). Green areas correspond to the error caused by Δ*i*1(*t*) in (3). Orange areas correspond to the error caused by Δ*i*2(*t*) in (4).

An estimation error happens when the actual infection status of the *i*th person, *xi*(*t*), is different than the estimate of the health care provider, *x*ˆ*i*(*t*), at time *t*. This could happen in two ways: when *xi*(*t*) = 1 and *x*ˆ*i*(*t*) = 0, i.e., when the *i*th person is sick, but remains undetected by the health care provider, and when *xi*(*t*) = 0 and *x*ˆ*i*(*t*) = 1, i.e., when the *i*th person has recovered, but the health care provider is unaware that the *i*th person has recovered.

We denote the error caused by the former case, i.e., when *xi*(*t*) = 1 and *x*ˆ*i*(*t*) = 0, by Δ*i*1(*t*) (green areas in Figure 2b),

$$\Delta\_{i1}(t) = \max\{\mathfrak{x}\_i(t) - \mathfrak{x}\_i(t), 0\},\tag{3}$$

and we denote the error caused by the latter case, i.e., when *xi*(*t*) = 0 and *x*ˆ*i*(*t*) = 1, by Δ*i*2(*t*) (orange areas in Figure 2b),

$$\Delta\_{i2}(t) = \max\{\pounds\_i(t) - x\_i(t), 0\}. \tag{4}$$

Then, the total estimation error for the *i*th person Δ*i*(*t*) is

$$
\Delta\_i(t) = \theta \Delta\_{i1}(t) + (1 - \theta)\Delta\_{i2}(t), \tag{5}
$$

where *θ* is the importance factor in [0, 1]. A large *θ* gives more importance to the detection of infected people, and a small *θ* gives more importance to the detection of recovered people. We define the long-term weighted average difference between *xi*(*t*) and *x*ˆ*i*(*t*) as

$$\Delta\_{\bar{i}} = \lim\_{T \to \infty} \frac{1}{T} \int\_0^T \Delta\_{\bar{i}}(t)dt. \tag{6}$$

Then, the overall average difference of all people Δ is

$$
\Delta = \frac{1}{n} \sum\_{i=1}^{n} \Delta\_i. \tag{7}
$$

Our aim is to track the infection status of all people. Due to limited resources, there is a total test rate constraint ∑*n i*=1 *si* + ∑*n i*=1 *ci* ≤ *C*. Thus, our aim is to find the optimal test rates *si* and *ci* to minimize Δ in (7) while satisfying this total test rate constraint. We formulate the following problem,

$$\begin{array}{ll}\min & \Delta\\ \{s\_i, c\_i\} & \\ & \text{s.t.} & \sum\_{i=1}^n s\_i + \sum\_{i=1}^n c\_i \le \mathbb{C} \\ & s\_i \ge 0, \quad c\_i \ge 0, \quad i = 1, \dots, n. \end{array} \tag{8}$$

We provide a summary of the list of the variables used in this work in Table 1. In the next section, we find the total average difference Δ.

**Table 1.** List of variables used in this work.

#### **3. Average Difference Analysis**

In this section, we provide a probabilistic analysis to characterize the average difference Δ. In Section 5.1, we give an alternative method to find Δ by analyzing the steady-state distribution of the Markov chain induced by the states (*xi*(*t*), *x*ˆ*i*(*t*)). Here, we first find analytical expressions for Δ*i*1(*t*) in (3) and Δ*i*2(*t*) in (4) when *si* > 0 and *ci* > 0. We note that Δ*i*1(*t*) can be equal to 1 when *x*ˆ*i*(*t*) = 0 and is always equal to 0 when *x*ˆ*i*(*t*) = 1. Assume that at time 0, both *xi*(0) and *x*ˆ*i*(0) are 0. After an exponentially distributed time with rate *λi*, which is denoted by *Wi*, the *i*th person is infected, and thus *xi*(*t*) becomes 1. At that time, since *x*ˆ*i*(*t*) = 0, Δ*i*1(*t*) becomes 1. Further, Δ*i*1(*t*) will be equal to 0 again either when the *i*th person recovers from the disease which happens after *Ri* which is exponentially distributed with rate *μi* or when the health care provider performs a test on the *i*th person after *Di*, which is exponentially distributed with rate *si*. We define *Tm*(*i*) as the earliest time at which one of these two cases happens, i.e., *Tm*(*i*) = min{*Ri*, *Di*} (which is shown by the green areas in Figure 3a). We note that *Tm*(*i*) is also exponentially distributed with rate *μi* + *si*, and we have P(*Tm*(*i*) = *Ri*) = *μi μi*+*si* and P(*Tm*(*i*) = *Di*) = *si μi*+*si* . If the *i*th person recovers from the disease before testing, we return to the initial case where both *xi*(*t*) and *x*ˆ*i*(*t*) are equal to 0 again. In this case, the cycle repeats itself, i.e., the *i*th person becomes sick again after *Wi* and Δ*i*1(*t*) remains as 1 until either the person recovers or the health care provider performs a test which takes another *Tm*(*i*) duration. If the health care provider performs a test before the person recovers, then *x*ˆ*i*(*t*) becomes 1. We denote the time interval for which *x*ˆ*i*(*t*) stays at 0 as *Ii*1 which is given by

$$I\_{i1} = \sum\_{\ell=1}^{K\_1} T\_{\mathfrak{m}}(i,\ell) + \mathcal{W}\_i(\ell), \tag{9}$$

where *K*1 is geometric with rate P(*Tm*(*i*) = *Di*) = *si μi*+*si* . Due to [103] (Prob. 9.4.1), ∑*K*1 =1 *Tm*(*i*, ) and ∑*K*1 =1 *Wi*() are exponentially distributed with rates *si* and *λisi μi*+*si* , respectively. As E[*Ii*1] = E[∑*K*1 =1 *Tm*(*i*, )] + E[∑*K*1 =1 *Wi*()], we have

$$\mathbb{E}[I\_{i1}] = \frac{1}{s\_i} + \frac{s\_i + \mu\_i}{s\_i \lambda\_i}. \tag{10}$$

When *x*ˆ*i*(*t*) = 1, the health care provider marks the *i*th person as infected. The *i*th person recovers from the virus after *Ri*. After the *i*th person recovers, either the health care provider performs a test after *Zi* which is exponentially distributed with rate *ci* or the *i*th person is reinfected with the virus which takes *Wi* time. We define *Tu*(*i*) as the earliest time at which one of these two cases happens, i.e., *Tu*(*i*) = min{*Wi*, *Zi*} (which is shown by the orange areas in Figure 3b). Similarly, we note that *Tu*(*i*) is exponentially distributed with rate *λi* + *ci*, and we have P(*Tu*(*i*) = *Wi*) = *λi λi*+*ci* and P(*Tu*(*i*) = *Zi*) = *ci λi*+*ci* . If the person is reinfected with the virus before a test is applied, this cycle repeats itself, i.e., the *i*th person recovers after another *Ri*, and then either a test is applied to the *i*th person, or the person is infected again which takes another *Tu*(*i*). If the health care provider performs a test to the *i*th person before the person is reinfected, the health care provider marks the *i*th person as healthy again, i.e., *x*ˆ*i*(*t*) becomes 0. We denote the time interval that *x*ˆ*i*(*t*) is equal to 1 as *Ii*2 which is given by

$$I\_{i2} = \sum\_{\ell=1}^{K\_2} T\_u(i,\ell) + R\_i(\ell),\tag{11}$$

where *K*2 is geometric with rate P(*Tu*(*i*) = *Zi*) = *ci λi*+*ci* . Similarly, ∑*K*2 =1 *Tu*(*i*, ) and ∑*K*2 =1 *Ri*() are exponentially distributed with rates *ci* and *ciμi λi*+*ci* , respectively. As E[*Ii*2] = E[∑*K*2 =1 *Tu*(*i*, )] + E[∑*K*2 =1 *Ri*()], we have

$$\mathbb{E}[I\_{i2}] = \frac{1}{c\_i} + \frac{c\_i + \lambda\_i}{c\_i \mu\_i}. \tag{12}$$

We denote the time interval between the *j*th and (*j* + 1)th times that *x*ˆ*i*(*t*) changes from 1 to 0 as the *j*th cycle *Ii*(*j*) where *Ii*(*j*) = *Ii*1(*j*) + *Ii*2(*j*). We note that Δ*i*1(*t*) is always equal to 0 during *Ii*2(*j*), i.e., *x*ˆ*i*(*t*) = 1, and Δ*i*1(*t*) is equal to 1 when *xi*(*t*) = 1 in *Ii*1(*j*). We denote the total time duration when Δ*i*1(*t*) is equal to 1 as *Te*,1(*i*, *j*) during the *j*th cycle where *Te*,1(*i*, *j*) = ∑*K*1 =1 *Tm*(*i*, ). Thus, we have E[*Te*,1(*i*)] = 1 *si* . Then, using ergodicity, similar to [80], Δ*i*1 is equal to

$$\Delta\_{i1} = \frac{\mathbb{E}[T\_{\mathfrak{e},1}(i)]}{\mathbb{E}[I\_{i}]} = \frac{\mathbb{E}[T\_{\mathfrak{e},1}(i)]}{\mathbb{E}[I\_{i1}] + \mathbb{E}[I\_{i2}]}.\tag{1.3}$$

Thus, we have

$$
\Delta\_{i1} = \frac{\mu\_i \lambda\_i}{\mu\_i + \lambda\_i} \frac{c\_i}{\mu\_i c\_i + \lambda\_i s\_i + c\_i s\_i}. \tag{14}
$$

Next, we find Δ*i*2. We note that Δ*i*2(*t*) is equal to 1 when *xi*(*t*) = 0 in *Ii*2(*j*) and is always equal to 0 during *Ii*1(*j*). Similarly, we denote the total time duration where Δ*i*2(*t*) is equal to 1 in the *j*th cycle *Ii*(*j*) as *Te*,2(*i*, *j*) which is equal to *Te*,2(*i*, *j*) = ∑*K*2 =1 *Tu*(*i*, ). Thus, we have E[*Te*,2(*i*)] = 1 *ci* . Then, similar to Δ*i*1 in (13), Δ*i*2 is equal to

$$
\Delta\_{i2} = \frac{\mu\_i \lambda\_i}{\mu\_i + \lambda\_i} \frac{s\_i}{\mu\_i c\_i + \lambda\_i s\_i + c\_i s\_i}. \tag{15}
$$

By using (5), (14), and (15), we obtain Δ*i* as

$$
\Delta\_i = \frac{\mu\_i \lambda\_i}{\mu\_i + \lambda\_i} \frac{\theta c\_i + (1 - \theta) s\_i}{\mu\_i c\_i + \lambda\_i s\_i + c\_i s\_i}. \tag{16}
$$

Then, by inserting (16) in (7), we obtain Δ. In the next section, we solve the optimization problem in (8).

**Figure 3.** A sample evolution of (**a**) Δ*i*1(*t*), and (**b**) Δ*i*2(*t*) in a typical cycle.

#### **4. Optimization of Average Difference**

In this section, we solve the optimization problem in (8). Using Δ*i* in (16) in (7), we rewrite (8) as

$$\begin{aligned} \min\_{\{s\_i, c\_i\}} & \quad \sum\_{i=1}^n \frac{\mu\_i \lambda\_i}{\mu\_i + \lambda\_i} \frac{\theta c\_i + (1 - \theta) s\_i}{\mu\_i c\_i + \lambda\_i s\_i + c\_i s\_i} \\ \text{s.t.} & \quad \sum\_{i=1}^n s\_i + \sum\_{i=1}^n c\_i \le \mathbb{C} \\ & \quad s\_i \ge 0, \quad c\_i \ge 0, \quad i = 1, \dots, n. \end{aligned}$$

We define the Lagrangian function [104] for (17) as

$$\mathcal{L} = \sum\_{i=1}^{n} \frac{\mu\_i \lambda\_i}{\mu\_i + \lambda\_i} \frac{\theta c\_i + (1 - \theta) s\_i}{\mu\_i c\_i + \lambda\_i s\_i + c\_i s\_i} + \beta \left( \sum\_{i=1}^{n} s\_i + c\_i - \mathcal{C} \right) - \sum\_{i=1}^{n} \nu\_i s\_i - \sum\_{i=1}^{n} \eta\_i c\_i \tag{18}$$

where *β* ≥ 0, *νi* ≥ 0, and *ηi* ≥ 0. The KKT conditions are

$$\frac{\partial \mathcal{L}}{\partial \mathbf{s}\_{i}} = \frac{\mu\_{i} \lambda\_{i} c\_{i}}{\mu\_{i} + \lambda\_{i}} \frac{(1 - \theta)\mu\_{i} - \theta(c\_{i} + \lambda\_{i})}{(\mu\_{i} c\_{i} + \lambda\_{i} s\_{i} + s\_{i} c\_{i})^{2}} + \beta - \nu\_{i} = 0,\tag{19}$$

$$\frac{\partial \mathcal{L}}{\partial \mathbf{c}\_{i}} = \frac{\mu\_{i}\lambda\_{i}\mathbf{s}\_{i}}{\mu\_{i} + \lambda\_{i}} \frac{\theta \lambda\_{i} - (1 - \theta)(\mu\_{i} + s\_{i})}{(\mu\_{i}\mathbf{c}\_{i} + \lambda\_{i}\mathbf{s}\_{i} + s\_{i}\mathbf{c}\_{i})^{2}} + \beta - \eta\_{i} = 0,\tag{20}$$

for all *i*. The complementary slackness conditions are

$$
\beta \left( \sum\_{i=1}^{n} s\_i + c\_i - \mathbb{C} \right) = 0, \quad \nu\_i s\_i = 0, \quad \eta\_i c\_i = 0. \tag{21}
$$

First, we find *si*. From (19), we have

$$\left(\mu\_i c\_i + \lambda\_i s\_i + s\_i c\_i\right)^2 = \frac{\mu\_i \lambda\_i c\_i}{\mu\_i + \lambda\_i} \frac{\theta(c\_i + \lambda\_i) - (1 - \theta)\mu\_i}{\beta - \nu\_i}.\tag{22}$$

When *θ*(*ci* + *λi*) ≥ (1 − *θ*)*μi*, we solve (22) for *si* as

$$s\_i = \frac{\mu\_i c\_i}{\lambda\_i + c\_i} \left( \sqrt{\frac{1}{\mu\_i c\_i} \frac{\lambda\_i}{\mu\_i + \lambda\_i} \frac{\theta(c\_i + \lambda\_i) - (1 - \theta)\mu\_i}{\beta}} - 1 \right)^+,\tag{23}$$

where we used the fact that we either have *si* > 0 and *νi* = 0, or *si* = 0 and *νi* ≥ 0, due to (21). Here, (·)+ = max(·, 0). On the other hand, when *θ*(*ci* + *λi*) < (1 − *θ*)*μi*, we have *∂*Δ*i ∂si* > 0, and thus it is optimal to choose *si* = 0 as our aim is to minimize Δ in (7). In this case, when *si* = 0, we have Δ*i* = *θλi μi*+*λi* which is independent of the value of *ci*. As we obtain the same Δ*i* for all values of *ci*, and the total update rate is limited, i.e., ∑*n i*=1 *si* + *ci* ≤ *C*, in this case, it is optimal to choose *ci* = 0 as well (i.e., when *si* = 0).

Next, we find *ci*. From (20), we have

$$(\mu\_i c\_i + \lambda\_i s\_i + s\_i c\_i)^2 = \frac{\mu\_i \lambda\_i s\_i}{\mu\_i + \lambda\_i} \frac{(1 - \theta)(\mu\_i + s\_i) - \theta \lambda\_i}{\beta - \eta\_i}. \tag{24}$$

When (1 − *θ*)(*μi* + *si*) ≥ *θλi*, we solve (24) for *ci* as

$$c\_{i} = \frac{\lambda\_{i}s\_{i}}{\mu\_{i} + s\_{i}} \left( \sqrt{\frac{1}{\lambda\_{i}s\_{i}}} \frac{\mu\_{i}}{\mu\_{i} + \lambda\_{i}} \frac{(1-\theta)(s\_{i} + \mu\_{i}) - \theta\lambda\_{i}}{\beta}} - 1 \right)^{+},\tag{25}$$

where we used the fact that we either have *ci* > 0 and *ηi* = 0, or *ci* = 0 and *ηi* ≥ 0, due to (21). Similarly, when (1 − *θ*)(*si* + *μi*) < *θλi*, we have *∂*Δ*i ∂ci* > 0. Thus, in this case, it is optimal to choose *ci* = 0. When *ci* = 0, we have Δ*i* = (1−*θ*)*μi μi*+*λi* which is independent of the value of *si*. Thus, it is optimal to choose *si* = 0 when *ci* = 0.

From (23), if 1 *μici λi μi*+*λi* (*θ*(*ci* + *λi*) − (1 − *θ*)*μi*) ≤ *β*, we must have *si* = 0. Thus, for a given *ci*, the optimal test rate allocation policy for *si* is a *threshold policy* where *si*'s with small 1 *μici λi μi*+*λi* (*θ*(*ci* + *λi*) − (1 − *θ*)*μi*) are equal to zero. Similarly, from (25), if 1 *λisi μi μi*+*λi* ((1 − *θ*)(*si* + *μi*) − *θλi*) ≤ *β*, we must have *ci* = 0. Thus, for a given *si*, the optimal policy to determine *ci* is a *threshold policy* where *ci*'s with small 1 *λisi μi μi*+*λi* ((1 − *θ*)(*si* + *μi*) − *θλi*) are equal to zero.

Next, we show that in the optimal policy, if *si* > 0 and *ci* > 0 for some *i*, then the total test rate constraint must be satisfied with equality, i.e., ∑*n i*=1 *si* + *ci* = *C*.

**Lemma 1.** *In the optimal policy, if si* > 0 *and ci* > 0 *for some i, then we have* ∑*n i*=1 *si* + *ci* = *C.*

**Proof of Lemma 1.** The derivatives of Δ*i* with respect to *si* and *ci* are

$$\frac{\partial \Delta\_{\dot{i}}}{\partial \mathbf{s}\_{\dot{i}}} = \frac{\mu\_{i}\lambda\_{i}\mathbf{c}\_{\dot{i}}}{\mu\_{i} + \lambda\_{\dot{i}}} \frac{\left(1 - \theta\right)\mu\_{i} - \theta\left(\mathbf{c}\_{\dot{i}} + \lambda\_{\dot{i}}\right)}{\left(\mathbf{c}\_{i}\mu\_{i} + \mathbf{s}\_{i}\mathbf{c}\_{\dot{i}} + \lambda\_{i}\mathbf{s}\_{\dot{i}}\right)^{2}}\,\mathrm{}\tag{26}$$

$$\frac{\partial \Delta\_i}{\partial c\_i} = \frac{\mu\_i \lambda\_i s\_i}{\mu\_i + \lambda\_i} \frac{\theta \lambda\_i - (1 - \theta)(s\_i + \mu\_i)}{\left(c\_i \mu\_i + s\_i c\_i + \lambda\_i s\_i\right)^2}. \tag{27}$$

We note that *si* > 0 in (23) implies that *θ*(*ci* + *λi*) > (1 − *θ*)*μi*. In this case, we have *∂*Δ*i ∂si* < 0. Similarly, *ci* > 0 in (25) implies that (1 − *θ*)(*si* + *μi*) > *θλi*. Thus, we have *∂*Δ*i ∂ci* < 0. Therefore, in the optimal policy, if we have *si* > 0 and *ci* > 0 for some *i*, then we must have ∑*n i*=1 *si* + *ci* = *C*. Otherwise, we can further decrease Δ in (7) by increasing *ci* or *si*.

Next, we propose an alternating minimization-based algorithm for finding *si* and *ci*. For this purpose, for given initial (*si*, *ci*) pairs, we define *φi* as

$$\phi\_{i} = \begin{cases} \frac{1}{\mu\_{i}c\_{i}} \frac{\lambda\_{i}}{\mu\_{i} + \lambda\_{i}} (\theta(c\_{i} + \lambda\_{i}) - (1 - \theta)\mu\_{i}), \ i = 1, \ldots, n, \\\frac{1}{\lambda\_{i}s\_{i}} \frac{\mu\_{i}}{\mu\_{i} + \lambda\_{i}} ((1 - \theta)(s\_{i} + \mu\_{i}) - \theta\lambda\_{i}), \ i = n + 1, \ldots, 2n. \end{cases} \tag{28}$$

Then, we define *ui* as

$$u\_{i} = \begin{cases} \frac{\mu\_{i}c\_{i}}{\lambda\_{i} + c\_{i}} \left(\sqrt{\frac{\Phi\_{i}}{\beta}} - 1\right)^{+}, & i = 1, \ldots, n, \\\frac{\lambda\_{i}s\_{i}}{\mu\_{i} + s\_{i}} \left(\sqrt{\frac{\Phi\_{i}}{\beta}} - 1\right)^{+}, & i = n + 1, \ldots, 2n. \end{cases} \tag{29}$$

From (23) and (25), *si* = *ui* and *ci* = *un*+*i*, for *i* = 1, . . . , *n*.

Next, we find *si* and *ci* by determining *β* in (29). First, assume that, in the optimal policy, there is an *i* such that *si* > 0 and *ci* > 0. Thus, by Lemma 1, we must have ∑*n i*=1 *si* + *ci* = *C*. We initially take random (*si*, *ci*) pairs such that ∑*n i*=1 *si* + *ci* = *C*. Then, given the initial (*si*, *ci*) pairs, we immediately choose *ui* = 0 for *φi* < 0. For the remaining *ui* with *φi* ≥ 0, we apply a solution method similar to that in [80]. By assuming *φi* ≥ *β*, i.e., by disregarding (·)+ in (29), we solve ∑2*n i*=1 *ui* = *C* for *β*. Then, we compare the smallest *φi* which is larger than zero in (28) with *β*. If we have *φi* ≥ *β*, then it implies that *ui* ≥ 0 for all remaining *i*. Thus, we have obtained *ui* values for given initial (*si*, *ci*) pairs. If the smallest *φi* which is larger than zero is smaller than *β*, then the corresponding *ui* is negative and we should choose *ui* = 0 for the smallest non-negative *φi*. Then, we repeat this procedure until the smallest non-negative *φi* is larger than *β*. After determining all *ui*, we obtain *si* = *ui* and *ci* = *un*+*i* for *i* = 1, ... , *n*. Then, with the updated values of (*si*, *ci*) pairs, we keep finding *ui*'s until the KKT conditions in (19) and (20) are satisfied.