# Thermodynamical limits for models of car-sharing systems: the Autolib' example

Guy Fayolle, Christine Fricker

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

## To cite this version:

Guy Fayolle, Christine Fricker. Thermodynamical limits for models of car-sharing systems: the Autolib' example. Markov Processes And Related Fields, 2025, 31 (1), pp.55-78. hal-04855046v2

HAL Id: hal-04855046

https://inria.hal.science/hal-04855046v2

Submitted on 25 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Thermodynamical limits for models of car-sharing systems: the Autolib' example

Guy Fayolle* Christine Fricker †

January 2025

## Abstract

We analyze mean-field equations obtained for models motivated by a large station-based car-sharing system in France called *Autolib'.* The main focus is on a version where users reserve a parking space when they take a car. In a first model, the reservation of parking spaces is effective for all users (see [4]) and capacity constraints are ignored. The model is carried out in thermodynamical limit, that is when the number *N* of stations and the number of cars *M*<sub>*N*</sub> tend to infinity, with *U* = $\lim_{N \to \infty} M*_{*N}*/*N*. This limit is described by Kolmogorov's equations of a two-dimensional time-inhomogeneous Markov process depicting the numbers of reservations and cars at a station. It satisfies a non-linear differential system. We prove analytically that this system has a unique solution, which converges, as *t* → ∞, to an equilibrium point exponentially fast. Moreover, this equilibrium point corresponds to the stationary distribution of a two queue tandem (reservations, cars), which is here always ergodic. The intensity factor of each queue has an explicit form obtained from an intrinsic mass conservation relationship. Two related models with capacity constraints are briefly presented in the last section: the simplest one with no reservation leads to a one-dimensional problem; the second one corresponds to our first model with finite total capacity *K*.

**Keywords:** Markov Process, Queueing Systems, Thermodynamical limit, Mean-field, Car sharing.

## 1 Introduction

The paper investigates asymptotic properties for some time-inhomogeneous Markov processes obtained as mean-field limits. More precisely, in a system with *N* sites and *M*<sub>*N*</sub> particles, they describe the limiting behavior of an arbitrary site in the system, when the numbers *N* and *M*<sub>*N*</sub> get large together at the same speed, i.e. when *M*<sub>*N*</sub>/*N* tends to a constant *U* > 0 as *N* → ∞. The limit thus obtained is generally called the *thermodynamic limit* in statistical physics.

Such a Markov process is described by its transient probability given by the forward Kolmogorov's equations. We prove the existence and uniqueness of a solution to these non-linear differential equations, as well as an exponential rate of convergence toward the stationary regime. This is a main issue, which appears to be crucial in the literature to obtain further results on the mean-field convergence as, for example, the so-called *refined mean-field* in [13, 1].

This question has also been investigated for finite state reversible dynamics [6, 5], involving the construction of Lyapounov functions by using relative entropy. Our study goes beyond this framework with a simple model emanating from car-sharing systems (CSS) (see [4]), which can also be represented as systems of interacting particles. Here, stations are viewed as sites and cars become particles. The main feature of these station-based CSS (e.g. *Autolib'* ) is that users can reserve a parking slot at their destination.

<footnote>*INRIA Paris, 48 rue Barrault, 75013 Paris, France. Email: Guy.Fayolle@inria.fr</footnote>

<footnote>†INRIA Paris, 48 rue Barrault, 75013 Paris, France. E.mail: Christine.Fricker@inria.fr</footnote>
## 1.1 A particle model with reservation

Let us describe the model as an homogeneous version of simplified car-sharing dynamics. Consider $M_N$ particles randomly spread among $N$ sites, each site having a capacity $K \le +\infty$, so that $M_N \le KN$ when $K$ is finite. Particles can be of two types, either a reserved parking space (reservation) or a car. At time $t=0$, the repartition of the particles among the sites is arbitrary. From site $i$, after a random time exponentially distributed with parameter $\lambda$, a car decides to go to site $j$ with probability $1/N$ and a reservation appears at site $j$. If the reservation is possible, i.e. when the total number of particles at site $j$ is less than $K$, then the car leaves $i$, and the reservation at $j$ remains for a travel time exponentially distributed with parameter $\mu$ (the case $i=j$ is possible, with the same parameters). Otherwise, no car departure occurs.

It should be noted that another dynamics would be to allow the particle of type car to randomly search for an unsaturated site, but we shall not consider this model.

## 1.2 A queueing description of the model

We will recast the mean field system introduced above in a queueing context, and it will be useful to recall the classical notation, which can be found e.g. in [14].

$M/M/1/\infty$ stands for a single server queue, with Poisson input, exponentially distributed service times, infinite capacity. In this study, the service discipline can be FIFO (first-in-first-out), LIFO (last-in-first-out) or RAND (random order). When the Poisson intensity is time-dependent, we shall write $M(t)/M/1/\infty$; otherwise, we shall write $1/M/1/\infty$ for an unspecified arrival point process.

Similarly, $M/M/\infty$ (resp. $M(t)/M/\infty$) denotes a queue with infinitely many servers, Poisson input with constant (resp. time dependent) intensity and exponentially distributed service times.

In the queueing systems just defined, arrivals and service times are independent.

The model proposed in Section 1.1 can alternatively be viewed as a closed Jackson network (see [3]), with $M_N$ customers. This means that each of the $N$ sites can be replaced by a tandem with total capacity $K$. The customers in the first infinite-server queue of tandem $i$ represent reservations, and customers in the second $M/M/1/\infty$ queue represent available cars. The service times of the tandem are, respectively, the travel times of the cars and the inter-arrival times of the users. Both are independent and exponentially distributed, with respective parameters $\mu$ and $\lambda$. Note that the routing considered here involves a blocking policy. Indeed, if the destination site (tandem) is saturated, the user request is rejected.

## 1.3 Main outcomes

* Section 2 is devoted to a mean-field model of the type presented in Section 1.1, in which reservation of parking spaces is effective for all users (see [4]) and capacity constraints are ignored ($K = \infty$). We prove existence and uniqueness of the equations describing the evolution of an underlying dynamical system. The related stationary distribution has a product form corresponding to a Jackson network of two queues in tandem, where the intensity factor in each queue is obtained as the unique admissible root of a second degree equation, derived from a mass conservation relationship. In addition, the speed of convergence toward equilibrium is exponential.

* Section 3 contents the proof of the main theorem stated in Section 2.

* In Section 4, two related models with capacity constraints are presented: one without reservation, which leads to a single queue problem, the second being merely the model of Section 2 with finite $K$.

* The Appendix provides some technical properties, in particular stochastic dominance, concerning the dependence of time-inhomogeneous birth and death processes on the parameters.
## 2 Thermodynamical limit in the unbounded capacity case

In this section and in Section 3, we focus on the case $K = \infty$, referring to Section 4 for variants of the model with finite $K$.

Although system dynamics may seem simpler when $K = \infty$, its study is still far from complete. For related topics, we refer for instance to [8, 12].

For any fixed $N$, let us denote by $R_i^N(t)$ and $V_i^N(t)$ the number of reserved slots and cars at site $i$ ($1 \le i \le N$), at time $t \ge 0$. Then the process

$$ (R^N(t), V^N(t)) = \{(R_i^N(t), V_i^N(t)), 1 \le i \le N\} $$

is Markovian and irreducible on the finite state-space

$$ S^N = \{(r, v) = (r_i, v_i)_{1 \le i \le N} \in \mathbb{N}^{2N}, \sum_{i=1}^{N} (r_i + v_i) = M_N\}, $$

with a unique invariant probability measure. Its infinitesimal generator is given by

$$ G^N f(r, v) = \frac{\lambda}{N} \sum_{i,j=1}^{N} (f(r + e_j, v - e_i) - f(r, v)) 1_{v_i > 0} \\ + \mu \sum_{i=1}^{N} r_i (f(r - e_i, v + e_i) - f(r, v)), \quad \forall (r, v) \in S^N, $$

for any function $f(\cdot, \cdot)$ with support on $S^N$. The so-called mean-field approach is used to investigate the large-scale behavior of the system. We start from the empirical measure given, for any pair $(j, k) \in \mathbb{N}^2$ and $t \ge 0$, by

$$ \alpha_{j,k}^N(t) = \frac{1}{N} \sum_{i=1}^{N} 1_{\{(R_i^N(t), V_i^N(t))=(j,k)\}}. $$

With this definition, $\alpha_{j,k}^N(t)$ represents the proportion of tandems (sites) at time $t$, with $j$ and $k$ customers in their respective queues.

Let $\bar{\alpha}^N(t) \stackrel{\text{def}}{=} \{\alpha_{j,k}^N(t)\}_{(j,k) \in \mathbb{N}^2}$. Note that, since there are $M_N$ particles in the system as defined in Section 1.1, we have the following conservation equality

$$ \sum_{j=1}^{N} (j + k) \alpha_{j,k}^{N}(t) = \frac{M_N}{N}, \quad t \ge 0. $$

In the sequel, we are mainly interested in the limit behavior of the system as $N \to \infty$, assuming

$$ U \stackrel{\text{def}}{=} \lim_{N \to \infty} \frac{M_N}{N}. $$

Letting $\Lambda$ denote the set of probability measures $\{\rho_{j,k}, j, k \ge 0\}$ on $\mathbb{N}^2$, subject to the constraint

$$ \sum_{j,k \ge 0} (j+k) \rho_{j,k} = U, \qquad (2.1) $$

the following proposition holds.

**Proposition 2.1.** Assume $\bar{\alpha}^N(0)$ converges weakly, as $N \to \infty$, to some fixed distribution $\bar{\alpha}_0 \in \Lambda$. Then the empirical measure $\bar{\alpha}^N(t)$ converges in distribution to a deterministic dynamical system denoted by $\bar{\alpha}(t)$, which satisfies the following infinite system of nonlinear forward Kolmogorov's equations

$$ \frac{d\alpha_{j,k}(t)}{dt} + [\lambda b(t) + \mu j + \lambda \mathbb{1}_{\{k>0\}}] \alpha_{j,k}(t) = \\ \lambda b(t) \alpha_{j-1,k}(t) \mathbb{1}_{\{k>0\}} + \mu (j+1) \alpha_{j+1,k-1}(t) \mathbb{1}_{\{k>0\}} + \lambda \alpha_{j,k+1}(t), \quad (j,k) \in \mathbb{N}^2, t \ge 0, \quad (2.2) $$
where $\vec{\alpha}(t) \in \Lambda$, $\vec{\alpha}(0) = \vec{\alpha}_0$ and $b(t) \stackrel{\text{def}}{=} 1 - \sum_{j \in \mathbb{N}} \alpha_{j,0}(t)$.

*Proof.* This result was demonstrated in [11]. It could also be proved along the lines proposed in [8], which heavily rely on standard theoretical tools involving the convergence of generators (see [9]).

**Remark 2.2.** It will be shown in Section 3.1 (see (3.2)), that the mass conservation equation

$$ \sum_{j,k \ge 0} (j+k)\alpha_{j,k}(t) = U, \quad \forall t \ge 0, $$

is in fact intrinsic to the system (2.2).

## 2.1 Mean-field equations viewed as a tandem of two queues

System 2.2 represents the joint distribution of the number of units in the tandem shown in Fig. 2.1. The first queue stands for the reservations and is of $M(t)/M/\infty$ type, with arrival (resp. service) rate $\lambda b(t)$ (resp. $\mu$), while the second queue containing the cars is of $./M/1/\infty$ type, with FIFO discipline and service rate $\lambda$. Here $b(t)$ is the probability that there is at least one car in the second queue.

<figure><img src="image_3.png" /><figcaption>Fig. 2.1: Dynamics of the tandem. The first queue is of $M(t)/M/\infty$ type and the second queue is a simple $./M/1/\infty$ queue.</figcaption></figure>

Introduce the generating function

$$ F(x, y, t) \stackrel{\text{def}}{=} \sum_{j,k \ge 0} \alpha_{jk}(t) x^j y^k, $$

where $x, y$ are complex variables such that $0 \le |x|, |y| \le 1$. Then, (2.2) is equivalent to the nonlinear functional equation

$$ \frac{\partial F(x, y, t)}{\partial t} + \left[ \lambda b(t)(1-x) + \lambda \left(1 - \frac{1}{y}\right) \right] F(x, y, t) = \mu(y-x) \frac{\partial F(x, y, t)}{\partial x} + \lambda \left(1 - \frac{1}{y}\right) F(x, 0, t), \quad (2.3) $$

where $b(t) = 1 - F(1, 0, t)$.

## 2.2 Results

Equation (2.2) can be viewed as the limit of mean-field equations subject to the capacity constraint $j+k \le K$ as $K \to \infty$. In this case, we will show that the solution of (2.2) always exist, and that the stationary regime is reached exponentially fast, for any $\lambda, \mu > 0$. The global situation is pictured in the next theorem.

**Theorem 2.3.**

(i) For each $t \ge 0$ there exists a unique distribution $\vec{\alpha}(t)$ satisfying the nonlinear system (2.2). In addition, for any $\lambda, \mu > 0$, as $t \to \infty$, the distribution $\vec{\alpha}(t)$ converges in $\Lambda$ (defined by (2.1)) to the unique equilibrium point $\vec{\pi} > \vec{0}$ given by

$$ \pi_{j,k} = e^{-\rho\beta} \frac{(\rho\beta)^j}{j!} (1-\beta)^{\beta^k}, \quad \forall j, k \in \mathbb{N}, \quad (2.4) $$
where

$$ \rho \stackrel{\text{def}}{=} \frac{\lambda}{\mu}, \quad \text{and} \quad \beta \stackrel{\text{def}}{=} \lim_{t \to \infty} b(t) = \frac{U + \rho + 1 - \sqrt{(U + \rho + 1)^2 - 4\rho U}}{2\rho} < 1, \qquad (2.5) $$

$U = \lim_{N \to \infty} M_N/N$ being the constant introduced in Section 1.

(ii) There is a mass conservation of the form

$$ U = \sum_{j,k \ge 0} (j+k) \alpha_{j,k}(t), \quad 0 \le t \le \infty, \qquad (2.6) $$

which, although somehow surprising, follows implicitly from the functional equation (2.3).

(iii) Assuming at $t = 0$ all particles are reserved parking spaces, i.e. $\bar{\alpha}(0)$ has its support in $\mathbb{N} \times \{0\}$, then, as $t \to \infty$,

$$ b(t) = \beta + \mathcal{O}(e^{-vt}), \qquad (2.7) $$

where

$$ v = \min[\mu, \lambda(\sqrt{1 - \sqrt{\beta}})^2], \qquad (2.8) $$

and $1/v$ being viewed as the relaxation time of the system.

This theorem shows in particular that the stationary distribution of (2.3) is given by (3.17), which has a product-form solution. This corresponds to the network shown in Figure 2.1, for $t = \infty$, with constant rates, whose joint-distribution is well-known to have a product-form (see [3]).

# 3 Proof of Theorem 2.3

The line of argument (both of analytic and probabilist nature) consists in viewing the second queue of our tandem as a $M(t)/M/1/\infty$ queue with service rate $\lambda$. Its arrival rate $\delta(t)$ defined in (3.6) is theoretically computable via system 3.12, and directly related to $b(t)$ introduced in (2.3). In a second step, we show via an iterative scheme the convergence of $\delta(t)$ to a finite constant $\delta$, whose explicit form is given in Lemma 3.6. The speed of convergence is obtained in Section 3.3.

## 3.1 Subsidiary equations

Beforehand, we establish three useful functional equations which follow immediately from (2.3), keeping in mind that $x, y$ are complex variables with $|x|, |y| \le 1$.

### 3.1.1 Equation for $F(x, x, t)$

Setting $G(x, t) \stackrel{\text{def}}{=} F(x, x, t)$, with $G(1, t) = 1$, equation (2.3) yields

$$ \frac{\partial G(x, t)}{\partial t} + \lambda \left(b(t) - \frac{1}{x}\right) (1-x) G(x, t) = \lambda \left(1 - \frac{1}{x}\right) F(x, 0, t). \qquad (3.1) $$

Then, dividing equation (3.1) by $1-x$ and letting $x \to 1$, we get

$$ \frac{\partial}{\partial t} \frac{\partial}{\partial x} G(1, t) = 0, $$

or, equivalently,

$$ \frac{\partial G(1, t)}{\partial x} = U, \qquad (3.2) $$

thus proving the mass conservation (2.6), stated in Section 2.2, and the point (ii) of the theorem.

**Remark 3.1.** It might be interesting to note that the conclusion (3.2) still holds if, in equation (2.3), $b(t)$ would be replaced by any meaningful function $b(x, y, t)$ satisfying $b(1, 1, t) = 1 - F(1, 0, t)$.
### 3.1.2 Equation for $F(1, y, t)$

Setting $A(y, t) \stackrel{\text{def}}{=} F(1, y, t)$, with $A(1, t) = 1$ and $A(0, t) = 1 - b(t)$, equation (2.3) yields

$$ \frac{\partial A(y,t)}{\partial t} + \lambda \left(1 - \frac{1}{y}\right) A(y,t) + \mu(1-y) \frac{\partial F(1, y, t)}{\partial x} = \lambda \left(1 - \frac{1}{y}\right) A(0, t). \quad (3.3) $$

We shall refer to this equation in Section 3.2.

### 3.1.3 Equation for $F(x, 1, t)$

Setting $B(x, t) \stackrel{\text{def}}{=} F(x, 1, t)$, with $B(1, t) = 1$, and $r(t) \stackrel{\text{def}}{=} \frac{\partial F(1, 1, t)}{\partial x}$, equation (2.3) yields immediately

$$ \frac{\partial B(x,t)}{\partial t} + \lambda b(t)(1-x)B(x,t) = \mu(1-x)\frac{\partial B(x,t)}{\partial x}. \quad (3.4) $$

Dividing now equation (3.4) by $1-x$ and letting $x \to 1$, we obtain

$$ \frac{dr(t)}{dt} = \lambda b(t) - \mu r(t), \quad (3.5) $$

noting that $r(t)$ represents the mean number of customers in the first queue of the tandem at time $t$.

**Remark 3.2.** Equation (3.5) seems quite natural, since it can be viewed as a fluid equation for the first queue of the tandem, which is a $M(t)/M/\infty$ queue with arrival rate $\lambda b(t)$ and service rate $\mu$. It can also be obtained directly from equation (2.3), by taking the partial derivative with respect to $y$, then substituting $x = y = 1$ and using (3.2).

The first queue is a $M(t)/M/\infty$ queue with arrival rate $\lambda b(t)$ and service rate $\mu$.

**Assumption H1** The first queue is supposed to be empty at time zero.

Then the departure process of this queue is a non-homogeneous Poisson process, see e.g. [14, p.99], with intensity

$$ \delta(t) \stackrel{\text{def}}{=} \int_{0}^{t} \lambda b(t-x) \mu e^{-\mu x} dx, \quad t > 0, $$

or, after solving (3.5) with the condition $r(0) = 0$ (which corresponds to **H1**),

$$ \delta(t) = \lambda \mu \int_{0}^{t} b(s) e^{-\mu(t-s)} ds = \mu r(t), \quad t > 0. \quad (3.6) $$

For the sake of clarity, let us emphasize that in the rest of Section 3 the variable $r(t)$ will not appear explicitly, and we will focus on $\delta(t)$.

## 3.2 A differential system for $\delta(t)$

We suppose that at time $t = 0$, the first queue $M(t)/M/\infty$ is empty. This state can be reached with probability one by a direct coupling argument, due to the underlying Markovian evolution given by (2.2). Then, from the previous section, it appears that the second queue of the tandem, denoted by $\mathcal{T}_2$, is of $M(t)/M/1$ type, with constant service rate $\lambda$ and Poisson arrival intensity $\delta(t)$, defined by (3.6).

Let $N(t)$ denote the random number of units (namely the number of cars) in $\mathcal{T}_2$ at time $t$, and its associated generating function

$$ N(y,t) \stackrel{\text{def}}{=} \sum_{k \ge 0} \mathbb{P}(N(t) = k) y^k. $$

The original system thus behaves in the same stochastic way as the two following queues in interaction:
• $\mathcal{T}_2$, just defined above;

• $\mathcal{T}_1$, which is of $M(t)/M/\infty$ type, which Poisson arrival rate $\lambda(1 - N(0, t))$ and service rate $\mu$.

We analyze $\mathcal{T}_2$ in more detail, starting from the forward Kolmogorov's equations, which give rise to the functional equation

$$ \frac{\partial N(y,t)}{\partial t} + \left( \lambda \left( 1 - \frac{1}{y} \right) + \delta(t)(1-y) \right) N(y,t) = \lambda \left( 1 - \frac{1}{y} \right) N(0,t), \quad (3.7) $$

where $\delta(t)$ is given by (3.6), keeping in mind the similarity between equations (3.3) and (3.7).

The direct integration of (3.7) leads to

$$ \varphi(y,t)N(y,t) = N(y,0) + \lambda \left(1 - \frac{1}{y}\right) \int_{0}^{t} \varphi(y,s)N(0,s)ds, \quad (3.8) $$

where

$$ \varphi(y,t) = \exp \left\{ \left( 1 - \frac{1}{y} \right) \int_{0}^{t} (\lambda - y\delta(s)) ds \right\}. \quad (3.9) $$

Setting, for notational convenience, $H(t) \stackrel{\text{def}}{=} N(0,t) = 1 - b(t)$, we have from Cauchy's formula $H(t) = \frac{1}{2i\pi} \int_{\mathcal{L}} \frac{N(z,t)}{z} dz$, where $\mathcal{L}$ stands for a simple closed contour around 0.

Then, we get from (3.8) the following Volterra integral equation of the second kind

$$ H(t) = \psi(t) + \lambda \int_{0}^{t} D(s,t)H(s)ds, \quad (3.10) $$

with

$$ \psi(t) = \frac{1}{2i\pi} \int_{\mathcal{L}} \frac{N(z,0)}{z\varphi(z,t)} dz, \quad D(s,t) = \frac{1}{2i\pi} \int_{\mathcal{L}} \frac{\varphi(z,s)}{\varphi(z,t)} \frac{1-1/z}{z} dz. $$

It turns out that $D(s,t)$ can be expressed in terms of the modified Bessel functions

$$ I_n(z) = \frac{1}{2i\pi} \int_{|\omega|=1} \frac{\exp[(z/2)(\omega + 1/\omega)]}{\omega^{n+1}} d\omega = \frac{1}{\pi} \int_0^\pi e^{z \cos(\theta)} \cos(n\theta) d\omega, \quad n \in \mathbb{Z}. $$

Setting $A_t \stackrel{\text{def}}{=} \int_0^t \delta(s)ds$, we have exactly

$$ D(s,t) = e^{-\lambda(t-s)-(A_t-A_s)} \left( I_0\left(2\sqrt{\lambda(t-s)(A_t-A_s)}\right) - \left(\frac{A_t-A_s}{\lambda(t-s)}\right)^{1/2} I_1\left(2\sqrt{\lambda(t-s)(A_t-A_s)}\right) \right). \quad (3.11) $$

Hence, upon combining equations (3.5), (3.6) and (3.10), we have proved the following result.

**Proposition 3.3.** The rate function $\delta(t)$ satisfies the differential system

$$ \left\{ \begin{array}{l} \frac{d\delta(t)}{dt} + \mu\delta(t) = \lambda\mu(1 - H(t)), \quad \delta(0) = 0, \\ H(t) = \psi_{\delta}(t) + \lambda \int_{0}^{t} D(s,t)H(s)ds. \end{array} \right. \quad (3.12) $$

**Remark 3.4.** The first equation of system (3.12) shows that $\frac{d\delta(t)}{dt} + \mu\delta(t) \le \lambda\mu$, with $\delta(0) = 0$, whence

$$ \delta(t) \le \lambda(1 - e^{-\mu t}). $$

In the sequel, and in particular in the proof of the forthcoming lemma, to emphasize the functional dependence on $\delta(t)$, we will add the subscript $\delta$ to the concerned functions, e.g. $H_\delta, \psi_\delta, D_\delta$, etc. This convention also applies to (3.12).
**Lemma 3.5.** The integro-differential system (3.12) has a unique solution $\delta(\cdot)$ with

$$ \delta(0) = 0, \quad \text{and} \quad \delta \stackrel{\text{def}}{=} \lim_{t \to \infty} \delta(t) \qquad (3.13) $$

exists and is a finite positive constant.

*Proof.* According to a method similar to that proposed in [8], we introduce the following iteration scheme:

$$ \left\{ \begin{array}{l} \delta_0(t) = 0, \quad \forall t \ge 0, \\ \displaystyle \frac{d\delta_{n+1}(t)}{dt} + \mu\delta_{n+1}(t) = \lambda\mu(1 - H_{\delta_n}(t)), \\ \delta_{n+1}(0) = 0, \quad \forall n \ge 0, \\ H_{\delta_n}(t) = \psi_{\delta_n}(t) + \lambda \int_0^t D_{\delta_n}(s,t)H_{\delta_n}(s)ds. \end{array} \right. \qquad (3.14) $$

We will show that this scheme is increasing. For each $n \ge 0$, let $Q_n$ denote the $M(t)/M/1/\infty$ queue with FIFO service discipline, arrival rate $\delta_n(t)$ and service intensity $\lambda$. Then, the probability for $Q_n$ to be empty at time $t$ is equal to $H_{\delta_n}(t)$, the unique solution of the Volterra integral equation appearing in system (3.14). Using the simple stochastic monotonicity proved in Appendix A.1, we shall argue by induction on $n$, assuming that all $Q_n$'s have the same initial conditions [this argument is rendered possible by the third equation in (3.14)].

Suppose $\delta_n(t) \ge \delta_{n-1}(t)$, which is in particular true for $n = 1$. Then, by property (ii) of Lemma A.1, we have the inequality $H_{\delta_n}(t) \le H_{\delta_{n-1}}(t)$, which, after rewriting the second equation of (3.14) as

$$ \delta_{n+1}(t) = \lambda\mu \int_0^t e^{\mu(s-t)} (1 - H_{\delta_n}(s)) ds, \qquad (3.15) $$

yields immediately $\delta_{n+1}(t) \ge \delta_n(t)$. Hence, the uniformly bounded sequence $\{\delta_n(t), n \ge 0\}$ (resp. $\{H_{\delta_n}(t), n \ge 0\}$) is, for each fixed $t$, non-decreasing (resp. non-increasing), and the functions

$$ \delta(t) = \lim_{n \to \infty} \delta_n(t) \quad \text{and} \quad H_{\delta}(t) = \lim_{n \to \infty} H_{\delta_n}(t) $$

do exist and satisfy (3.12).

After having shown the existence of a solution to system (3.12), we are left with the problem of uniqueness. In fact, considering the first equation of (3.12) and using standard results on integral equations, this uniqueness is straightforward, since by (iii) of Lemma A.1, $H_{\delta}(\cdot)$ satisfies a Lipschitz condition with respect to $\delta(\cdot)$.

**Lemma 3.6.** Any reachable stationary distribution of the tandem driven by equation (2.3) has a product form $F(x, y) \stackrel{\text{def}}{=} \lim_{t \to \infty} F(x, y, t)$ in the sense that

$$ F(x, y) = e^{p(x-1)} \frac{1-q}{1-qy}, \qquad (3.16) $$

with $p > 0, 0 \le q < 1$ and $\mu p = \lambda q$.

*Proof.* Letting $t \to \infty$ in (2.3), we obtain the following stationary Kolmogorov's equations of a quasi-reversible system, which, after setting $F(x, y) \stackrel{\text{def}}{=} F(x, y, \infty)$, have the form

$$ \left( \delta(1-x) + \lambda \left(1 - \frac{1}{y}\right) \right) F(x, y) = \mu(y-x) \frac{\partial F(x, y)}{\partial x} + \lambda \left(1 - \frac{1}{y}\right) F(x, 0). \qquad (3.17) $$

Then, inserting (3.16) into (3.17), we obtain immediately the necessary relations

$$ \delta = \mu p = \lambda q = \lambda(1 - F(1, 0)), \qquad (3.18) $$

and the proof of Lemma 3.6 is concluded.

■
To proceed further with the proof of Lemma 3.5, we use the mass conservation equation (3.2), which says that the total mean number of units in the system is equal to some given bounded constant $U$, for all $t \ge 0$. Hence, any possible reachable stationary regime corresponds necessarily to a positive recurrent process. Although we do not yet know the behavior of $\delta(t)$ when $t \to \infty$, we see that, for any positive increasing sequence $\{t_k, k \ge 0\}$, with $\lim_{k \to \infty} t_k = \infty$ and $\lim_{k \to \infty} \delta(t_k) = \bar{\delta}$, the equilibrium equation (3.17) is obtained as the limit of (2.3), where $\delta$ is replaced by $\bar{\delta}$. In addition, $\bar{\delta}$ must satisfy (3.18), together with the relation

$$U = \frac{\partial F(1, 1)}{\partial x} + \frac{\partial F(1, 1)}{\partial y} = p + \frac{q}{1-q},$$

or, equivalently,

$$U = \frac{\bar{\delta}}{\mu} + \frac{\bar{\delta}}{\lambda - \bar{\delta}},$$

which yields a second degree equation for $\bar{\delta}$ having the unique admissible root

$$\bar{\delta} = \frac{\mu(U + \rho + 1 - \sqrt{(U + \rho + 1)^2 - 4\rho U})}{2}. \qquad (3.19)$$

Therefore $\bar{\delta}$ does not depend on the choice of the sequence $\{t_k\}$, hence the existence of the limit

$$\delta = \lim_{t \to \infty} \delta(t) = \bar{\delta},$$

given by (3.19). The proof of Lemma 3.5 is concluded. ■

The preceding arguments show that we are entitled to analyze the equations (2.2), (2.3), treating $\lambda b(t)$ as if it were the exogenous function (see (3.5), (3.6))

$$\lambda b(t) = \delta(t) + \frac{1}{\mu} \frac{d\delta(t)}{dt},$$

where $\delta(t)$ satisfies (3.12). After doing this, equations (2.2) —describing the evolution of the process $\xi(t)$ associated with the original tandem queue network shown in Figure 2.1)— become standard forward Kolmogorov’s equations, whose existence and uniqueness of solution follow at once from the general theory (see e.g. [10, 15]), and they satisfy the condition

$$\sum_{j,k} \alpha_{j,k}(t) = 1, \quad \forall t < \infty.$$

To study the steady-state behavior of $\xi(t)$ when $t \to \infty$, we start by claiming, on the basis of equations (2.6) and (3.2), that $\xi(t)$ is necessarily ergodic, which implies in particular $\delta < \lambda$.

Moreover, the invariant measure of $\xi(t)$ is expected to coincide with the invariant measure of the random walk, say $\tilde{\xi}(t)$, obtained just replacing the function $\delta(t)$ by the constant $\delta$. This statement can be established by a simple coupling argument, relying on the continuity of $\delta(t)$. Indeed, for all $\varepsilon > 0$, there exists $T_{\varepsilon}$ such that

$$|\delta(t) - \delta| < \varepsilon, \quad \forall t \ge T_{\varepsilon}. \qquad (3.20)$$

Taking now $\tilde{\xi}(0) = \xi(T_{\varepsilon})$, the result follows directly from inequality (A.5) and Section A.1.

On the other hand, $\tilde{\xi}(t)$ corresponds to a standard Jackson network known to have the product-form stationary distribution (see Lemma 3.6)

$$\pi_{j,k} = e^{-(\delta/\mu)} \frac{(\delta/\mu)^j}{j!} (1 - \beta)^k, \quad \forall j, k \in \mathbb{N},$$

and one checks directly from (3.19) that $\beta = \frac{\delta}{\lambda}$ is given by (2.5).

This concludes the proof of the point (i) of Theorem 2.3.
### 3.3 Speed of convergence of $\delta(t)$

The following proposition provides an estimate of the speed of convergence of $\delta(t)$ toward $\delta$.

**Proposition 3.7.** As $t \to \infty$, $\delta(t) = \delta + \mathcal{O}(e^{-vt})$, with $v = \min(\mu, (\sqrt{\lambda} - \sqrt{\delta})^2)$.

*Proof.* Let $\tilde{Q}_n, n \ge 0$, be the following sequence of ergodic $M/M/1/\infty$ queues with the following parameters for each $\tilde{Q}_n$:

* the service rate is equal to $\lambda$;

* the intensity of the Poisson arrival process has the constant value $\gamma_n < \lambda$.

Let $\tilde{H}_n(t)$ be the probability for $\tilde{Q}_n$ to be empty at time $t$. It is well-known that $\tilde{H}_n(t)$ can be expressed in terms of Bessel integrals (see e.g. [16, p. 23]). Moreover, setting

$$c_n \stackrel{\text{def}}{=} (\sqrt{\lambda} - \sqrt{\gamma_n})^2,$$

we have, as $t \to \infty$, the classical estimate, valid here since $\lambda > \alpha$ (see e.g. [2, p. 107]),

$$\tilde{H}_n(t) = 1 - \frac{\gamma_n}{\lambda} + \tilde{h}_n t^{-3/2} e^{-c_n t} + o(t^{-3/2} e^{-c_n t}), \quad (3.21)$$

where $\tilde{h}_n$ is a uniformly bounded constant depending on the initial conditions.

**Lemma 3.8.** Let $f(\cdot)$ be an integrable function of $t$ on $\mathbb{R}_+$, such that $l = \lim_{t \to \infty} f(t)$ exists. Then

$$\lim_{t \to \infty} \int_0^t \mu e^{\mu(s-t)} f(s) ds = l. \quad (3.22)$$

*Proof.* Just splitting the integral into the sum $\int_0^T + \int_T^t$, where, for a given $\varepsilon$, $T$ is such that

$$|l - f(t)| < \varepsilon, \quad \forall t \ge T,$$

the proof of the lemma is immediate by letting $t \to \infty$, and then $\varepsilon \to 0$.

Let us consider now the following scheme

$$\left\{
\begin{array}{l}
\gamma_0(t) = \gamma \in ]0, \lambda[, \quad \forall t \ge 0, \\
\\
\displaystyle\frac{d\gamma_{n+1}(t)}{dt} + \mu\gamma_{n+1}(t) = \lambda\mu(1 - H_{\gamma_n}(t)), \\
\\
\gamma_{n+1}(0) = 0, \quad \forall n \ge 0, \\
\\
H_{\gamma_n}(t) = \psi_{\gamma_n}(t) + \lambda \int_0^t D_{\gamma_n}(s,t)H_{\gamma_n}(s)ds,
\end{array}
\right.
\quad (3.23)$$

which differs from (3.14) only by its first equation, but this is not a minor point!

The last equation of (3.23) can be seen as emanating from a $M(t)/M/1$ ergodic queue, which we suppose to be empty at time 0. Then, by Takács formula (see [16, Th.1, p. 23]), we can write

$$\tilde{H}_n(t) \ge 1 - \frac{\gamma_n}{\lambda}, \quad \forall \gamma_n \le \lambda.$$

Using now Lemma 3.8 and (3.15) with $\delta_n$ replaced by $\gamma_n$, we get at once $\gamma_1(t) \le \gamma(1 - e^{-\mu t}) \le \gamma$, which implies by induction that $\gamma_n(t)$ is monotone decreasing in $n$, and

$$\lim_{n \to \infty} \gamma_n(t) = \delta(t), \quad \text{with} \quad \delta_n(t) \le \delta(t) \le \gamma_n(t), \quad \forall t, n \ge 0.$$
Suppose, and this is true for $n = 0$,

$$ \gamma_n(t) = \gamma_n + \mathcal{O}(e^{-a_n t}). $$

Using the monotonicity properties already mentioned and the estimate (3.21), we can write

$$ H_{\gamma_n}(t) = 1 - \frac{\gamma_n}{\lambda} + \mathcal{O}(\max(e^{-a_n t}, t^{-3/2}e^{-c_n t})). $$

Then, the integrated form

$$ \gamma_{n+1}(t) = \lambda\mu \int_{0}^{t} e^{\mu(s-t)}(1 - H_{\gamma_n}(s))ds \quad (3.24) $$

yields (some details are omitted),

$$ \gamma_{n+1}(t) = \gamma_n + \mathcal{O}(e^{-a_{n+1}t}) $$

where

$$ a_{n+1} = \min(a_n, v_n), \quad \text{with } v_n = \min(\mu, c_n). $$

Choosing now $\gamma_0(t) = \delta, \forall t \ge 0$, in the first equation of (3.23), we have

$$ \lim_{n \to \infty} a_n = \lim_{n \to \infty} \min(\mu, (\sqrt{\lambda} - \sqrt{\gamma_n})^2) = (\sqrt{\lambda} - \sqrt{\delta})^2, $$

concluding the proof of Proposition 3.7.

Applying now Lemma 3.8 to the integral form obtained from (3.12),

$$ \delta(t) = \lambda\mu \int_{0}^{t} e^{\mu(s-t)}(1 - H_{\delta}(s))ds, $$

we get $\delta = \lambda\beta$. Hence, Proposition 3.7 is clearly equivalent to point (iii) of Theorem 2.3, the proof of which is complete.

# 4 Two other related systems

In this section, we present two systems in the same context with capacity constraints. Both can be analyzed via the same methods, and the main results are stated without detailed proofs.

## 4.1 Model 2: mean-field limit viewed as a single $M(t)/M/1/K$ queue with finite capacity $K$

This model is the simplest one among the car-sharing class (see e.g. [12]), with no parking space reservation. It is particularly meaningful for bike-sharing systems, which do not permit reservation. There are $M_N$ cars moving among $N$ stations of capacity $K$. A car leaves the station at rate $\lambda$ and travels for a random time (exponentially distributed with parameter $\mu$). Then, it returns to a uniformly chosen station, if possible, or else sets off on a new journey.

As $N \to \infty$ with $U = \lim_{N \to \infty} M_N/N$, the usual empirical distribution tends to a deterministic dynamical system $\vec{\alpha}(t)$ described by the following equations, for $0 \le j \le K, t \ge 0$,

$$ \frac{d\alpha_j(t)}{dt} + [\lambda \mathbb{1}_{\{j>0\}} + \mu(U - a(t))\mathbb{1}_{\{j<K\}}] \alpha_j(t) = \mu(U - a(t))\alpha_{j-1}(t) + \lambda\alpha_{j+1}(t)\mathbb{1}_{\{j<K\}}, \quad (4.1) $$

where $a(t) = \sum_{j=0}^{K} j\alpha_j(t)$. The underlying time-inhomogeneous Markov process is a $M(t)/M/1/K$ queue with service rate $\lambda$ and Poisson arrival rate $\mu(U - a(t))$, where $a(t)$ is the mean number of customers in the queue at time $t$.
Setting $Q(z, t) = \sum_{j=0}^{K} \alpha_j(t)z^j$, for an arbitrary complex variable $z$ and $t \ge 0$, we get the nonlinear functional equation

$$ \frac{\partial Q(z, t)}{\partial t} + \left[ \mu(U - a(t))(1 - z) + \lambda \left(1 - \frac{1}{z}\right) \right] Q(z, t) = \lambda \left(1 - \frac{1}{z}\right) Q(0, t) + \mu(U - a(t)) \alpha_K(t) z^K (1 - z), \quad (4.2) $$

with $a(t) = \frac{\partial Q(1, t)}{\partial z}$.

**Theorem 4.1.** For $t \ge 0$, there exists a unique distribution $\tilde{\alpha}(t)$ satisfying the nonlinear system (4.1). In addition, for any $\lambda, \mu > 0$, there is a unique equilibrium point $\lim_{t \to \infty} \tilde{\alpha}(t) = \vec{\pi} > \vec{0}$ given by

$$ \pi_j = \frac{(1 - \beta)\beta^j}{1 - \beta^{K+1}}, \quad 0 \le j \le K, $$

where

$$ \beta = \lim_{t \to \infty} \frac{\mu}{\lambda} (U - a(t)) = \frac{\mu}{\lambda} \left( U - \sum_{j=0}^{K} j \pi_j \right) $$

is the unique solution of the fixed point equation

$$ \beta = \frac{\mu}{\lambda} \left( U - \frac{\beta}{1 - \beta} + \frac{(K + 1)\beta^{K+1}}{1 - \beta^{K+1}} \right). $$

*Proof.* System (4.1) rewrites in the matrix form

$$ \frac{d\tilde{\alpha}(t)}{dt} = \tilde{\alpha}(t)M(a(t)), \quad (4.3) $$

where $\tilde{\alpha}(t)$ is a row vector and $M(a(t))$ a tridiagonal matrix, denoted in this way to emphasize the dependence on $a(t)$. Now we introduce the iterative vector scheme

$$ \left\{ \begin{array}{l} a_0(t) = \tilde{\alpha}_0(t) = 0, \quad t \ge 0, \\ a_n(t) = \sum_{j=0}^{K} j \alpha_{n,j}(t), \\ \tilde{\alpha}_{n+1}(0) = \tilde{\alpha}(0), \\ \frac{d\tilde{\alpha}_{n+1}(t)}{dt} = \tilde{\alpha}_{n+1}(t)M(a_n(t)), \quad n \ge 0, \end{array} \right. \quad (4.4) $$

which can be shown, using Appendix A.1, to be increasing in the sense that $a_n(t) < a_{n+1}(t) < U$.

## 4.2 Model 3: mean-field limit viewed as a system of two queues in tandem with total capacity bounded by $K < \infty$

The model is the same as that described in Sections 1.1 and 1.2, with a finite capacity constraint $K$. The first (resp. second) queue of the tandem is still an $M(t)/M/\infty$ (resp. $M(t)/1$) queue, but the total number of particles remains bounded by $K$. Then, for $j, k \in \mathbb{N}$, $0 \le k + j \le K$, the following mean-field equations hold (see [4]), for all $t \ge 0$,

$$ \frac{d\alpha_{j,k}(t)}{dt} + \left[ \lambda d(t) \mathbb{1}_{\{j+k<K\}} + \mu j + \lambda c(t) \mathbb{1}_{\{k>0\}} \right] \alpha_{j,k}(t) = \lambda d(t) \alpha_{j-1,k}(t) + \mu(j+1) \alpha_{j+1,k-1}(t) + \lambda c(t) \alpha_{j,k+1}(t) \mathbb{1}_{\{j+k<K\}}, \quad (4.5) $$

where

$$ d(t) = 1 - \sum_{j=0}^{K} \alpha_{j,0}(t), \quad \text{and} \quad c(t) = 1 - \sum_{j=0}^{K} \alpha_{j,K-j}(t). \quad (4.6) $$
Setting

$$F(x, y, t) \stackrel{\text{def}}{=} \sum_{0 \le j+k \le K} \alpha_{jk}(t) x^j y^k, \quad \forall t \ge 0,$$

for $x, y$ arbitrary complex variables, the following functional equation holds.

$$\begin{aligned} \frac{\partial F(x, y, t)}{\partial t} + & \left[ \lambda d(t)(1-x) + \lambda c(t) \left(1 - \frac{1}{y}\right) \right] F(x, y, t) \\ = & \mu(y-x) \frac{\partial F(x, y, t)}{\partial x} + \lambda c(t) \left(1 - \frac{1}{y}\right) F(x, 0, t) + \lambda d(t)(1-x) F_K(x, y, t), \end{aligned} \quad (4.7)$$

where $F_K(x, y, t) = \sum_{j=0}^{K} \alpha_{j,K-j}(t)x^j y^{K-j}$, $d(t) = 1 - F(1, 0, t)$ and $c(t) = 1 - F_K(1, 1, t)$.

The analog of Theorem 2.3 can be stated as follows.

**Theorem 4.2.** For each $t \ge 0$ there exists a unique distribution $\bar{\alpha}(t)$ satisfying the nonlinear system (4.1). In addition, for any $\lambda, \mu > 0$, there is a unique equilibrium point $\lim_{t \to \infty} \bar{\alpha}(t) = \bar{\pi} > \bar{0}$ given by

$$\pi_{jk} = \frac{1}{Z} \frac{\rho_R^j}{j!} \rho_V^k, \quad 0 \le j+k \le K, \quad (4.8)$$

where $Z$ is a normalizing constant and $(\rho_R, \rho_V)$ is obtained as the unique solution of the system

$$\left\{ \begin{array}{l} \rho_R = \frac{\lambda}{\mu} \left(1 - \sum_{j=0}^{K} \pi_{j0}\right), \\ U = \sum_{j,k} (j+k)\pi_{jk}. \end{array} \right. \qquad (4.9)$$

As previously announced, we only sketch out a few arguments. To this end, as in Section 3.1, we write for Model 3 the following subsidiary equations, which emanate directly from (4.7).

(1) Letting $G(x, t) \stackrel{\text{def}}{=} F(x, x, t)$, with $G(1, t) = 1$, we have

$$\frac{\partial G(x, t)}{\partial t} + \lambda \left(d(t) - \frac{c(t)}{x}\right) (1-x) G(x, t) = \lambda c(t) \left(1 - \frac{1}{x}\right) F(x, 0, t) + \lambda d(t) (1-x) x^K(t) (1-c(t)). \quad (4.10)$$

Then, dividing (4.10) by $1-x$ and letting $x \to 1$, we get

$$\frac{\partial}{\partial t} \frac{\partial}{\partial x} G(1, t) = 0,$$

or, equivalently,

$$\frac{\partial G(1, t)}{\partial x} = U, \qquad (4.11)$$

which depicts a mass conservation valid for $t \ge 0$.

(2) For $A(y, t) \stackrel{\text{def}}{=} F(1, y, t)$, with $A(1, t) = 1$ and $A(0, t) = 1 - d(t)$, equation (3.3) still holds, just replacing $\lambda$ by $\lambda c(t)$.

$$\frac{\partial A(y, t)}{\partial t} + \lambda c(t) \left(1 - \frac{1}{y}\right) A(y, t) = \mu(1 - y) \frac{\partial F(1, y, t)}{\partial x} + \lambda c(t) \left(1 - \frac{1}{y}\right) A(0, t). \quad (4.12)$$
(3) Analogously, letting $B(x, t) \stackrel{\text{def}}{=} F(x, 1, t)$, with $B(1, t) = 1$,

$$ \frac{\partial B(x, t)}{\partial t} + \lambda d(t)(1-x)B(x, t) = \mu(1-x)\frac{\partial B(x, t)}{\partial x} + \lambda d(t)(1-x)F_K(x, 1, t), $$

which yields the equivalent of (3.5), just replacing $\lambda$ by $\lambda c(t)$, i.e.

$$ \frac{dr(t)}{dt} = \lambda c(t)d(t) - \mu r(t), \qquad (4.13) $$

where $r(t) \stackrel{\text{def}}{=} \frac{\partial F(1, 1, t)}{\partial x}$ stands for the mean number of customers in the first queue of the tandem, at time $t$.

**Remark 4.3.** As equation (3.5), equation (4.13) seems quite natural. Indeed, it can be viewed as a fluid equation for the first queue of the tandem, which is a $M(t)/M/\infty$ queue with arrival rate $\lambda c(t)d(t)$ and service rate $\mu$, due to the additional capacity constraint, which translates into a probability of acceptance $d(t)$.

The following result is the analog of Lemma 3.6.

**Lemma 4.4.** Any reachable stationary distribution of the tandem driven by equation (4.7) has the generating function $F(x, y) \stackrel{\text{def}}{=} \lim_{t \to \infty} F(x, y, t)$,

$$ F(x, y) = \frac{1}{Z} \sum_{j+k \le K} \frac{(px)^j}{j!} (qy)^k, \qquad (4.14) $$

where $Z = \sum_{j+k \le K} \frac{p^j}{j!} q^k$ is the normalizing constant and where $p > 0$ and $0 \le q < 1$ are uniquely determined by $\mu p = \lambda d$ with $d = 1 - Z^{-1} \sum_{j=0}^{K} p^j / j!$, and by the mass conservation equation

$$ U = \frac{1}{Z} \sum_{j+k \le K} (j+k) \frac{p^j q^k}{j!}. \qquad (4.15) $$

**Proof.** Letting $t \to \infty$ in equation (4.7), we obtain the following stationary Kolmogorov's equations of a quasi-reversible system, which, after setting $F(x, y) \stackrel{\text{def}}{=} F(x, y, \infty)$, have the form

$$ \left( \lambda d(1-x) + \lambda c \left( 1 - \frac{1}{y} \right) \right) F(x, y) = \mu(y-x) \frac{\partial F(x, y)}{\partial x} + \lambda c \left( 1 - \frac{1}{y} \right) F(x, 0) + \lambda d(1-x) F_K(x, y), \qquad (4.16) $$

where $c = c(\infty)$ and $d = d(\infty)$.

Then, inserting (4.14) into (4.16), we obtain (after some easy algebra) the necessary relation

$$ \mu p = \lambda d, $$

which, combined with the mass conservation (4.15), yields a system of two equations with two unknowns $p$ and $q$, having a unique solution (see [11, Theorem 3] for details). The proof of Lemma 4.4 is concluded. ■

From the previous lemma, we get directly (4.8), which corresponds to the last point of Theorem 4.2.
# A Queues with time-dependent arrival rates

## A.1 The $M(t)/M/1/\infty$ queue

This section briefly presents some basic features of the operator (birth and death process type) describing the evolution of the $M(t)/M/1/\infty$ queue, with time varying arrival rate $\beta(\cdot) \in C^+[0, \infty]$ and constant service rate $\mu$.

The probabilities $p_n(t) \stackrel{\text{def}}{=} \mathbb{P}(Z(t) = n)$, $n \ge 0$, where $Z(t)$ denotes the number in the queue at time $t$, obey the following set of forward differential equations:

$$ 
\left\{
\begin{aligned}
\frac{dp_0(t)}{dt} &= -\beta(t)p_0(t) + \mu p_1(t), \\
\frac{dp_n(t)}{dt} &= \beta(t)p_{n-1}(t) - (\beta(t) + \mu)p_n(t) + \mu p_{n+1}(t), \quad n \ge 1,
\end{aligned}
\right.
\qquad (A.1) $$

which we rewrite in operator form

$$ \frac{d\mathbf{P}_{\beta}(t)}{dt} = \mathbf{P}_{\beta}(t)\mathbf{K}_{\beta}(t), \qquad (A.2) $$

where $\mathbf{K}_{\beta}(t)$ is a generator (represented by an infinite matrix) and $\mathbf{P}_{\beta}(t)$ is an infinite row vector belonging to the Banach space $\ell_1$ of absolutely summable sequences. It is known either from a probabilistic point of view (e.g. [10, 15]) or by an analytic argument (e.g. [7]), that (A.2) has, for all $t \ge 0$, a unique solution in $\ell_1$. In addition, the generator $\mathbf{K}_{\beta}(t)$ has a continuous spectrum of eigenvalues, located on the negative real line.

Similarly, the distribution function

$$ s_n(t) \stackrel{\text{def}}{=} \mathbb{P}(Z(t) \le n), \quad \forall n \ge 0, $$

satisfies the system

$$ 
\left\{
\begin{aligned}
\frac{ds_0(t)}{dt} &= -(\beta(t) + \mu)s_0(t) + \mu s_1(t), \\
\frac{ds_n(t)}{dt} &= \beta(t)s_{n-1}(t) - (\beta(t) + \mu)s_n(t) + \mu s_{n+1}(t), \quad n \ge 1,
\end{aligned}
\right.
\qquad (A.3) $$

which will be written as

$$ \frac{d\mathbf{S}_{\beta}(t)}{dt} = \mathbf{S}_{\beta}(t)\mathbf{L}_{\beta}(t), \qquad (A.4) $$

where $\mathbf{S}_{\beta}(t)$ denotes the row vector $\mathbf{S}_{\beta}(t) \stackrel{\text{def}}{=} (s_0(t), s_1(t), \dots)$.

**Lemma A.1.** Let $\mathbf{P}_{\beta}(t)$ and $\tilde{\mathbf{P}}_{\beta}(t)$ be the solutions of (A.2) corresponding to respective initial conditions $\mathbf{P}_{\beta}(0)$ and $\tilde{\mathbf{P}}_{\beta}(0)$. The following properties hold:

(i) If $\mathbf{P}_{\beta}(0) \ge \tilde{\mathbf{P}}_{\beta}(0) \ge 0$, then $\mathbf{P}_{\beta}(t) \ge \tilde{\mathbf{P}}_{\beta}(t) \ge 0$, $\forall t \ge 0$.

(ii) Let $\beta(\cdot), \gamma(\cdot) \in C^+[0, \infty]$, such that

$$ \beta(t) \le \gamma(t), \quad \forall t \ge 0, \quad \text{and} \quad \mathbf{S}_{\beta}(0) \ge \mathbf{S}_{\gamma}(0). $$

Then (stochastic dominance)

$$ \mathbf{S}_{\beta}(t) \ge \mathbf{S}_{\gamma}(t), \quad \forall t \ge 0. $$

In particular, with the notation of equation (3.12), $H_{\beta}(t) \ge H_{\gamma}(t)$, $\forall t \ge 0$.

(iii) For any $\beta(\cdot), \gamma(\cdot) \in C^+[0, \infty]$, with $\mathbf{P}_{\beta}(0) = \mathbf{P}_{\gamma}(0)$, we have the Lipschitz condition

$$ |\mathbf{P}_{\beta}(t) - \mathbf{P}_{\gamma}(t)| \le K \|\beta(\cdot) - \gamma(\cdot)\|, \qquad (A.5) $$

where $K$ is an absolute constant independent of $t$ and $|\cdot|$ denotes, for an arbitrary vector $\mathbf{X}$, the usual $\ell_1$-norm $|\mathbf{X}| = \sum_i |x_i|$.
(iv) The $M(t)/M/1/\infty$ queue with arrival rate $\delta(t)$ and the $M/M/1/\infty$ queue with arrival rate $\delta = \lim_{t\to\infty} \delta(t)$ have the same stationary regime.

*Proof.* As for point (i), it is not difficult to see that $\mathbf{K}_\beta(t)$ and $\mathbf{L}_\beta(t)$ are positive operators. For instance, making in A.3 the change of functions

$$s_n(t) = w_n(t) \exp(-\mu t - \int_0^t \beta(s)ds), \quad n \ge 0,$$

leads to the system

$$\begin{cases} \frac{dw_0(t)}{dt} = \mu w_1(t), \\ \frac{dw_n(t)}{dt} = \beta(t)w_{n-1}(t) + \mu w_{n+1}(t), \quad n \ge 1, \end{cases}$$

which has the form

$$\frac{d\mathbf{W}_\beta(t)}{dt} = \mathbf{W}_\beta(t)\tilde{\mathbf{L}}_\beta(t),$$

where $\tilde{\mathbf{L}}_\beta(t)$ has only positive coefficients. A similar argument can be used for the positivity of $\mathbf{K}_\beta(t)$.

The stochastic dominance in (ii) follows now from (i). Indeed, setting

$$\mathbf{R}(t) = \mathbf{S}_\beta(t) - \mathbf{S}_\gamma(t),$$

the row vector $\mathbf{R}(t)$ satisfies the non homogeneous differential equation

$$\frac{d\mathbf{R}(t)}{dt} = \mathbf{R}(t)\mathbf{L}_{\gamma}(t) + (\gamma(t) - \beta(t))\mathbf{D}(t), \quad (A.6)$$

where $\mathbf{D}(t) = (d_0(t), d_1(t), \dots)$, with $d_n(t) = [s_{n+1}(t) - s_n(t)]_\beta$. By (i), the vector $\mathbf{D}(t)$ has non-negative components and the operator $\mathbf{L}_\gamma(t)$ is positive, whence it follows that the solutions of (A.6) are also non negative.

To prove (iii), we shall use differential calculus in Banach spaces. In this framework, most of the classical results for the real line or the complex plane apply without substantial modification.

For any $\beta(\cdot) \in C^+[0, \infty]$, with $0 \le \|\beta\| \le B$, take an arbitrary perturbation function $\Delta(\cdot)$, with $\beta(t) + \Delta(t) \in C^+[0, \infty]$. When it exists, the partial derivative with respect to $\beta(\cdot)$ of a differentiable mapping

$$g: C^+[0, \infty] \times [0, \infty] \to \ell_1$$

is a functional (see [7]) written $\frac{\partial g(t)}{\partial \beta}$. With this notation, one sees easily that $\mathbf{Q}_\beta \stackrel{\text{def}}{=} \frac{\partial \mathbf{P}_\beta(t)}{\partial \beta}$, where $P_\beta$ satisfies (A.2), must be a solution of the following non-homogeneous linear differential equation

$$\frac{d\mathbf{Q}_\beta(t)}{dt} = \mathbf{Q}_\beta(t)\mathbf{K}_\beta(t) + \mathbf{P}_\beta(t)\mathbf{M}, \quad \mathbf{Q}(0) = 0, \quad (A.7)$$

where $\mathbf{M}$ is a constant infinite matrix given by

$$\begin{pmatrix} -1 & 1 & 0 & 0 & \dots \\ 0 & -1 & 1 & 0 & \dots \\ 0 & 0 & -1 & 1 & \dots \\ 0 & 0 & 0 & -1 & \dots \\ \vdots & \vdots & \vdots & \vdots & \dots \end{pmatrix}.$$

Then the solution of (A.7) writes in the form

$$\mathbf{Q}_{\beta}(t) = \int_{0}^{t} \mathbf{P}_{\beta}(s) \mathbf{M} \Phi(t, s) ds,$$
where $\Phi(t, s)$ is the so-called fundamental solution (see [7]) of the homogeneous equation of type (A.2). The preceding argument yields directly the rough estimate

$$ \sup_{\beta} \|Q_{\beta}\| = \sup_{\beta} \sup_{t \ge 0} \|Q_{\beta}(t)\| \le K, \qquad (A.8) $$

where $K$ is a bounded constant. Since $Q_{\beta}$ is the derivative with respect to $\beta(\cdot)$ of the function $P_{\beta}$, defined on the Banach space $C^{+}[0, \infty]$, (A.8)) gives the Lipschitz condition (A.5).

The last property (iv) in the lemma can be viewed as an immediate consequence of the stochastic ordering contained in (i) and (ii), and details will be omitted. The proof of the Lemma is concluded.

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

## References

[1] ALLMEIER, S., AND GAST, N. Mean field and refined mean field approximations for heterogeneous systems: It works! *Proceedings of the ACM on Measurement and Analysis of Computing Systems* 6, 1 (2022), 1–43.

[2] ASMUSSEN, S. *Applied probability and queues, Second Edition*, vol. 51 of *Applications of mathematics*. Springer, 2003.

[3] BASKETT, F., CHANDY, K. M., MUNTZ, R. R., AND PALACIOS, F. G. Open, closed, and mixed networks of queues with different classes of customers. *Journal of the ACM (JACM)* 22, 2 (1975), 248–260.

[4] BOURDAIS, C., FRICKER, C., AND MOHAMED, H. A mean field analysis of a stochastic model for reservation in car-sharing systems. *ACM SIGMETRICS Performance Evaluation Review* 48, 2 (2020), 18–20.

[5] BUDHIRAJA, A., DUPUIS, P., FISCHER, M., AND RAMANAN, K. Limits of relative entropies associated with weakly interacting particle systems. *Electronic Journal of Probability* 20, 80 (2015), 1–22.

[6] BUDHIRAJA, A., DUPUIS, P., FISCHER, M., AND RAMANAN, K. Local stability of Kolmogorov forward equations for finite state nonlinear Markov processes. *Electronic Journal of Probability* 20, 81 (2015), 1–30.

[7] CARTAN, H. *Calcul différentiel*. Hermann, Paris, 1967.

[8] DELCOIGNE, F., AND FAYOLLE, G. Thermodynamical limit and propagation of chaos in polling systems. *Markov Processes and Related Fields* 5, 1 (1999), 89–124.

[9] ETHIER, S. N., AND KURTZ, T. G. *Markov Processes. Characterization and Convergence*. John Wiley & Sons, 2009.

[10] FELLER, W. *An Introduction to Probability Theory and its Applications*, second ed., vol. Volume II. Wiley, 1971.

[11] FRICKER, C., AND MOHAMED, H. Mean-field analysis of stochastic networks with reservation. *Journal of Applied Probability* (2025), 1–27.

[12] FRICKER, C., AND TIBI, D. Equivalence of ensembles for large vehicle-sharing models. *Annals of applied probability: an official journal of the Institute of Mathematical Statistics* 27, 2 (2017), 883–916.

[13] GAST, N., AND VAN HOUDT, B. A refined mean field approximation. *Proceedings of the ACM on Measurement and Analysis of Computing Systems* 1, 2 (2017), 1–28.

[14] GNEDENKO, B. W., AND KÖNIG, D. *Handbuch der Bedienungstheorie II : Formeln und andere Ergebnisse*. De Gruyter, Berlin, Boston, 1984.

[15] REUTER, G. E. H., AND LEDERMANN, W. On the differential equations for the transition probabilities of Markov processes with enumerably many states. *Proc. Cambridge Philos. Soc.* 49 (1953), 247–262.

[16] TAKÁCS, L. *Introduction to the theory of queues*. University Texts in the Mathematical Sciences. Oxford University Press, New York, 1962.