# Efficient reinforcement learning with Fleming-Viot particle systems: application to stochastic networks with rarely observed rewards

Daniel Mastropietro, Urtzi Ayesta, Matthieu Jonckheere, Szymon Majewski

► To cite this version:

Daniel Mastropietro, Urtzi Ayesta, Matthieu Jonckheere, Szymon Majewski. Efficient reinforcement learning with Fleming-Viot particle systems: application to stochastic networks with rarely observed rewards. 2023. [hal-04129885](https://hal.science/hal-04129885)

HAL Id: hal-04129885

https://hal.science/hal-04129885

Preprint submitted on 21 Jul 2023

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# Fast-exploring reinforcement learning with applications to stochastic networks

Daniel Mastropietro, Urtzi Ayesta, Matthieu Jonckheere, Szymon Majewski

► To cite this version:

Daniel Mastropietro, Urtzi Ayesta, Matthieu Jonckheere, Szymon Majewski. Fast-exploring reinforcement learning with applications to stochastic networks. Queueing Systems, inPress, 109 (3), pp.23. 10.1007/s11134-025-09950-5 . hal-04129885v2

HAL Id: hal-04129885

https://hal.science/hal-04129885v2

Submitted on 10 Sep 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
# Efficient reinforcement learning with Fleming-Viot particle systems: application to stochastic networks with rarely observed rewards

Daniel Mastropietro¹, Urtzi Ayesta², Matthieu Jonckheere³, and Szymon Majewski⁴

¹CNRS-IRIT, Université de Toulouse INP, 31071 Toulouse, France

²CNRS-IRIT, 31071 Toulouse, France / Ikerbasque - Basque Foundation for Science, 48009 Bilbao, Spain / UPV/EHU, University of the Basque Country, 20018 Donostia, Spain

³CNRS-LAAS, CNRS-IRIT, 31071 Toulouse, France

⁴Ecole Polytechnique, Paris, France

July 14, 2023

## Abstract

We consider reinforcement learning control problems under the expected reward criterion in which non-zero rewards are both sparse and rare, that is, they occur in very few states and have a very small stationary probability under all policies. In this context, usual discovery techniques including importance sampling are inapplicable because no policy exists that increases the visit frequency of the rare states. Using renewal theory and Fleming-Viot particle systems, we propose a novel approach that exploits prior knowledge on the sparse structure of the reward landscape to boost exploration of the rare non-zero rewards and achieve an accurate estimation of their stationary probability. We also demonstrate how to combine the methodology with policy gradient learning to construct the FVRL algorithm that efficiently solves control problems under these scenarios. We provide theoretical guarantees of the convergence of both the stationary probability estimator and the policy gradient learner, and illustrate the method on two optimisation problems to maximize the expected reward: a simple $M/M/1/K$ queue system where the blocking threshold $K$ is optimised, and a two-job-class loss network where a threshold-type rejection policy is optimised. Our results show that FVRL learns the optimum thresholds much more efficiently than vanilla Monte-Carlo reinforcement learning.

**Keywords:** semi-Markov decision process, parameterised policy, policy gradient, queues, loss network

## 1 Introduction

Reinforcement learning methods, by being able to take advantage of large amounts of computational resources, have been very successful at solving very complex problems with large state dimensions and sparse rewards, obtaining super-human performance particularly in games [23]. Many of these successes have been obtained in an episodic setting in which, even under sparse rewards conditions, there is certainty that the episode will eventually finish, at which moment a reward will be observed.

In several application domains, in particular in networking or robotics, the environment is not episodic (i.e. there is no notion of progression as in games) and non-zero rewards are received only in a handful of states which are visited very rarely. We refer to these as environments with sparse and rare rewards. For example in networking, a fundamental problem is how to dimension the system in order to optimise the performance, bearing in mind that the blocking probability (i.e. the probability that the system cannot accept a new data packet, call, or computation task) can be extremely small. As pointed out in [11], efficiently managing the exploration task when non-zero

rewards are very rarely observed remains a challenge, and this provides the main motivation for the present work. Our starting point is the structural knowledge on the underlying Markov Decision Process (MDP) with which the system is modelled, that in many cases can be leveraged to drastically improve exploration.

In this paper, we focus on model-free approaches (i.e. without previous learning of a specific model), and take as performance criterion the expected reward (or expected cost, when more natural for the problem at hand). We assume that we have access to a simulator or emulator of the system, or in its absence, to a large amount of data to construct experience replay of the system. Our approach relies on the identification of sets $\mathcal{A}$ of states with zero reward, which can be obtained through previous knowledge on the underlying MDP, or from information gained from exploration already performed.

We first show that the expected reward of the original problem can be expressed in terms of modified trajectories that are absorbed in $\mathcal{A}$. We then show that the expected reward of trajectories outside $\mathcal{A}$ can be efficiently estimated via the so-called Fleming-Viot particle system (FV). We finally introduce FVRL, a reinforcement learning algorithm that uses Fleming-
# Fast-exploring reinforcement learning with applications to stochastic networks*

Daniel Mastropietro<sup>1,2,*</sup>, Urtzi Ayesta<sup>2,3,4</sup>, Matthieu Jonckheere<sup>2</sup>, and Szymon Majewski<sup>5</sup>

<sup>1</sup>Université de Toulouse INP, 31071 Toulouse, France

<sup>2</sup>CNRS-IRIT, 31071 Toulouse, France

<sup>3</sup>Ikerbasque - Basque Foundation for Science, 48009 Bilbao, Spain

<sup>4</sup>UPV/EHU, University of the Basque Country, 20018 Donostia, Spain

<sup>5</sup>Ecole Polytechnique, Paris, France

September 9, 2025

## Abstract

We introduce FVRL (Fleming-Viot reinforcement learning), a reinforcement learning algorithm for optimisation problems where a long-term objective is largely influenced by states that are very rarely observed under all policies. In this context, usual discovery techniques including importance sampling are inapplicable because no alternative policy exists that increases the observed frequency of the rare states. We instead propose a novel approach that uses Fleming-Viot particle systems, a family of stochastic processes evolving simultaneously under the same law, that exploits prior knowledge of the environment to boost exploration of the rare states. A renewal theory argument allows us to consistently estimate the stationary probability of the rare states from excursions that have considerably lower sample complexity than usual Monte-Carlo explorations. We demonstrate how to combine this estimator with policy gradient learning to construct the FVRL algorithm, which is suited to efficiently solve problems where the optimisation function is expressed as a long-run expectation, such as the long-run expected reward. We show that the FVRL algorithm converges to a local optimizer of the parameterised objective function, and illustrate the method on two optimisation problems that aim at minimizing the long-run expected cost under admission control policies of threshold type: a simple $M/M/1$ queue system and a two-job-class loss network. Our experimental results show that, under the same sample complexity, FVRL outperforms a vanilla Monte-Carlo reinforcement learning method by converging to the optimum thresholds considerably faster.

<footnote>*Corresponding author: daniel.mastropietro@gmail.com</footnote>

<footnote>Research partially supported by the French National Research Agency (ANR) under the programs "Investments for the Future" with reference code CIMI ANR-11-LABX-0040 and ANR-22-CE25-0013-02 (ANR EPLER), and by project STIC-AMSUD LAGOON.</footnote>
Viot to solve control problems with sparse and rare rewards and gradients.

To illustrate the main ideas, we shall consider an $M/M/1/K$ queue, where the objective is to minimize a convex expected cost on blocking events. This is a fundamental model in queueing theory and an interesting archetype of dynamics where the rewards (here seen as costs) can only be observed during unlikely excursions for system loads smaller than 1. More concretely, the stationary probability of being in state $K$ decreases exponentially with $K$ at rate $\rho^K$, where $\rho$ is the load of the system, assumed smaller than 1. The rewards are thus sparse, since blocking only occurs when the system is full, and rare since this happens with very low probability for a wide range of values of $K$ and $\rho$ that are common in practice.

Using this example, the intuition behind our proposal relies on the fact that the probability of visiting state $K$ conditioned on not having visited state 0, which decays as $\sqrt{K}\rho^{K/2}$ when $K$ increases [10], is considerably larger than the stationary probability of visiting $K$. For example, for $K = 40$ and $\rho = 0.7$, the conditioned probability is about 8 thousand times larger than the unconditioned one, and about 6 million times larger when $\rho = 0.5$. Hence, by shifting the estimation toward a conditioned dynamics, which we propose to achieve using Fleming-Viot particle systems, we expect a lower sample complexity, both for estimation and for control. Our numerical results show that the Fleming-Viot approach overcomes the difficulty of vanilla Monte-Carlo to properly estimate key quantities for reinforcement learning algorithms, such as value functions and gradients.

In summary, we propose a methodology that can be used to boost reinforcement learning in environments with the following characteristics: (i) they present sparse non-zero rewards and sparse gradients of parameterised policies occurring very rarely at unknown states; (ii) thanks to prior knowledge, it is possible to define a set of states presenting zero rewards. The learning goal considered here is to maximize the expected reward under stationary regime operation but this could be easily generalised to contexts with discounts. In environments where the occurrence of non-zero rewards has an extremely low probability, this methodology enables learning at admissible learning times, whereas vanilla Monte-Carlo methods fail.

In our main contribution, we develop FVRL, a method that combines FV and RL in problems with sparse and rare rewards that are ubiquitous in the control of stochastic networks. In our main theoretical contributions for FVRL, we establish convergence of the estimator as the number of particles tends to infinity, and convergence to the optimal policy. Our numerical results on a simple queueing model and a multi-class loss network illustrate the potential benefits of the approach, and we will investigate in further research its applicability to a wider variety of examples, including classical RL environments.

We would also like to note that, to the best of our knowledge, there is no adequate solution to this problem in the state of the art. Even though we do not focus here on deep learning tools, our proposal could definitely be used in combination with them. This is left for future work. Initially, one might think that the problem of estimating rare events could be tackled using Importance Sampling (IS), an area in which there exists a large literature. However, it is important to observe that the problem

at hand impedes the application of IS methods. Firstly, we here look at scenarios where **there is no policy** allowing to explore rare states. In other words, states with non-zero rewards are very rarely explored under all policies, so IS is not an option, as it is designed to change the exploration policy. Secondly, it is also not possible to invoke the original IS principle when looking purely at the evaluation task for a fixed policy, as this requires a change of measure for Monte-Carlo, which is not possible in our case as the transition rates of the Markov process are unknown.

On the other hand, the problem of sparse rewards has been of central interest for a long time and the most intuitive solution to sparse (but not necessarily rare) reward problems is reward shaping. Mataric formulated the idea back in 1994 [17]. However, these methods have a few drawbacks: expertise is needed to shape the rewards, and only very few and quite arbitrary policies might be reached as a consequence of shaping. Our proposal fits more into the idea of curiosity-driven methods as explained for instance in [18]. The idea is basically to encourage the exploration of unvisited states in the environment. While existing mechanisms are based on including a bonus term in the loss function to favour exploration ([18, 6]), we define a new (more radical) curiosity mechanism adapted to sparse and rare reward environments that forces exploration outside a set of states that have already been explored or are known to be uninformative.

The rest of the paper is organized as follows. Section 2 describes the mathematical setting of the problem, Section 3 the general methodology, in Section 4 we show its applicability in the cases of an $M/M/1/K$ queue system and of a loss network serving two classes of jobs.

## 2 Problem description

We consider a continuous-time MDP $(S, A, q, \mathcal{R})$ with a finite state space $S$, action space $A$, jump rates $q$, and rewards $\mathcal{R}$, under the expected reward criterion. Throughout the paper we assume that for each policy $\pi$ the continuous-time Markov process $X_t^\pi$ obtained by following the policy $\pi$ is irreducible. We denote by $p^\pi$ the stationary distribution of $X_t^\pi$, and by $\mathbb{E}^\pi(\eta)$ the expectation with respect to $p^\pi$ of a function of interest, $\eta: S \to \mathbb{R}$, such as the occupation measure or the reward.

We will be interested in computing $\mathbb{E}^\pi(\eta)$ under the assumption that the function $\eta$ is zero outside of a small set of states $C \subset S$. For example, if the rewards are assumed to be sparse, the reward function $r: S \times A \to \mathbb{R}$ is zero outside a small set $C \times A$, from which we can define the function $\eta$ as $\eta(x) = \sum_a r(x, a)\pi(a|x)$, whose expected value can be reduced to a calculation over the states in set $C$ as $\mathbb{E}^\pi(\eta) = \sum_{x \in C} p^\pi(x)\eta(x)$.

The objective of computing $\mathbb{E}^\pi(\eta)$ will be achieved as follows: we will first choose a set $A \subset S$, such that $C \cap A = \emptyset$. Then, we will use Fleming-Viot particle systems to compute the expectation truncated to $A^c$, $\sum_{x \in A^c} p^\pi(x)\eta(x)$, which is equal to $\mathbb{E}^\pi(\eta)$ given that $\eta$ is zero in $A$.
# 1 Introduction

Reinforcement learning methods (RL) [25] are very versatile in solving stochastic optimisation problems thanks to their little model assumptions. In the model-free paradigm, they only assume that the stochastic system can be modeled as a Markov decision process (MDP), but no knowledge of the transition probabilities or rates are required. From interactions between an agent and the system (environment), RL methods are able to learn near optimal policies by trial and error. The agent tries different possible actions on many different environment states and uses the signals received in the form of rewards or costs to incrementally adjust its decisions towards a better policy. Although this approach is appealing by the potential to tackle very complex optimisation problems whose modeling can be challenging, it makes RL algorithms a computational intensive task. Depending on the environment size and characteristics, the learning process can be prohibitively slow due to lack of sufficient coverage of all possible states and actions, some of which can be highly informative. This makes the design of efficient exploration schemes one of the key components of RL algorithms, particularly in environments where the reward landscape is sparse, which usually leads to situations where informative rewards are rarely observed.

Environments with rarely observed states (under all policies) are the object of study of the present work, in particular in the context of problems where the optimisation function is largely influenced by visiting those rare states. More specifically, we look into problems where the optimisation function is based on a long-run expectation criterion, such as the long-run expected reward. For example, a network system will have its long-run quality of service significantly degraded if the cost of interruption is large even when such interruptions are rare.

We propose to tackle this problem by leveraging prior or learned knowledge about the MDP in order to *force* exploration away from frequently visited states that are uninformative for the optimisation task at hand, while keeping the current policy in place. We integrate this prior knowledge into the parameterisation of a population evolution model known as Fleming-Viot in order to efficiently estimate the stationary probability distribution of the rare high-reward states, which is a fundamental piece of the long-run expected criterion being optimised.

To fix ideas, consider the MDP defined by a simple $M/M/1$ queue system and an admission control policy that is solely a function of the queue size, $x$. When the system load $\rho$ is less than 1, we know from standard queueing theory that most of the excursion of the Markov chain defined by the queue size, takes place at lower values of $x^1$. If the policy rejects an incoming job when the queue is at a predefined length $K$ at which point a negative reward is accrued, observing a non-zero reward may be very rare. For instance, by PASTA (Poisson arrivals see time averages), the rejection probability is of the order of $10^{-7}$ when $K = 40$ and $\rho = 0.7$. Given a rejection cost function, reinforcement learning could be used to obtain the optimal $K$ threshold that minimizes the long-run expected cost estimated from the job rejection events observed during a Monte-Carlo simulation of the MDP. This, however, may be a prohibitively time consuming task because Monte-Carlo would require a very long simulation when the rejection probability is very small. To tackle this problem, we instead propose a Fleming-Viot simulation (explained in the *Contributions* paragraph below) that forces the Markov chain excursions away from the high probability and non-informative lower occupation states, to obtain an estimation of the long-run expected cost objective with considerably lower sample complexity.

**Related work** The problem of sparse and rare rewards has been of central interest for a long time, as it is directly linked to the core topic of efficient exploration of RL environments. We cite a few approaches that have been proposed to tackle this problem: importance sampling, reward shaping, reward-free exploration, and curiosity-driven methods.

**Importance sampling** is an offline method that proposes exploring the environment using an alternative exploration policy specially designed to increase the occurrence probability of the rare events under the original policy. An appropriate adjustment factor is then applied to obtain a correct estimate of expected values that depend on the observation frequency of these events [27]. An important constraint of the alternative exploration policy is that it must have the same state coverage as the original policy, i.e. all states reachable by the original policy must also

<footnote n="1">¹The stationary probability at $x$ is $\rho^x(1 - \rho)$, an exponentially decreasing function of $x$ when $\rho < 1$.</footnote>
# 3 Methodology

In this section we present a new method to estimate $\mathbb{E}^{\pi}(\eta)$ tailored to Markov decision processes with sparse and rare rewards. We then discuss how to use this algorithm to improve the estimation of gradients in the context of the policy gradient methodology to solve optimal control problems with sparse and rare rewards. For simplicity of exposition, we assume that we have direct access to the function $\eta$. We also assume throughout this section that the policy $\pi$ and the chosen set $A \subset S$, whose intersection with $C$ is empty, are fixed.

## 3.1 Definitions

The following definitions will be instrumental in our discussion. We denote by $\vec{\partial}A$ the entrance boundary of $A$, i.e. the set of states $x \in A$ for which there exists at least a state $y \in A^c$ with positive jump rate to $x$, i.e. $q(y, x) > 0$. The entrance boundary of $A^c$, $\vec{\partial}A^c$, is defined analogously. We define the first time of entry into $A$ as

$$ \mathcal{T}_{A,0} \doteq \inf\{t > 0 : X_t^\pi \in A \text{ and } X_{t-}^\pi \notin A\}, $$

and we denote the entrance state distribution into $A$ under stationarity as:

$$ p_{\vec{\partial}A}^{\pi}(x) \doteq \mathbb{P}(X_{T_{A,0}}^{\pi} = x | X_0^{\pi} \sim p^{\pi}), \forall x \in \vec{\partial}A. $$

Using Figure 1 as a visual aid, we further define the first time of entry into $A^c$ following $\mathcal{T}_{A,0}$ as $T_{A^c} \doteq \inf\{t > \mathcal{T}_{A,0} : X_t^\pi \in A^c\}$, the first time of entry into $A$ following $T_{A^c}$ as $T_A \doteq \inf\{t > T_{A^c} : X_t^\pi \in A\}$, and their difference as $T_K \doteq T_A - T_{A^c}$, also referred to as the killing time. The stopping time $T_A$ will be regarded in the sequel as a cycle return time to $A$.

Finally, for any measurable subset $B$, we define $\mathbb{P}_{\vec{\partial}A}(B) \doteq \mathbb{P}(B|X_0^\pi \sim p_{\vec{\partial}A}^\pi)$. For the complement set $A^c$, we define the entrance state distribution into $A^c$ under stationarity as:

$$ p_{\vec{\partial}A^c}^{\pi}(x) \doteq \mathbb{P}_{\vec{\partial}A}(X_{T_{A^c}}^{\pi} = x), \forall x \in \vec{\partial}A^c $$

The two state probability distributions defined above, $p_{\vec{\partial}A}^{\pi}(x)$ and $p_{\vec{\partial}A^c}^{\pi}(x)$, will be thoroughly used in the development of our proposed methodology to condition the start state of the Markov decision process $X_t^\pi$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Stopping times defined in Section 3.1.</figcaption></figure>

## 3.2 Estimation with Fleming-Viot particle systems

We now introduce the Fleming-Viot method as a way to tackle the problem of sparse and rare rewards and estimate $\mathbb{E}^{\pi}(\eta)$ efficiently. Our approach penalizes trajectories that enter $A$ in order to boost exploration of the subset $C$ of the state space that is relevant for the estimation of the expectation, which is in $A^c$. The penalisation consists in immediately replacing trajectories that enter $A$ by trajectories outside $A$. To this end, we use the dynamics of a particle system known as Fleming-Viot (FV) [7] which has been used in the literature to simulate quasi-stationary distributions [7, 2].

More specifically, the Fleming-Viot $N$-particle system with driving process $X_t^\pi$ and absorption set $A$ is a continuous-time Markov process $(\xi_t^\nu)_{t \ge 0}$ defined on the state space $(A^c)^N$ as follows: given a probability distribution $\nu$ on $A^c$, an $N$-dimensional vector $\xi_t^\nu(k)_{k=1,...,N}$ defines the initial state of the particles in the system and is obtained as i.i.d. samples from $\nu$. Each particle $\xi_t^\nu(k)$ then evolves independently according to the dynamics of $X_t^\pi$, but whenever it hits a state in $A$, it immediately jumps to the position of one of the other particles chosen uniformly at random. This mechanism allows us to only explore trajectories outside $A$, which is where the informative rewards are located.

In order to exploit the Fleming-Viot particle system for the estimation of $\mathbb{E}^{\pi}(\eta)$, we leverage the renewal theory characterization of the stationary probability of a state in terms of return cycles [1, chapter 6, Theorem 1.2], as follows: if we use the entrance to $A$ as the event defining the beginning and end of a cycle, the stationary expectation of an arbitrary function $\eta$ can be written in terms of the cycle time $T_A$ as

$$ \mathbb{E}^{\pi}(\eta) = \frac{\mathbb{E}_{\vec{\partial}A} \left( \int_{0}^{T_A} \eta(X_t^{\pi}) dt \right)}{\mathbb{E}_{\vec{\partial}A} T_A}. \quad (1) $$

Note that, if we used Monte-Carlo to estimate this expectation by simulating the Markov process $X_t^\pi$ starting at $\vec{\partial}A$ and observing cycle times $T_A$, it might take a very long time before observing a non-zero contribution from $\eta$, as this function is assumed to be sparse with non-zero values rarely observed. The following proposition is key to define the Fleming-Viot estimation method.

**Proposition 1.** Given a set $A \subset S$ and a function $\eta: S \to \mathbb{R}$ that is zero on $A$, the following holds:

$$ \mathbb{E}^{\pi}(\eta) = \frac{\int_{0}^{\infty} \mathbb{E}_{\vec{\partial}A^c} \left( \eta(X_t^{\pi}) \mathbf{1}_{T_K > t} \right) dt}{\mathbb{E}_{\vec{\partial}A} T_A}, \quad (2) $$

which, for estimation purposes, can be further simplified as:

$$ \mathbb{E}^{\pi}(\eta) = \int_{0}^{\infty} f_{\eta}(t)g(t)dt, \quad (3) $$

where

$$ f_{\eta}(t) \doteq \sum_{x \in A^c} \eta(x) \phi_t^{\vec{\partial}A^c}(x); g(t) \doteq \frac{\mathbb{P}_{\vec{\partial}A^c}(T_K > t)}{\mathbb{E}_{\vec{\partial}A} T_A}, $$

and $\phi_t^{\vec{\partial}A^c}(x) \doteq \mathbb{P}_{\vec{\partial}A^c}(X_t^\pi = x | T_K > t)$ is the probability that the process $X_t^\pi$, started at a state in $\vec{\partial}A^c$ chosen with probability $p_{\vec{\partial}A^c}^\pi$, is in $x$ provided it has not been absorbed.
be reachable by the alternative policy. However, this is not always possible. For instance, in the simple $M/M/1$ queue example mentioned above, there is no admission control policy that is able to increase the visit frequency of the rejection event while maintaining the state coverage of the original policy: if the rejection threshold $K$ is large and the alternative proposed policy decreases its value in order to increase the frequency of the rejection event, all the states between the original and the new smaller threshold will never be visited.

**Reward shaping** was initially proposed by Mataric back in 1994 [18], as a method where an expert proposes or shapes the reward landscape of an environment in order to obtain policies that are able to visit certain states of interest. We believe this type of approach has two main drawbacks: domain knowledge is needed to shape the rewards, and only a few (due to handcraft shaping) and uncontrolled policies may be reached as a result of shaping.

**Reward-free exploration methods** propose exploration mechanisms that are based exclusively on the probability of state visit, regardless of the reward landscape –and thus of policies–, as their goal is to collect enough data that can later be used to optimise policies under *any* reward context [16]. These methods propose reducing the sample complexity by identifying a set of states worth visiting based on a minimum visit frequency criterion, assuming the states with very small visit probability have negligible impact on reward optimisation. Although, like ours, these methods aim at reducing sample complexity, they ignore the states at the tail of the visit frequency distribution, which nonetheless could bring relevant information for the optimisation process if their rewards happen to be sufficiently large. The method does not provide a sample complexity reduction mechanism for these states, which is the focus of our work.

**Curiosity-driven methods** encourage the exploration of less visited states by including a bonus term in the objective function [20, 6]. Our method could be considered a type of curiosity-driven method, with the important difference that, as we will see below, the exploration of the rare states is not just encouraged, but *forced*, by prohibiting excursions from visiting a subset of more frequently visited and uninformative states.

**Reinforcement learning for stochastic networks:** Several works have considered RL in the context of stochastic network models. Dai and Gluzman [15] propose a deep reinforcement learning algorithm for queueing network control. This study adapts the proximal policy optimisation (PPO) algorithm to Markov Decision Processes (MDPs) with infinite state spaces, unbounded costs, and long-run average cost objectives, which are common in queueing networks. Our primary goal here is actually to address a challenge explicitly mentioned by Dai and Gluzman as future work: developing more sample-efficient simulation methods by incorporating learned knowledge about the environment. Additionally, we test our method on multi-class stochastic networks, a benchmark also suggested by Dai and Gluzman for evaluating RL methods.

On the other hand, the work in [13] establishes convergence of Natural Policy Gradient algorithm, when initialized with the MaxWeight policy, in infinite-state queueing systems under the average-reward criterion. While deep reinforcement learning algorithms often fail to produce stable policies in unbounded state spaces, [21] introduces a Lyapunov-based cost-shaping technique and state transformations that improve stability and performance in queueing networks. We also refer to [4] which proposes a machine learning approach to the optimal control of multi-class fluid queueing networks. In contrast to the aforementioned papers, in our case, the optimal policy heavily depends on a state that is rarely visited, a problem that to the best of our knowledge has not yet been addressed in the queueing and stochastic network literature.

**Contributions** Our main contribution is the proposal of an exploration strategy based on a population evolution model known as Fleming-Viot particle system (FV), designed with the goal of increasing the discovery rate of rare states compared to the natural exploration of the environment by the MDP. The ultimate goal is the design of an efficient policy learning mechanism for optimisation problems with long-run objectives that may be affected by large rewards observed on those rare states.

Since learning a policy requires the ability to measure the function being optimised, we start by proposing FVEE (Fleming-Viot expectation estimator), an estimation scheme of the stationary probability distribution, based on renewal theory, that will allow us to estimate long-run expectations. In the FV context, individuals are called
Proof. Since at the beginning of a cycle the process starts at a state in $\partial_A$, in order to reach a state in the set of interest $C \in A^c$ (the set of states with non-zero values of $\eta$), the process needs to go through $\tilde{\partial}_A A^c$.

Thus, using the definition of the sojourn times $T_{A^c}$ and $T_A$ and recalling that $T_A > T_{A^c}$, expression (1) can be written as

$$ \mathbb{E}^{\pi}(\eta) = \frac{\mathbb{E}_{\tilde{\partial}_A}(\int_{0}^{T_{A^c}} \eta(X_t^{\pi})dt + \int_{T_{A^c}}^{\infty} \eta(X_t^{\pi})\mathbf{1}_{T_A>t}dt)}{\mathbb{E}_{\tilde{\partial}_A}T_A} $$

Since by assumption $\eta$ is zero in $A$, the first integral is zero. After changing the integration variable to $u = t - T_{A^c}$ in the second integral and using the notation of the killing time $T_K = T_A - T_{A^c}$, we get

$$ \mathbb{E}^{\pi}(\eta) = \frac{\mathbb{E}_{\tilde{\partial}_A}(\int_{0}^{\infty} \eta(X_{u+T_{A^c}}^{\pi})\mathbf{1}_{T_K>u}du)}{\mathbb{E}_{\tilde{\partial}_A}T_A} $$

At $u = 0$, the Markov process is at $X_t^{\pi}$ which, as stated in Section 3.1, is distributed according to $p_{\tilde{\partial}_A}^{\pi}$ when $X_0^{\pi} \sim p_{\tilde{\partial}_A}^{\pi}$, as is the case above. This allows us to redefine the time origin of the Markov process at $u = 0$ and replace $\mathbb{E}_{\tilde{\partial}_A}$ with $\mathbb{E}_{\tilde{\partial}_A}^{c}$ to obtain

$$ \mathbb{E}^{\pi}(\eta) = \frac{\mathbb{E}_{\tilde{\partial}_A}^{c}(\int_{0}^{\infty} \eta(X_{u}^{\pi})\mathbf{1}_{T_K>u}du)}{\mathbb{E}_{\tilde{\partial}_A}T_A} $$

Expression (2) follows from replacing $u \to t$ and inter-changing the order of the integral and the expectation, and expression (3) is directly obtained from (2) by conditioning on the event $\mathbf{1}_{T_K>t}$. $\square$

An estimator $\hat{\mathbb{E}}^{\pi}(\eta)$ of (3) is constructed by estimating each function inside the integral, as follows: the denominator of $g(t)$, $\mathbb{E}_{\tilde{\partial}_A}(T_A)$, is estimated using regular Monte-Carlo from observations of the stopping time $T_A$ coming from the simulation of $X_t^{\pi}$, $f_{\eta}(t)$ and the numerator of $g(t)$, $\mathbb{P}_{\tilde{\partial}_A}^{c}(T_K > t)$, are estimated from the simulation of the FV N-particle system driven by $X_t^{\pi}$ with absorption set $A$. The estimation details are given in Appendix A.

## 3.3 Bound on the estimation error

It has been proved that for finite state spaces [8, 9], uniform in time propagation of chaos holds for Fleming-Viot particle systems. We let $m(\cdot, \xi) : A^c \to [0, 1]$ denote the empirical distribution of the $N$ particles with positions described by vector $\xi$, defined as the empirical mean $m(x, \xi) = \frac{1}{N} \sum_{i=1}^{N} \mathbf{1}_{\xi(i)=x}, \forall x \in A^c$. It then follows that if $\nu$ is a probability measure on $S$, then, under assumption of proper initialization, [8, Theorem 1.4] shows the following bound on the speed of convergence w.r.t. the number of particles $N$ of the empirical mean $m(\cdot, \xi_t^\nu)$ towards $\phi_t^\nu$:

$$ \sup_{x \in S} \sup_{t \ge 0} \mathbb{E} \left| [m(x, \xi_t^\nu)] - \phi_t^\nu(x) \right| \le \frac{C_{FV}}{\sqrt{N}}, \quad (4) $$

where $C_{FV}$ is a positive constant depending on the characteristics of the driving process.

Using this result, we can show the error bound stated in the following theorem for the estimator $\hat{\mathbb{E}}^{\pi}(\eta)$ of (3) described at

the end of Section 3.2. The theorem is valid in an idealized case where the simulation used to estimate the denominator of $g(t)$, $\mathbb{E}_{\tilde{\partial}_A}(T_A)$, is started according to $p_{\tilde{\partial}_A}^{\pi}$, and the Fleming-Viot simulation used to estimate $f_{\eta}(t)$ and the numerator of $g(t)$, $\mathbb{P}_{\tilde{\partial}_A}^{c}(T_K > t)$, is started from i.i.d. samples of $p_{\tilde{\partial}_A}^{\pi}$.

**Theorem 2.** Assume that we start the simulation of $X_t^{\pi}$ for the estimator of $\mathbb{E}_{\tilde{\partial}_A}(T_A)$ according to the distribution $p_{\tilde{\partial}_A}^{\pi}$, that $M$ return cycles to $A$ under stationarity are observed during that simulation, and that we compute the estimator of $f_{\eta}$ in (3) and of $\mathbb{P}_{\tilde{\partial}_A}^{c}(T_K > t)$ using the FV particle system started at the positions of $N$ i.i.d samples from $p_{\tilde{\partial}_A}^{\pi}$. Let $\eta$ be a bounded state function such that $\eta(x) = 0$ for $x \in A$. Then, there exists a constant $C > 0$ such that

$$ \mathbb{E} |\hat{\mathbb{E}}^{\pi}(\eta) - \mathbb{E}^{\pi}(\eta)| \le C \left( \frac{1}{\sqrt{M}} + \frac{1}{\sqrt{N}} \right). $$

The proof is given in Appendix B.

**Remark 1.** Note that the estimator $\hat{\mathbb{E}}^{\pi}(\eta)$ is in general a biased estimator of $\mathbb{E}^{\pi}(\eta)$.

**Remark 2.** Theorem 2 ensures that the estimation error is of the same order as Monte-Carlo, which gives minimal guarantees for FV. However, the estimation error should not be our only focus to evaluate the difference between FV and Monte-Carlo, especially when the ultimate goal is the convergence of a reinforcement learning algorithm. Indeed, in the control problem, a noisy but still informative signal might be very useful compared to no signal at all. Our main idea when replacing MC by FV is to trade the observation of a very rare event in the original problem by the observation of a more common event for FV particles. Although a fully rigorous analysis of the probability to observe a non-zero reward is out of the scope of this paper, we can give rough estimations using the results of [8, 12]. The dynamics of a tagged particle of the FV process converges when $N \to \infty$ to a one-dimensional Markov process having as stationary measure the quasi-stationary distribution $\nu_{QS}$ of the original process (see [12]), where

$$ \nu_{QS}(B) = \lim_{t \to \infty} \phi_t^{\tilde{\partial}_A^c}(B), $$

for any set of states $B$. If the state space is finite, this one dimensional process in turn converges in distribution exponentially fast to its stationary distribution, $\nu_{QS}$. Hence, the probability of finding a non-zero signal (by visiting states in set $C$) for the FV process in a finite time interval is of the order of $\nu_{QS}(C)$, which can be significantly larger than the original $p^{\pi}(C)$ since the dynamics of FV particles within the restricted state space $A^c$ are the same as those of the original process $X_t^{\pi}$. Thus, FV allows a shift of the estimation target, from $\mathbb{E}^{\pi}(\eta)$ to $\mathbb{E}^{\pi}_{\nu_{QS}}(\eta)$, where the latter can be significantly larger than the former.

## 3.4 FVRL: Policy gradient learning with Fleming-Viot particle systems

In this subsection we show how the Fleming-Viot estimation introduced in Section 3.2 can be combined with the policy gradient theorem to solve optimal control problems in environments with sparse and rare rewards under the expected reward criterion. For example, there are many MDPs (see for instance [22, 21, 14, 3, 15]) whose optimal policy is known to be of
particles which evolve according to the MDP dynamics, and are said to be absorbed when reaching the condition of visiting a predefined subset A of the environment's state space, called the absorption set. When a particle is absorbed, it is instantly regenerated with the characteristic given by the position (state) of a particle chosen uniformly at random among all the other particles that are still outside A. Thus, all particles are forced to stay outside A, a set that intuitively should contain frequently visited states, as we are interested in boosting the exploration of states visited rarely.

Armed with this estimator of the stationary probability of rare states, we then introduce FVRL (Fleming-Viot reinforcement learning), an RL control algorithm that integrates FVEE to estimate the gradient of a policy gradient approach, that may also be non-zero in a handful of states rarely observed. This allows the policy gradient algorithm to learn a policy with a smaller sample complexity than Monte-Carlo.

We provide a proof of consistency of FVEE and guarantees of convergence of FVRL to a local optimizer of the parameterised objective function, valid for finite state and action spaces. We showcase their application to optimisation problems on stochastic networks, namely a one-dimensional problem on the $M/M/1$ queueing system, and a two-dimensional problem on a two-job-class loss network. The results are compared with a vanilla Monte-Carlo benchmark, used both for expectation estimation and policy learning.

The rest of the paper is organized as follows. Section 2 describes the mathematical setting of the problem, Section 3 presents the general methodology, divided into the expectation estimation problem tackled by FVEE in Section 3.2, and the control problem tackled by FVRL in Section 3.3. Section 4 presents the results of applying FVEE and FVRL on an $M/M/1$ queue system (Sections 4.1.1 and 4.1.2) and on a loss network serving two classes of jobs (Sections 4.2.1 and 4.2.2). Section 5 concludes, summarizing the assumptions and limitations, and ideas for future work.

## 2 Problem description

We consider a continuous-time MDP $(S, U, q, \mathcal{R}, \Pi)$ with a finite state space $S$, finite action space $U$, jump rates $q \in \mathbb{R}^+$, bounded rewards $\mathcal{R}$, and family of policies $\Pi$, under the average reward criterion. Throughout the paper, as customary in the literature (e.g. [22, Chapter 11]), we assume the continuous-time Markov process $X_t^\pi$ obtained by following policy $\pi \in \Pi$ is irreducible, making the MDP ergodic. We denote by $p^\pi$ the stationary probability distribution of $X_t^\pi$, and by $\mathbb{E}^\pi(\eta)$ the long-run expectation (i.e. with respect to $p^\pi$, or equivalently $\mathbb{E}^\pi(\eta) = \lim_{t \to \infty} \mathbb{E}^\pi(\eta(X_t^\pi)))$ of a bounded function of interest, $\eta: S \to \mathbb{R}$, such as the occupation measure or the expected one-step reward.

We will be interested in estimating $\mathbb{E}^\pi(\eta)$ under the assumption that the function $\eta$ is zero everywhere except in a set of states $C \subset S$. For example, if the rewards are assumed to be sparse, the reward function $r: S \times U \to \mathbb{R}$ is allowed to be non-zero only in $C \times U$, from which we can define the expected one-step reward function $\eta$ as $\eta(x) = \sum_{a \in U} r(x, a)\pi(a|x)$. Thus, the expected value of $\eta$ can be reduced to a calculation over the states in set $C$ as

$$ \mathbb{E}^\pi(\eta) = \sum_{x \in C} p^\pi(x)\eta(x), \qquad (1) $$

where we have used the ergodicity of the MDP under $\pi$ to write $\mathbb{E}^\pi(\eta)$ as a $p^\pi$-expectation. The objective of estimating $\mathbb{E}^\pi(\eta)$ will be achieved as follows: we will first choose an absorption set $A \subset S$, such that $C \cap A = \emptyset$. Then, we will use Fleming-Viot particle systems to construct an estimator $\hat{p}_{FV}^\pi(x)$ of the stationary probability in $C$, and finally compute $\hat{\mathbb{E}}^\pi(\eta) = \sum_{x \in C} \hat{p}_{FV}^\pi(x)\eta(x)$ when $\eta(x)$ is known, or $\hat{\mathbb{E}}^\pi(\eta) = \sum_{x \in C} \hat{p}_{FV}^\pi(x)\hat{\eta}(x)$, for some appropriate estimator $\hat{\eta}(x)$, when only the support of $\eta(x)$ is known.

In this work we will consider three different $\eta$ functions for illustration of the methodology: (i) the occupation measure, $\eta(x) = \mathbf{1}_{X_t^\pi=x}$ (where $\mathbf{1}_B$ is the indicator function of set $B$), which will allow us to estimate the stationary probability $p^\pi(x)$ of states of interest; (ii) the expected one-step reward function, $\eta(x) = \bar{r}(x) = \sum_{a \in U} r(x, a)\pi(a|x)$,
threshold-type (because of underlying monotonicity properties); when the rewards structure in these MDPs is sparse and rare, learning the optimal thresholds with policy gradient algorithms becomes very slow because the gradient is zero except in very few and rarely observed states. We provide here a general description of the FVRL method we propose to speed up learning in those types of scenarios, which we will later apply to two network optimisation problems described in Section 4: an M/M/1/K queue and a loss network system serving two classes of jobs.

Let us consider the case in which an agent interacts with an environment (available in practice either through a simulator, an emulator, or through experience replay of historical data), with the aim of maximizing the expected reward. As customary in the literature, we let $S_n, R_n$ and $A_n$ denote the state, reward and action at the n-th discrete time step, respectively. We denote by $\pi_\theta$ the policy function parameterised by $\theta$. We recall that we assume that a non-zero reward is accrued only in states in C.

It follows from classical MDP theory [20] that the optimal policy for the expected reward criterion is also optimal for the state value (bias) function defined as $v^{\pi_\theta}(x) \doteq E^{\pi_\theta} [\sum_{n=0}^{\infty} (R_n - \bar{R}^{\pi_\theta}) | S_0 = x]$, where $\bar{R}^{\pi_\theta} \doteq \lim_{T \to \infty} \frac{1}{T} E^{\pi_\theta} \sum_{n=0}^{T} R_n$ is the expected reward. We also define the corresponding state-action value function by

$$Q^{\pi_\theta}(x, a) = E^{\pi_\theta} \left[ \sum_{n=0}^{\infty} (R_n - \bar{R}^{\pi_\theta}) \mid S_0 = x, A_0 = a \right]. \quad (5)$$

It thus follows that $v^{\pi_\theta}(x) = \sum_a \pi_\theta(a|x) Q^{\pi_\theta}(x, a)$. We seek to find the policy $\pi_\theta$ that maximizes the average state value, $v^{\pi_\theta} = \sum_x p^{\pi_\theta}(x) v^{\pi_\theta}(x) = \sum_x p^{\pi_\theta}(x) \sum_a \pi_\theta(a|x) Q^{\pi_\theta}(x, a)^1$.

We propose to use a gradient-based algorithm to learn the parameter $\theta$ that maximizes $v^{\pi_\theta}$. Denoting by $X$ the random variable associated to the stationary distribution under $\pi$, it follows from the Policy Gradient Theorem [24] that the gradient of $v^{\pi_\theta}$ can be written as

$$\begin{aligned} \nabla_{\theta} v^{\pi_{\theta}} &= E^{\pi_{\theta}}[Q^{\pi_{\theta}}(X, a) \nabla_{\theta} \pi_{\theta}(a|X)] \\ &= \sum_{x \in S} p^{\pi_{\theta}}(x) \sum_{a} Q^{\pi_{\theta}}(x, a) \nabla_{\theta} \pi_{\theta}(a|x) \end{aligned} \quad (6)$$

We note that the gradient to learn $\theta$ is a linear combination of the policy gradients weighted by the stationary probability $p^{\pi_\theta}$. Considering the assumed sparse and rare structure of the rewards in the environment to control, it is critical to observe the rare states sufficiently often in order to obtain a non-zero estimation of the stationary probability at those states. Otherwise, any non-zero policy gradient $\nabla_{\theta}\pi_{\theta}(a|x)$ coming from these rarely observed states (which are the most informative for learning $\theta$ in the context of the threshold-type policies being discussed) will not contribute to the gradient of the average state value, thus making learning difficult.

In this context, the FV particle system can be leveraged to obtain informative gradients of the average state value through

an accurate estimation of the stationary probability $p^{\pi_\theta}$, in particular at the rarely observed states. The combination of FV and RL leads to the FVRL algorithm, which is used to estimate the gradient of the average state value and learn the optimum $\theta$.

The FVRL algorithm is derived from the following two-step estimation of the gradient of the average state value:

1. For each $x$ where the policy gradient is non-zero, estimate the function of interest $\eta$ introduced in Section 2, here defined as $\eta(x) = \sum_a Q^{\pi_\theta}(x, a) \nabla_{\theta} \pi_\theta(a|x)$.

2. Given an estimate $\hat{\eta}$ of $\eta$, use the FV procedure to estimate its expectation: $\sum_{x \in S} p^{\pi_\theta}(x) \hat{\eta}(x)$.

This means that, besides estimating $p^{\pi_\theta}$ with Fleming-Viot, we need to estimate $\eta$ for every $x$ where the policy gradient is not zero for at least one possible action $a$. Thus, an estimate of $\eta$ is constructed from estimates of $Q(x, a)$ for all actions $a$ whose policy gradient at $x$ is not zero. Although these $Q(x, a)$ values can be estimated by simulating separate Markov processes $X_t^\pi$ independently –each starting at one of the different $(x, a)$ contributing to $\eta$– here we propose a method based on coupled trajectories that leverages the sum-to-one property of the policy over all actions at state $x$. More precisely, the estimation method is as follows: if we let $A_x$ be the set of possible actions when the system is at state $x$, we run $|A_x|$ copies of the Markov process on the extended state-action space, each starting at $(x, a_i)_{i=1, \dots, |A_x|}$. When two such chains meet at the same state-action, they continue evolving together forever. Thus, after all chains meet, the contribution to $\sum_a \hat{Q}^{\pi_\theta}(x, a) \nabla_{\theta} \pi_\theta(a|x)$ is zero, because at that point all values contributing to $\hat{Q}^{\pi_\theta}(x, a)$ are the same, and the partial derivatives of the policy sum up to zero for any given $\theta$. Therefore, chains need to be run only until they all meet. Also, because of the sum-to-one property of the policy, the contribution to $\hat{\eta}$ from the term $\bar{R}^{\pi_\theta}$ in (5) is zero and therefore its value does not need to be computed, only the observed rewards $R_n$ need to be recorded.

The following proposition states that the FVRL algorithm converges with probability 1 to the optimum parameter $\theta^*$.

**Proposition 3.** Given a continuously differentiable parameterisation of the average state value $v^{\pi_\theta}$, the policy gradient FVRL algorithm converges with probability 1 to the optimum parameter $\theta^*$.

The proof is a consequence of [13, Section 5.2, Theorem 2.1] and [4], that consider the convergence of stochastic approximation algorithms with bias, as is the case with the FV estimator used in the estimation of the gradient of the average state value.

# 4 Application to stochastic networks with blocking

In this section we apply the methodologies outlined in Sections 3.2 and 3.4 to two different systems:

<footnote n="1">¹With the same arguments presented in [20], it can easily be proven that the optimal policy for the expected reward criterion is optimal, not only for the state value (bias) function, but also for the average state value (i.e. its expectation over all states).</footnote>
which will allow us to estimate the long-run expected reward $g \doteq \lim_{t \to \infty} \mathbb{E}^{\pi}(\bar{r}(X_t^{\pi}))$ and (iii) its gradient, $\nabla_{\theta}g(\pi_{\theta})$, presented in Section 3.3, which will allow us to update a policy $\pi_{\theta}$ parameterised by $\theta \in \mathbb{R}^d$ in control problems.

In the case of control problems, in Section 4 we will consider queueing problems where the optimal admission control policy is of threshold type. To solve the problem, we will choose a parameterised policy whose gradient is non-zero only near the threshold, making the gradient $\nabla_{\theta}g(\pi_{\theta})$ sparse. We will then use Fleming-Viot particle systems to estimate the sparse gradient as described in Section 3.3.

**Remark 1.** As briefly mentioned in the Introduction, the absorption set $\mathcal{A}$ can be defined following an initial exploration of the environment, by adding states that are visited frequently enough, e.g. above an appropriate relative frequency threshold. For ease of illustration, $\mathcal{A}$ will be assumed known throughout the paper.

**Remark 2.** We make the assumption of zero reward in the absorption set $\mathcal{A}$ to concentrate the attention on our proposed methods. This assumption however can be relaxed, considering that traditional expectation estimation methods (e.g. Monte-Carlo) can be used to estimate the contribution to $\mathbb{E}^{\pi}(\eta)$ from frequently visited states.

# 3 Methodology

In this section we present a new method to estimate $\mathbb{E}^{\pi}(\eta)$ tailored to Markov decision processes with sparse and rare rewards. We then discuss how to use this algorithm to improve the estimation of gradients in the context of the policy gradient methodology to approximately solve optimal control problems with sparse and rare gradients. For simplicity of exposition and unless specified otherwise, we assume that we have direct access to the function $\eta$. We also assume throughout this section that policy $\pi$ and the chosen set $\mathcal{A} \subset \mathcal{S}$, whose intersection with $\mathcal{C}$ is empty, are fixed.

## 3.1 Definitions

The following definitions will be instrumental in our discussion. We denote by $\tilde{\partial}\mathcal{A}$ the entrance boundary of $\mathcal{A}$, i.e. the set of states $x \in \mathcal{A}$ for which there exists at least a state $y \in \mathcal{A}^c$ with positive jump rate to $x$, i.e. $q(y, x) > 0$. The entrance boundary of $\mathcal{A}^c$, $\tilde{\partial}\mathcal{A}^c$, is defined analogously. We define the first time of entry into $\mathcal{A}$ as

$$ \mathcal{T}_{\mathcal{A},0} \doteq \inf\{t > 0 : X_t^{\pi} \in \mathcal{A} \text{ and } X_{t-}^{\pi} \notin \mathcal{A}\}, $$

and we denote the entrance state distribution into $\mathcal{A}$ under stationarity as:

$$ p_{\tilde{\partial}\mathcal{A}}^{\pi}(x) \doteq \mathbb{P}(X_{\mathcal{T}_{\mathcal{A},0}}^{\pi} = x | X_0^{\pi} \sim p^{\pi}), \forall x \in \tilde{\partial}\mathcal{A}. $$

Using Figure 1 as a visual aid, we further define the first time of entry into $\mathcal{A}^c$ following $\mathcal{T}_{\mathcal{A},0}$ as $T_{\mathcal{A}^c} \doteq \inf\{t > \mathcal{T}_{\mathcal{A},0} : X_t^{\pi} \in \mathcal{A}^c\}$, the first time of entry into $\mathcal{A}$ following $T_{\mathcal{A}^c}$ as $T_{\mathcal{A}} \doteq \inf\{t > \mathcal{T}_{\mathcal{A}^c} : X_t^{\pi} \in \mathcal{A}\}$, and their difference as $T_K \doteq T_{\mathcal{A}} - \mathcal{T}_{\mathcal{A}^c}$, also referred to as the killing time. The stopping time $T_{\mathcal{A}}$ will be regarded in the sequel as a cycle return time to $\mathcal{A}$.

Finally, for any measurable subset $B$, we define $\mathbb{P}_{\tilde{\partial}\mathcal{A}}(B) \doteq \mathbb{P}(B | X_0^{\pi} \sim p_{\tilde{\partial}\mathcal{A}}^{\pi})$. For the complement set $\mathcal{A}^c$, we define the entrance state distribution into $\mathcal{A}^c$ under stationarity as:

$$ p_{\tilde{\partial}\mathcal{A}^c}^{\pi}(x) \doteq \mathbb{P}_{\tilde{\partial}\mathcal{A}}(X_{\mathcal{T}_{\mathcal{A}^c}}^{\pi} = x), \forall x \in \tilde{\partial}\mathcal{A}^c. $$

The two state probability distributions defined above, $p_{\tilde{\partial}\mathcal{A}}^{\pi}(x)$ and $p_{\tilde{\partial}\mathcal{A}^c}^{\pi}(x)$, will be thoroughly used in the development of our proposed methodology to condition the start state of the Markov decision process $X_t^{\pi}$.
1. an M/M/1/K queue, which is a unidimensional state system that will help us illustrate the basic concepts;

2. a loss network system with $R$ servers receiving two classes of jobs at exponentially distributed inter-arrival times served with exponentially distributed service times, which we will use to illustrate the application to a multi-dimensional state system.

In both cases, we assume there exists a small subset $C$ of the state space where costs may be generated due to the rejection of an incoming job (blocking), and that there is no reward for job acceptance. For each system we will consider two goals:

a) the accurate estimation of the expected cost, $\mathbb{E}^\pi(\eta)$, defined in Section 2;

b) the efficient learning of the blocking sizes that minimize the expected cost.

To simplify the exposition and to be able to do theoretical computations against which results are compared, we fix the parameters of the underlying Markov process. Nevertheless, the method is designed to be applied in practice to systems with possibly unknown parameters, the only condition being that the system can be simulated or emulated.

In all the simulations run to present the results of the method, either to estimate the expected cost or to learn the optimum blocking sizes, a burn-in period of 10 system transitions is considered before assuming that the stationary regime has been reached for the computation of estimators. In addition, the stationary probability estimate $p^\pi(x)$ –appearing in the calculation of the expected cost and of the policy gradient– is computed only when the number of observed return cycles to the absorption set $A$, used in the estimation of $\mathbb{E}_{\bar{\delta}A}(T_A)$, is at least 5 once the burn-in period has been completed$^2$. If this condition is not satisfied the estimate of the stationary probability $p^\pi(x)$ is set to 0, except for the $M/M/1/K$ case where it is left undefined, reducing the number of successful estimations.

The algorithms for the estimation and for the optimisation problems in both the $M/M/1/K$ queue system and the loss network are presented in Algorithms 1 and 2 in Appendix D, respectively.

In each application problem, the method's performance is compared with vanilla Monte-Carlo (MC) which is used as benchmark. In order to assure a fair comparison, the MC estimation is based on the same number of events as the FV estimation, and the MC learning is based on the same number of events per learning step as the average number of events per learning step observed in FVRL. The details of the estimation and learning processes by MC are described in each subsection below.

## 4.1 M/M/1/K system

We start by defining the Markov process $X_t^\pi$ that, under policy $\pi$, describes the dynamics of the $M/M/1/K$ system, namely the queue size that measures the number $x$ of jobs waiting to be served in the buffer at a given time including the job being served. Jobs arrive following a Poisson process and their

service times are assumed exponentially and independently distributed. We consider that the policy $\pi$ belongs to the family of deterministic accept/reject policies that reject an incoming job at just one state, when $X_t^\pi = K$, in which case a unitary rejection cost is accrued.

The process is thus a continuous-time discrete-state stochastic process living in $\{0, 1, \dots, K\}$ with upward jump rate $\lambda$ from any state $x$ to $x+1$ (except at $x=K$), and downward jump rate $\mu$ from $x$ to $x-1$ (except at $x=0$). We denote by $\rho = \lambda/\mu$ the load of the system which is typically smaller than 1 in real applications. Thus, the system can be represented by an MDP ($S = \{x: 0 \le x \le K\}$, $A = \{0, 1\}$, $q = \{\lambda, \mu\}$, $\mathcal{R} = r(x, a) = \mathbb{1}_{\{x=K, a=0\}}$) where actions 0 and 1 represent "reject" and "accept" an incoming job, respectively.

All presented results correspond to an $M/M/1/K$ system that serves, at rate $\mu = 1$, jobs arriving at rate $\lambda = 0.7$.

### 4.1.1 Estimation of the expected cost

We apply the methodology outlined in Section 3.2 to efficiently estimate the expected cost when blocking is a rare event.

Since the rejection cost is 1, estimating the expected cost is tantamount to estimating the rejection probability, which can be quite small depending on the system parameters. For example with $\rho = 0.5$ and $K = 20$, the rejection probability becomes of the order of $10^{-6}$, a result obtained from the PASTA property of Poisson arrivals and the calculation of the stationary probability where rejection occurs, which is $\mathbb{P}(X_t^\pi = K) = \frac{(1-\rho)\rho^K}{1-\rho^{K+1}}$ whenever $\rho < 1$.

In order to estimate the blocking probability using Fleming-Viot, we define the function $\eta$ introduced in Section 2 as $\eta(x) = \mathbb{1}_{\{x=K\}}$, which is sparse and non-zero in the single-state set $C = \{K\}$. As a consequence, the expectation $\mathbb{E}^\pi(\eta)$ in (3) becomes the blocking probability, which can be written as

$$ \mathbb{E}^\pi(\eta) = p^\pi(K) = \frac{\int_0^\infty \phi_t^J(K) \mathbb{P}_J(T_K > t) dt}{\mathbb{E}_{J-1}T_A}, \quad (7) $$

where we have used that any set $A = \{0, 1, \dots, J-1\}$ with $J \le K$ is a valid absorption set, making $\bar{\delta}A$ and $\bar{\delta}A^c$ two single-state entrance boundary sets equal to $\{J-1\}$ and $\{J\}$, respectively, which makes it possible to simplify (3) into (7).

According to Theorem 2, the estimator of $\mathbb{E}^\pi(\eta)$ converges to its true value as both the number $N$ of particles of the FV system, and the number $M$ of return cycles to $A$ increase. For simulation purposes, and since the number $M$ of return cycles to $A$ is random, we will refer to the number of arrival events $T$ (which directly impacts $M$) when defining the hyperparameters of experiments. The three quantities in (7) contributing to the blocking probability are estimated following the steps described in Algorithm 1 in Appendix D, which implements the estimation methodology described in Appendix A.

**Remark 3.** There is a trade-off between choosing a small or a large value of $J$, the state that defines the size of the absorption set $A$: for smaller $J$, the return times to $A$ will be smaller,

<footnote>$^2$The return cycles to $A$ are measured from the first time the system enters $A$ after the burn-in period has been overcome.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1: Sketch showing the stopping times defined in Section 3.1.</figcaption></figure>

## 3.2 FVEE: Expectation estimation with Fleming-Viot particle systems

We now introduce the Fleming-Viot method as a way to tackle the problem of sparse and rare rewards and obtain a consistent estimator of $E^\pi(\eta)$, which we call FVEE (Fleming-Viot expectation estimator). Our approach proposes a hard penalisation of trajectories that enter $\mathcal{A}$ in order to boost exploration of the subset $\mathcal{C}$ of the state space that is relevant for the estimation of the expectation, which is outside $\mathcal{A}$. The hard penalisation consists in immediately replacing trajectories that enter $\mathcal{A}$ by trajectories outside $\mathcal{A}$. To this end, we leverage the dynamics of the Fleming-Viot particle system [7] which has been used in the literature to simulate quasi-stationary distributions [7, 3].

More specifically, the Fleming-Viot $N$-particle system with driving process $X_t^\pi$ and absorption set $\mathcal{A}$ is a continuous-time Markov process $(\xi_t^\nu)_{t \ge 0}$ defined on the state space $(\mathcal{A}^c)^N$ as follows: given a probability distribution $\nu$ on $\mathcal{A}^c$, an $N$-dimensional vector $\xi_0^\nu(\bar{k})_{k=1,...,N}$ defines the initial state of the particles in the system and is obtained as i.i.d. samples from $\nu$. Each particle $\xi_t^\nu(k)$ then evolves independently according to the dynamics of $X_t^\pi$, but whenever it hits a state in $\mathcal{A}$, it immediately jumps to the position of one of the other particles chosen uniformly at random. This mechanism allows us to only explore trajectories outside $\mathcal{A}$, which is where the informative rewards are located.

In order to exploit the Fleming-Viot particle system for the estimation of $E^\pi(\eta)$, we leverage the renewal theory characterization of the stationary probability of a state in terms of return cycles [2, Chapter 6, Theorem 1.2], as follows: if we use the entrance to $\mathcal{A}$ as the event defining the beginning and end of a cycle, the stationary expectation of an arbitrary function $\eta$ can be written in terms of the cycle time $T_A$ as

$$E^\pi(\eta) = \frac{\mathbb{E}_{\vec{\partial}\mathcal{A}}\left(\int_0^{T_A} \eta(X_t^\pi) dt\right)}{\mathbb{E}_{\vec{\partial}\mathcal{A}}T_A} \quad (2)$$

Note that, if we used Monte-Carlo to estimate this expectation by simulating the Markov process $X_t^\pi$ starting at $\vec{\partial}\mathcal{A}$ and observing cycle times $T_A$, it might take a very long time before observing a non-zero contribution from $\eta$, as this function is assumed to be sparse with non-zero values rarely observed. The following proposition is key to defining the Fleming-Viot estimation method leading to FVEE.
requiring fewer arrival events *T* for an accurate estimation of the denominator in (7), but at the same time visiting the rare blocking state *K* will be rarer, requiring a larger number of particles *N* for an accurate estimation of the numerator in (7). The opposite is true for larger values of *J*. A detailed analysis on this trade-off in terms of appropriately choosing *N* and *T* for accurate estimations of the numerator and the denominator in (7) is presented in Appendix C.

We now study the convergence of the FV estimator as either *N* or *T* increases, and compare it with the benchmark MC estimator, obtained from a direct application of expression (1), i.e. as the fraction of the time spent at state *K* and the total

time of return cycles to the initial state *x* = *J* − 1 observed during the simulation. To guarantee a fair comparison between the two methods, we start the simulation at *x* = *J* − 1, so that both methods start at the same "distance" from the blocking state *K*, and let the simulation run until the same number of events observed in the FV estimator is reached.

Figure 2 compares the convergence of the FV estimator with the MC estimator both as *N* increases and as *T* increases. We considered the cases *K* = 20 and *K* = 40, which are regarded to represent moderate and large capacities based on their blocking probabilities at the considered value for *ρ* = 0.7 of order 10⁻⁴ and 10⁻⁷, respectively. The size *J* of the absorption set *A* is held fixed at *J* = 12.

<figure><img src="image_3.png" /><figcaption>Figure 2: Violin plots showing the convergence properties of the Fleming-Viot (green) and Monte-Carlo (red) estimators of the blocking probability for an M/M/1/K queue system with $\lambda = \rho = 0.7$. Left plots (a) and (c): convergence as the number of particles *N* increases with *T* fixed. Right plots (b) and (d): convergence as the number of arrival events *T* increases with *N* fixed. The absorption set size is *J* = 12 in all cases. For the left set of plots (a) and (c), the number of particles considered for the convergence analysis are *N* = 264, 1055, 4220 for *K* = 20, and *N* = 763, 1356, 3051 for *K* = 40. For the right set of plots (b) and (d), the number of arrival events considered for the convergence analysis are *T* = 462, 1040, for both *K* = 20 and *K* = 40. For more details on these choices, see Appendix C. The middle horizontal line in each violin plot represents the median. In the corresponding MC experiments, the average number of observed return cycles to the initial state *J* − 1 is used on the horizontal axis while the top horizontal axes show the average number of events observed in the experiments run in each set, which by design coincide between each paired FV-MC execution, as described in the text.</figcaption></figure>
**Proposition 1.** Given a set $\mathcal{A} \subset S$ and a function $\eta: S \to \mathbb{R}$ that is non-zero on $C \in \mathcal{A}^c$, the expectation $\mathbb{E}^\pi(\eta)$ is given by (1), with

$$p^\pi(x) = \frac{\int_0^\infty \mathbb{P}_{\vec{\partial}\mathcal{A}^c}(T_K > t) \phi_t^{\vec{\partial}\mathcal{A}^c}(x) dt}{\mathbb{E}_{\vec{\partial}\mathcal{A}} T_A} \quad \forall x \in C, \qquad (3)$$

where $\phi_t^{\vec{\partial}\mathcal{A}^c}(x) \doteq \mathbb{P}_{\vec{\partial}\mathcal{A}^c}(X_t^\pi = x | T_K > t)$ is the probability that the process $X_t^\pi$ (started at a state in $\vec{\partial}\mathcal{A}^c$ chosen with probability $p_{\vec{\partial}\mathcal{A}^c}^\pi$) is in $x$ provided it has not been absorbed into $\mathcal{A}$.

*Proof.* Since at the beginning of a cycle the process starts at a state in $\vec{\partial}\mathcal{A}$, in order to reach a state in the set of interest $C \in \mathcal{A}^c$ (the set of states with non-zero values of $\eta$), the process needs to go through $\vec{\partial}\mathcal{A}^c$.

Thus, using the definition of the sojourn times $T_{\mathcal{A}^c}$ and $T_A$ and recalling that $T_A > T_{\mathcal{A}^c}$, expression (2) can be written as

$$\mathbb{E}^\pi(\eta) = \frac{\mathbb{E}_{\vec{\partial}\mathcal{A}}\left(\int_0^{T_{\mathcal{A}^c}} \eta(X_t^\pi) dt + \int_{T_{\mathcal{A}^c}}^\infty \eta(X_t^\pi) \mathbf{1}_{T_A > t} dt\right)}{\mathbb{E}_{\vec{\partial}\mathcal{A}} T_A}.$$

Since by assumption $\eta$ is zero in $\mathcal{A}$, the first integral is zero. After changing the integration variable to $u = t - T_{\mathcal{A}^c}$ in the second integral and using the notation of the killing time $T_K = T_A - T_{\mathcal{A}^c}$, we get

$$\mathbb{E}^\pi(\eta) = \frac{\mathbb{E}_{\vec{\partial}\mathcal{A}}\left(\int_0^\infty \eta(X_{u+T_{\mathcal{A}^c}}^\pi) \mathbf{1}_{T_K > u} du\right)}{\mathbb{E}_{\vec{\partial}\mathcal{A}} T_A}.$$

At $u = 0$, the Markov process is at $X_{T_{\mathcal{A}^c}}^\pi$ which, as stated in Section 3.1, is distributed according to $p_{\vec{\partial}\mathcal{A}^c}^\pi$ when $X_0^\pi \sim p_{\vec{\partial}\mathcal{A}^c}^\pi$, as is the case above. This allows us to redefine the time origin of the Markov process at $u = 0$ and replace $\mathbb{E}_{\vec{\partial}\mathcal{A}}$ with $\mathbb{E}_{\vec{\partial}\mathcal{A}^c}$ to obtain

$$\mathbb{E}^\pi(\eta) = \frac{\mathbb{E}_{\vec{\partial}\mathcal{A}^c}\left(\int_0^\infty \eta(X_u^\pi) \mathbf{1}_{T_K > u} du\right)}{\mathbb{E}_{\vec{\partial}\mathcal{A}} T_A}.$$

Then, expression (1) with $p^\pi(x)$ given by (3) follows from interchanging the order of the integral and the expectation (the state space is assumed finite), conditioning on the event $\mathbf{1}_{T_K > t}$, and taking the summation over $x$ (coming from the expectation in the numerator) out of the integral. ■

When $\eta(x)$ is a known function, the FVEE of $\eta$ is constructed as

$$\text{FVEE}(\eta) \doteq \hat{\mathbb{E}}_{FV}^\pi(\eta) \doteq \sum_{x \in C} \hat{p}_{FV}^\pi(x) \eta(x), \qquad (4)$$

where $\hat{p}_{FV}^\pi(x)$ is the FV estimator of the stationary probability for all states $x$ in $C$. This estimator is constructed by estimating the denominator and the numerator in (3) using the following two-step procedure: (i) the denominator $\mathbb{E}_{\vec{\partial}\mathcal{A}} T_A$ is estimated using a standard Monte-Carlo estimator from observations of the return cycle time $T_A$ coming from the simulation of $X_t^\pi$, and (ii) for each $x \in C$, the numerator is estimated by summing the product of the finite-support piecewise-constant estimators of the two functions in the integral, $\mathbb{P}_{\vec{\partial}\mathcal{A}^c}(T_K > t)$ and $\phi_t^{\vec{\partial}\mathcal{A}^c}(x)$, obtained from the simulation of the FV $N$-particle system driven by $X_t^\pi$ with absorption set $\mathcal{A}$. The estimation details are given in Appendix A.

If only the support of $\eta(x)$ is known but not its values, an estimator $\hat{\eta}(x)$ must be obtained separately from $\hat{p}_{FV}^\pi(x)$ and then plugged into the expression for the FVEE($\eta$) in (4). This will be the case when estimating gradients for the FVRL algorithm described in Section 3.3.1.
We observe the following in terms of convergence of the FV and MC estimators to the true blocking probability, which in the plots is represented by a horizontal gray line:

1. Both the FV and the MC estimators converge to the true blocking probability when $K = 20$ (Figures 2(a,b)). MC presents a smaller variability than FV in the convergence analysis with $N$ (Figure 2(a)) but a larger variability than FV in the convergence analysis with $T$ (Figure 2(b)). This is due to the fact that Figure 2(a) is obtained from simulations whose events largely outnumber the events in Figure 2(b) by as much as 10 times ($\sim 20,000$ vs. $\sim 2,000$), which allows MC to observe the blocking event at $x = K$ frequently enough for an accurate estimation of the blocking probability. On the other hand, Figure 2(b) tells us that a much smaller number of events ($\sim 2,000$) is sufficient for FV to estimate the blocking probability accurately enough, but is not sufficient for MC. This demonstrates the higher efficiency of FV than MC in discovering the rare event at $K$.

2. When $K = 40$ (Figures 2(c,d)) the MC estimator basically fails as it almost never observes the blocking state $K$.

3. Increasing $N$ has a larger impact in increasing the FV estimator's computational complexity than increasing $T$. This conclusion is obtained by observing that the number of events at the leftmost violin plot in Figure 2(a) ($\sim 20,000$ events), where $N \sim 4,000$ and $T \sim 100$, is about 10 times larger than the number of events at the rightmost violin plot in Figure 2(b) ($\sim 3,000$ events), where $N \sim 100$ and $T \sim 4,000$ (i.e. the $N$ and $T$ values are switched w.r.t. the previous case).

### 4.1.2 Learning the optimum blocking size using FVRL

We illustrate the FVRL learning algorithm by defining a cost function of rejecting an incoming job that is exponentially increasing with the queue size at the time of the job arrival. Note that the blocking cost needs to be exponentially increasing with the queue size in order to obtain an optimisation problem with a non-trivial solution (i.e. where the optimum $K$ is finite). This is due to the fact that the stationary probability is exponentially decreasing with the queue size, as mentioned at the end of Section 4.1, and such exponentially increasing cost function guarantees that the expected cost to minimize is a convex function of $K$. More precisely, given the queue size $x$, the cost function is defined as $r(x, a) = B(1 + b^{x-x_{ref}})1_{\{a=0\}}$, where $B, b$ and $x_{ref}$ are positive constants which, for a given load $\rho$, are carefully chosen to make the expected cost to minimize a convex function. From [14] we know that the policy optimising such convex cost functions is of threshold type, that is, there exists a state $0 < K < \infty$ such that $a = 1$ is optimal for all $x < K$, and $a = 0$ is optimal for $x = K$. The value of $x_{ref}$ is the reference queue size that is instrumental in defining the optimum threshold, $K^*$, which is the closest integer to $x_{ref}$.

We now use the policy gradient methodology presented in Section 3.4 to learn the optimum threshold $K^*$. Following [16] and using the framework of Semi-Markov Decision Processes (SMDP) to describe the continuous-time Markov process $X_t^\pi$ by a discrete-time process $X_n^\pi$ [19, 16], we propose a parameterised acceptance policy $\pi(a = 1|x)$ that is a linear step function of the state $x$, which is deterministic for $x$ outside the interval $(\theta, \theta+1)$ and decreases linearly from 1 to 0 within the interval. That is, the acceptance policy parameterised by the positive real-valued $\theta$, is defined as:

$$ \pi_{\theta}(a = 1|x) = \begin{cases} 1 & \text{if } x \le \theta, \\ x - \theta + 1 & \text{if } \theta < x < \theta + 1, \\ 0 & \text{if } x \ge \theta + 1. \end{cases} \quad (8) $$

Note that the policy is deterministic for integer-valued $\theta$, in which case the blocking size is $K = \theta + 1$. Otherwise, $K$ is defined as $\lceil \theta \rceil + 1$.

We use a gradient descent algorithm to learn the optimum parameter $\theta$ that minimizes the average state value $v^{\pi\theta}$, where "value" here is thought of as cost. Using expression (6), the gradient of $v^{\pi\theta}$ becomes

$$ \frac{\partial v^{\pi\theta}}{\partial \theta} = p^{\pi\theta}(K-1)[Q^{\pi\theta}(K-1, 1) - Q^{\pi\theta}(K-1, 0)], \quad (9) $$

where $K-1$ is the smallest integer that is larger than or equal to $\theta$. We observe that this parameterisation leads, as expected, to gradients being 0 for $x \ne K-1$, that is, the policy gradient is sparse. We note also that the gradient is discontinuous at $\theta$ and $\theta+1$, making the assumptions of Proposition 3 not fully satisfied. However, these two points have measure zero and therefore, with probability 1, no discontinuity is observed$^3$.

Details of the learning algorithm are given in Algorithm 2 in Appendix D.

To illustrate the FVRL algorithm, we consider an MDP with the following characteristics: the system load is $\rho = 0.7$, the blocking cost function $r(x, a)$ is defined with the parameters $b = 3$ ($> 1/\rho$ so that the expected cost function is convex), $B = 5$, and $x_{ref} = 18$, giving $K^* = 18$.

The setup of the learning experiments is as follows: we choose the value $(x_{ref} + 10)$ for the initial blocking size guess, so that, already at the onset, blocking occurs rarely. Since the value of $K$ is no longer fixed (as was the case when estimating the blocking probability) but is now learned by the algorithm, it is not possible to choose a fixed value $J$ for the size of the absorption set $A$. Instead, we consider a fixed $J/K$ fraction that adapts $J$ to each value of $K$ at the start of each learning step. In order to experiment with different sizes of the absorption set $A$, we consider two different scenarios: $J = \lceil 0.3K \rceil$ and $J = \lceil 0.5K \rceil$. Following the conclusions of the choice of $J, N$ and $T$ and their impact in the FV estimator accuracy outlined in Appendix C, for each value of $J$, we adjust $N$ and $T$ at each learning step to obtain approximate expected relative errors of $\epsilon_\phi \sim 100\%$ and $\epsilon_{ET} \sim 20\%$, respectively for

<footnote>$^3$A special case occurs when $\theta$ is integer, in which case the discontinuities would be observed with non-zero probability in the gradient descent algorithm under the following scenario: an integer value is chosen for the initial guess of $\theta$, and integer-valued clipping (e.g. to $\pm 1$) is used for the next $\theta$ estimated by the algorithm. This problem is solved by simply not choosing an integer-valued initial guess of $\theta$.</footnote>
**Consistency of FVEE(η)** It has been proved that for finite state spaces [8, 9], uniform in time propagation of chaos holds for Fleming-Viot particle systems. We let $m(\cdot, \xi) : \mathcal{A}^c \to [0, 1]$ denote the empirical distribution of the $N$ particles with positions described by vector $\xi$, defined as the empirical mean $m(x, \xi) \doteq \frac{1}{N} \sum_{i=1}^{N} \mathbf{1}_{\xi(i)=x}, \forall x \in \mathcal{A}^c$. When $\nu$ is a probability measure on $\mathcal{S}$, under the assumption of proper initialization [8, Theorem 1.4] shows the following bound on the speed of convergence w.r.t. the number of particles $N$ of the empirical mean $m(\cdot, \xi_t^\nu)$ towards $\phi_t^\nu$:

$$ \sup_{x \in \mathcal{S}} \sup_{t \ge 0} \mathbb{E} \left| m(x, \xi_t^\nu) - \phi_t^\nu(x) \right| \le \frac{C_{\mathrm{FV}}}{\sqrt{N}}, \quad (5) $$

where $C_{\mathrm{FV}}$ is a positive constant depending on the characteristics of the driving process.

Using this result, we can show the following bound for the expected absolute error for the FV estimator of the stationary probability, $\hat{p}_{\mathrm{FV}}^{\pi}(x)$. The theorem is valid in an idealized case where the simulation used to estimate the denominator $\mathbb{E}_{\vec{\partial}_{\mathcal{A}}T_{\mathcal{A}}}$ is started according to $p_{\vec{\partial}_{\mathcal{A}}}^{\pi}$, and the Fleming-Viot simulation used to estimate $\mathbb{P}_{\vec{\partial}_{\mathcal{A}}^{c}}(T_{\mathcal{K}} > t)$ and $\phi_t^{\vec{\partial}_{\mathcal{A}}^{c}}(x)$, is started from i.i.d. samples of $p_{\vec{\partial}_{\mathcal{A}}^{c}}^{\pi}$.

**Theorem 2** (Bound for the expected absolute error of $\hat{p}_{\mathrm{FV}}^{\pi}(x)$). Assume that we start the simulation of $X_t^\pi$ for the estimator of $\mathbb{E}_{\vec{\partial}_{\mathcal{A}}T_{\mathcal{A}}}$ in (3) according to the distribution $p_{\vec{\partial}_{\mathcal{A}}}^{\pi}$, that $M$ return cycles to $\mathcal{A}$ under stationarity are observed during that simulation, and that we compute the estimator of $\mathbb{P}_{\vec{\partial}_{\mathcal{A}}^{c}}(T_{\mathcal{K}} > t)$ and $\phi_t^{\vec{\partial}_{\mathcal{A}}^{c}}(x)$ in (3) using the FV particle system started at the positions of $N$ i.i.d samples from $p_{\vec{\partial}_{\mathcal{A}}^{c}}^{\pi}$. Then, there exists a constant $C > 0$ such that

$$ \mathbb{E} |\hat{p}_{\mathrm{FV}}^{\pi}(x) - p^{\pi}(x)| \le C \left( \frac{1}{\sqrt{M}} + \frac{1}{\sqrt{N}} \right), \quad \forall x \in \mathcal{A}^c. $$

The proof is given in Appendix B.

Theorem 2 shows that the FV estimator of the stationary probability, $\hat{p}_{\mathrm{FV}}^{\pi}(x)$, is consistent. And when function $\eta(x)$ is known, it implies that the FVEE($\eta$) defined in (4) is consistent as well, since it is an $\eta$-weighted sum of the $\hat{p}_{\mathrm{FV}}^{\pi}(x)$ values on the finite set $\mathcal{C}$.

A few additional comments on the theorem are in order. On the one hand, it is important to observe that Theorem 2 ensures the estimation error is of the same order as Monte-Carlo, which gives minimal convergence properties to FVEE. The constant $C$ depends in an intricate manner on the state and action space sizes and on the precise underlying Markovian dynamics (see [8] for more details). It remains an open problem to characterize for which type of dynamics/rewards, this constant and hence the global error will be smaller than the corresponding ones for vanilla Monte-Carlo.

On the other hand, the mean square estimation error should not be our only focus to evaluate the difference between FVEE and Monte-Carlo, especially when the ultimate goal is the convergence rate of a reinforcement learning algorithm. Indeed, in the control problem, a noisy but still informative signal might be very useful compared to no signal at all. Our main idea when replacing MC by FV is to trade the observation of a very rare event in the original problem by the observation of a more common event for the FV particle system.

Although a fully rigorous analysis of the large deviation probability to observe a non-zero reward is out of the scope of this paper (see [11] for a recent analysis in the context of random walks and Lévy processes), we can give some intuition using the results of [8, 12]. The dynamics of a tagged particle of the FV process converges when $N \to \infty$ to a one-dimensional Markov process having the quasi-stationary distribution $\nu_{QS}$ of the original absorbed process as stationary measure (see [12]), where

$$ \nu_{QS}(B) = \lim_{t \to \infty} \phi_t^{\vec{\partial}_{\mathcal{A}}^{c}}(B), $$

for any set of states $B$. If the state space is finite, this one dimensional process in turn converges in distribution exponentially fast to its stationary distribution, $\nu_{QS}$. Hence, the probability of finding a non-zero signal (by visiting
the estimations of $\phi_t^J(K-1)$ and $\mathbb{E}_{J-1}(T_A)$, i.e. we set the approximate expected error for $\hat{\phi}_t^J(K-1)$ much larger than the approximate expected error for $\hat{\mathbb{E}}_{J-1}(T_A)$. For comparison purposes and to further check the conclusions of the error analysis in Appendix C, we also considered a scenario where the pre-defined relative errors are inverted, namely $\epsilon_\phi \sim 20\%$ and $\epsilon_{ET} \sim 100\%$.

For each setup we ran the FVRL policy learner on 100 learning steps. At each learning step, as long as a positive estimate of the stationary probability $p^{\pi^\theta}(K-1)$ has been obtained, we also estimate the sparse $\eta(x)$ function defined in Section 3.4—which in this case is simply the difference between two $Q$ values as seen in (9)—using the coupling procedure described in Section 3.4 and allowing up to 250 arrival events until mixing is observed in each of 100 replications. The value of $\hat{\eta}(x)$ is finally computed as the average $Q$-difference over these replications<fnref n="4" /> and multiplied with the stationary probability estimate $\hat{p}^{\pi^\theta}(K-1)$. Parameter $\theta$ is then updated by gradient descent using a constant learning rate of $\alpha = 10$, and the number of observed events at each learning step is recorded. The estimated optimum threshold is set to $\hat{K}^* = \text{round}(\hat{\theta}^*) + 1$, where $\hat{\theta}^*$ is the value of the $\theta$ parameter obtained after the last learning step.<fnref n="5" />

The benchmark MC learning algorithm uses the same policy gradient approach as FVRL with the only difference that it estimates the probability $p^{\pi^\theta}(K-1)$ in expression (9) using Monte-Carlo instead of Fleming-Viot, i.e. based on a single trajectory of $X_t^\pi$ as the ratio between the continuous time that the system spends at $K-1$ and the total return cycle time to $J-1$ under stationarity<fnref n="6" />, where $J$ is defined as a function of $K$ as in FVRL. Each replication of the MC learner is started at $J-1$ and is stopped when the average number of observed events over all learning steps in the respective FVRL replication is observed. These two conditions allow a fair comparison

between the two methods.

The results of the above procedure run on 20 replications are shown in Figure 3. We see that FVRL clearly outperforms the benchmark MC learning in three out of the four scenarios considered, depicted in Figures 3(a), 3(b) and 3(d). In the scenario depicted in Figure 3(c), learning by the two methods is similar because the average number of events is about 10 times larger than the one observed in the other three scenarios (~ 40, 000 vs. ~ 5, 000), which allows MC to also observe the rare blocking state and thus learn almost as fast as FVRL. This is the same situation previously observed about the probability estimation results presented in Figure 2.

Of the two setups, $\epsilon_\phi \sim 100\%$ and $\epsilon_{ET} \sim 20\%$ used for Figures 3(a) and 3(b), and $\epsilon_\phi \sim 20\%$ and $\epsilon_{ET} \sim 100\%$ used for Figures 3(c) and 3(d), the setup that is both the safest and the least computationally intensive is the former. This is consistent with the conclusions of the analysis of the choice of $J$, $N$ and $T$ in FV presented in Appendix C, namely that, for an accurate estimation of the stationary probability, more importance should be given to achieve an accurate estimation of $\mathbb{E}_{J-1}(T_A)$ than an accurate estimation of $\phi_t^J(K)$. In terms of safety we note that, for the least convenient setup of $\epsilon_\phi = 20\%$ and $\epsilon_{ET} = 100\%$ of Figure 3(d), in one replication the $\theta$ parameter suddenly increases from 28 to 140 (not shown but made apparent by the mean learning curve being significantly above the median learning curve). This overshoot impedes further learning because the blocking probability becomes extremely small at $K=140$. This clearly illustrates the risks of choosing a too small $T$ value (associated to the large error $\epsilon_{ET} = 100\%$) which may considerably overestimate the blocking probability (due to an underestimation of $\mathbb{E}_{J-1}(T_A)$) and generate such out-of-control excursion coming from a noisy estimation of the $Q$-difference contributing to the gradient expression.

## 4.2 Loss network system

Consider now a loss network that processes jobs of $I$ different classes which, as in the $M/M/1/K$ system, arrive following independent Poisson processes and are served with independently exponentially distributed times by one of $R$ available servers. We are interested in analyzing the number of jobs of each class being served by the system at a given time $t$, $X_t^\pi = (X_{t,1}^\pi, \dots, X_{t,I}^\pi)$, whose dynamics is governed by the system characteristics and the policy $\pi$ being applied. As before, we also consider a policy $\pi$ that belongs to the family of deterministic accept/reject policies that reject an incoming job of a given class when the system is already serving a predefined number of jobs of that class, in which case a class-dependent rejection cost is accrued.

Let us denote by $\lambda_i, \mu_i, C_i, K_i$ the job arrival rate, the service rate, the rejection cost, and the blocking size of each job class $i = 1, \dots, I$. Thus, the system

can be represented by an MDP ($S = \{X \in \mathbb{R}^I : 0 \le X_i \le K_i, i = 1, \dots, I \text{ s.t. } X^T1 \le R, A = \{0,1\}, q = \{\lambda_1, \dots, \lambda_I; \mu_1, \dots, \mu_I\}, R = r(x, a) = \sum_{i=1}^I C_i 1_{\{x \in S_i, \mathcal{L}_i, a=0\}}$), where $S_i$ is the set of states where blocking occurs following an $i$-class job arrival, an event that is denoted by $\mathcal{L}_i$.

For the purposes of simplifying illustration, in our experiments we consider the smallest multi-class loss network, one that serves just two classes of job.

### 4.2.1 Estimation of the expected cost

Contrary to the $M/M/1/K$ case, estimating the expected cost in the loss network system is not equivalent to estimating the blocking probability, because blocking occurs at more than one state whose cost in general depends on the class of the arriving job being rejected.

The expression of the expected cost $\mathbb{E}^\pi(\mathcal{R})$ under the

<footnote n="4">⁴The number of replications on which the $Q$-difference is averaged may be less than 100 because a replication is excluded from the average when mixing does not occur, but this occurs very rarely.</footnote>

<footnote n="5">⁵Note that the estimated optimum threshold is not set to $[\hat{\theta}^*] + 1$, as in the definition of the parameterised policy, so that $\hat{K}^*$ is more naturally chosen to be e.g. 5 when $\hat{\theta}^* = 4.01$ rather than 6.</footnote>

<footnote n="6">⁶Stationary is assumed after the burn-in period whose value is defined in Section 4.</footnote>
states in set $\mathcal{C}$) for the FV process in a finite time interval is of the order of $\nu_{QS}(\mathcal{C})$, which can be considerably larger than the original $p^{\pi}(\mathcal{C})$.

As an illustration, we bring back the $M/M/1$ queue example with load $\rho < 1$ and admission control policy that rejects an incoming job when the queue size is $K$. Thus $\mathcal{C} = \{K\}$. If the absorption set $\mathcal{A}$ is chosen as the singleton containing the high probability state of zero jobs in the queue, the quasi-stationary distribution of the absorbed process has $\nu_{QS}(K) \sim \sqrt{K}\rho^{K/2}$ [10], which is considerably larger than the stationary probability of the original process at $K$, given by $p^{\pi}(K) \sim \rho^K$. For example, $\nu_{QS}(K)$ is about 8 thousand times larger than $p^{\pi}(K)$ when $K = 40$ and $\rho = 0.7$, and 6 million times larger when $K = 40$ and $\rho = 0.5$.

### 3.3 FVRL: Policy gradient learning with Fleming-Viot particle systems

In this subsection we show how the Fleming-Viot expectation estimator (FVEE) introduced in Section 3.2 can be combined with the policy gradient algorithm to approximately solve optimal control problems in environments with sparse and rare rewards under the average reward criterion. This leads to FVRL, the Fleming-Viot reinforcement learning method.

We provide here a general description of the FVRL method to speed up learning in those types of scenarios. The methodology will be illustrated on two network optimisation problems described in Section 4: an $M/M/1$ queue, and a loss network system serving two classes of jobs.

Let us consider in detail the case in which an agent interacts with an environment (available in practice through a simulator or emulator), with the aim of maximizing the long-run expected reward. Following the customary RL framework, we consider a discrete-time representation $X_n^\pi$ of the continuous-time process $X_t^\pi$ such that both processes have the same stationary probability distribution. This can be achieved by different methods, such as uniformization [22, Chapter 11.5] or, in the case of queueing systems with Poisson arrivals, by defining $X_n^\pi$ as the chain embedded at arrival times [5] (using PASTA). In our queueing examples we will rely on the latter.

As customary in the literature, we let $X_n$, $A_n$ and $R_{n+1}$ respectively denote the state, action and reward at the $n$-th discrete time step, $n \ge 0$. We denote by $\pi_\theta$ the policy function parameterised by $\theta \in \mathbb{R}^d$, and we propose to use a gradient-based algorithm to learn the parameter $\theta$ that maximizes $g(\pi_\theta)$, the long-run expected reward [25, Chapter 13],

$$ 
\begin{aligned}
g(\pi_{\theta}) &= \lim_{H \to \infty} \frac{1}{H} \sum_{n=1}^{H} \mathbb{E}^{\pi_{\theta}}[R_n | X_0, A_{0:n-1} \sim \pi_{\theta}] \\
&= \sum_{x \in S} p^{\pi_{\theta}}(x) \sum_{a \in \mathcal{A}} r(x, a) \pi_{\theta}(a|x),
\end{aligned} 
\quad (6) 
$$

where the last equality comes from the ergodicity assumption of the MDP under all parameterised $\pi_\theta$ policies considered, making the long-run expected reward independent of the start state $X_0$.

It follows from the policy gradient theorem for continuing learning problems [25, 26] that the gradient of $g(\pi_\theta)$ can be written as:

$$ 
\begin{aligned}
\nabla_{\theta} g(\pi_{\theta}) &= \sum_{x \in S} p^{\pi_{\theta}}(x) \sum_{a \in \mathcal{A}} Q^{\pi_{\theta}}(x, a) \nabla_{\theta} \pi_{\theta}(a|x) \\
&= \sum_{x \in S} p^{\pi_{\theta}}(x) \nabla_{\theta} g_x(\pi_{\theta}),
\end{aligned} 
\quad (7) 
$$

where we have defined the contribution to the gradient by state $x$ as

$$ \nabla_{\theta} g_x(\pi_{\theta}) \doteq \sum_{a \in \mathcal{A}} Q^{\pi_{\theta}}(x, a) \nabla_{\theta} \pi_{\theta}(a|x), \quad (8) $$
<figure><img src="image_4.png" /><figcaption>Figure 3: Comparison of FVRL learning (left subplot in green) with Monte-Carlo learning (right subplot in red) of the optimum parameter $\theta^* = 17$ (horizontal dashed gray line) of the parameterised acceptance policy $\pi_\theta$ of incoming jobs to the $M/M/1/K$ queue system with $\lambda = \rho = 0.7$. Each plot caption indicates the following learning settings for the FVRL approach: the approximate expected relative error wished for the estimation of $\phi_t(K)$ (which defines the number of particles $N$ at each learning step), the approximate expected relative error wished for the estimation of $\hat{E}(T_A)$ (which defines the number of arrival events $T$ at each learning step), and the fraction $J/K$ defining the size of the absorption set. In all cases the learning parameter is kept constant at $\alpha = 10$, and the initial $\theta$ guess is 28.1, which defines a "large" blocking size of $K = 30$ whose small stationary probability ($\sim 10^{-6}$) makes blocking a rare event already at the onset of learning.</figcaption></figure>

threshold policy $\pi$ and rewards $\mathcal{R}$ is obtained from the total probability theorem applied on all possible job arrival classes and all possible accept/reject actions on those arrivals, as:

$$E^\pi(\mathcal{R}) = \sum_{x \in S} p^\pi(x)c(x), \quad (10)$$

where $p^\pi(x) \doteq Pr(\mathbf{X}_n^\pi = x)$ is the stationary probability that the SMDP $\mathbf{X}_n^\pi$ (defined in Section 4.1.2) is at state $x$, and $c(x) = \sum_{i=1}^I C_i \lambda_i / \Lambda [1_{\{\sum_{j=1}^I x_j = R\}} + 1_{\{\sum_{j=1}^I x_j < R\}} 1_{\{x_i = K_i\}}]$ is the expected cost of rejection at $x$ over all possible arriving job classes, where $\Lambda \doteq \sum_{i=1}^I \lambda_i$ is the total job arrival rate.

The true expected cost is computed from (10) using the stationary probability of the stochastic knapsack as $p^\pi(x)$ [21, Chapter 4], against which the estimated expected cost is compared.

The convergence of the estimated expected cost to its true value is analyzed via the violin plots shown in Figure 4, where, following the conclusions about the choice of $J$, $N$ and $T$ for the $M/M/1/K$ queue system described in Appendix C, is analyzed in terms of increasing $T$ (as opposed to increasing $N$), since a large value of $T$ is more crucial than a large value of $N$ for an accurate estimation. In order to analyze the impact on the estimation accuracy of the distribution of the start states, two scenarios are considered in this analysis: (a) one where the $N$ initial states for the FV simulation used to estimate $P_{\tilde{\partial}A^c}(T_K > t)$ and $\phi_t^{\tilde{\partial}A^c}$ are chosen uniformly at random out of the states in $\tilde{\partial}A^c$, and (b) one where the initial states are chosen according to the stationary distribution of the states in $\tilde{\partial}A^c$, which is known for the loss network considered.
and

$$Q^{\pi_{\theta}}(x, a) \doteq \mathbb{E}^{\pi_{\theta}} \left[ \sum_{n=1}^{\infty} (R_n - g(\pi_{\theta})) \mid X_0 = x, A_0 = a \right] \quad (9)$$

is the action value function for continuing learning problems under the average reward criterion.

We note that the gradient in (7) has the form of the expectation in (1) if we let $\eta(x) = \nabla_{\theta}g_x(\pi_{\theta})$ and define the set $C = \{x \in S : \nabla_{\theta}g_x(\pi_{\theta}) \neq 0\}$. Therefore, we can leverage FVEE to obtain an estimate of the gradient in (7) by combining the FV estimator of the stationary probability presented in Section 3.2, for states $x \in C$, with an estimate of $\nabla_{\theta}g_x(\pi_{\theta})$. This leads to the FVRL algorithm described in detail below. We note that the advantage of our approach lies on the fact that in a wide class of queueing problems, a smart choice of the parameterised policy leads to a sparse $\nabla_{\theta}g_x(\pi_{\theta})$ as a function of $x$. This is the case for instance in problems in which the optimal policy is of threshold type where, as we will see in Section 4.1.2, an appropriate parameterisation yields $\nabla_{\theta}g_x(\pi_{\theta}) = 0$ in all states $x$ except in a neighbourhood of the thresholds<fnref n="2" />.

### 3.3.1 FVRL algorithm

The FVRL algorithm is a gradient ascent algorithm that uses estimates of the gradient defined in (7) to learn a policy that is locally optimal in the parameter space. Assuming a policy parameterisation $\nabla_{\theta}\pi(a|x)$ with sparse gradients in $x$ and known support $C$, the gradient in (7) will also be sparse in $x$. Thus, we propose to estimate the gradient $\nabla_{\theta}g(\pi_{\theta})$ using (4), i.e. by estimating the stationary probability $p^{\pi_{\theta}}(x)$ for $x \in C$ using an FV particle system, as described in Section 3.2. However, contrary to the estimation problem presented in that section, the $\eta$ function in this case, $\eta(x) = \nabla_{\theta}g_x(\pi_{\theta})$ defined in (8), is only known up to its support $C$, as $Q^{\pi_{\theta}}(x, a)$ is unknown. This requires that $\nabla_{\theta}g_x(\pi_{\theta})$ be estimated as well (as explained below) to obtain the final estimator of the gradient as:

$$\hat{\nabla}_{\theta}g(\pi_{\theta}) = \sum_{x \in C} \hat{p}_{FV}^{\pi_{\theta}}(x) \hat{\nabla}_{\theta}g_x(\pi_{\theta}).$$

**Estimation of $\nabla_{\theta}g_x(\pi_{\theta})$** We propose an estimation procedure that leverages the sum-to-zero property of the policy gradient in combination with a coupling argument, as described in the following lemma.

**Lemma 3.** Given a state $x \in S$, let $\{X_n^{(x,a)}\}_{a \in U, n \in \mathbb{N}}$ be a family of $|U|$ identical Markov processes driven by policy $\pi_{\theta}(\cdot|x)$, all starting at $x$, and each taking one of the different possible actions $a \in U$ at time $n = 0$. Define the coupling time $\tau$ when all such chains meet, i.e. $\tau = \inf\{n > 0 : X_n^{(x,a)} = X_n^{(x,a')}, \forall a, a' \in U\}$. Then, the function $\nabla_{\theta}g_x(\pi_{\theta})$ of $x$, defined in (8), can be written as:

$$\nabla_{\theta}g_x(\pi_{\theta}) = \sum_{a \in U} \mathbb{E}^{\pi_{\theta}} \left[ \sum_{n=1}^{\tau} R_n |X_0 = x, A_0 = a \right] \nabla_{\theta}\pi_{\theta}(a|x), \quad \forall x \in S.$$

*Proof.* After the stopping time $\tau$, the contribution to $Q^{\pi_{\theta}}(x, a)$ in (9) does not depend anymore on the initial action $a$ taken by each chain $X_n^{(x,a)}$ because they all have the same expected trajectory under the fixed policy $\pi_{\theta}(\cdot|x)$. Hence,

<footnote n="2">²Even if the gradient is not sparse, it may be the case that important contributions to its value comes, via $\nabla_{\theta}g_x(\pi_{\theta})$, from states $x$ that are rarely observed. Failing to observe these states frequently enough may hinder the convergence speed of the optimisation process to their actual optimum.</footnote>
The latter is still not the correct distribution to use for the initial states, but is a better approximation than the uniform to the actual distribution that should be used according to Theorem 2, namely the entrance stationary distributions of the states in $\bar{\partial}A^c$. In practice, when the stationary distribution of the states in $\bar{\partial}A^c$ is unknown, one should use the entrance stationary distribution of the states in $\bar{\partial}A^c$, estimated using the procedure described in Appendix A and implemented by Algorithm 1, which is based on the single simulation of the Markov process that estimates $E_{\bar{\partial}A}(T_A)$. Note that, as stated in Appendix A, the initial state for the estimation of $E_{\bar{\partial}A}(T_A)$ is always chosen uniformly at random among the states in $\bar{\partial}A$, and after that the burn-in period is used so that the system is closer to stationarity before collecting the data for the estimation process.

Finally, to guarantee a fair comparison with the benchmark, the simulation for the MC estimator is started at a uniformly randomly chosen state in $\bar{\partial}A$–so that both methods start at a similar "distance" from the set of blocking states– and is let run

until the number of events of the FV simulation is observed.

A loss network with the following characteristics was considered for the estimation of the expected cost: capacity $R = 6$ servers, $\lambda = [1, 5]$, $\mu = [3.33, 50.0]$, hence $\rho = [0.3, 0.1]$, $C = [2.5 \times 10^3, 4.9 \times 10^6]$, $K = [4, 6]$. The choice of the rejection costs allows us to illustrate the benefits of the FV approach over MC as it makes a few of the smaller probability cost-generating states as important as those with larger probability in terms of their contribution to the system's expected cost. These contributions are listed in Table 1, together with their respective stationary probabilities which have been computed following the product form of the stationary probability of the stochastic knapsack described in [21, chapter 4]. We observe that the top three states in terms of probability contribute to about 65% of the expected cost, while the bottom four, with probability smaller than $10^{-6}$, contribute to as much as 35%, therefore obtaining accurate estimates of the probability of a few of those smaller probability states is important for an accurate estimation of the expected cost.

<table><caption>Table 1: Contribution to the expected cost $E^\pi(\mathcal{R})$ by each blocking state in the loss network system with $R = 6$, $\lambda = [1, 5]$, $\mu = [3.33, 50.0]$ ($\rho = [0.3, 0.1]$), $C = [2.5 \times 10^3, 4.9 \times 10^6]$, blocking sizes $K = [4, 6]$ defining policy $\pi$, sorted by decreasing probability $p^\pi(\mathbf{x})$.</caption><thead><tr><th>State x</th><th>Prob. p<sup>π</sup>(x)</th><th>Expected cost c(x)</th><th>p<sup>π</sup>(x)c(x)</th><th>% Exp. cost E<sup>π</sup>(R)</th></tr></thead><tbody><tr><td>[4, 0]</td><td>2.2 &times; 10<sup>-4</sup></td><td>0.416 &times; 10<sup>3</sup></td><td>0.094</td><td>1.3%</td></tr><tr><td>[4, 1]</td><td>2.3 &times; 10<sup>-5</sup></td><td>0.416 &times; 10<sup>3</sup></td><td>0.009</td><td>0.2%</td></tr><tr><td>[4, 2]</td><td>1.1 &times; 10<sup>-6</sup></td><td>4084 &times; 10<sup>3</sup></td><td>4.525</td><td>62.7%</td></tr><tr><td>[3, 3]</td><td>5.0 &times; 10<sup>-7</sup></td><td>4084 &times; 10<sup>3</sup></td><td>2.011</td><td>27.8%</td></tr><tr><td>[2, 4]</td><td>1.3 &times; 10<sup>-7</sup></td><td>4084 &times; 10<sup>3</sup></td><td>0.503</td><td>7.0%</td></tr><tr><td>[1, 5]</td><td>1.7 &times; 10<sup>-8</sup></td><td>4084 &times; 10<sup>3</sup></td><td>0.067</td><td>0.9%</td></tr><tr><td>[0, 6]</td><td>9.3 &times; 10<sup>-10</sup></td><td>4084 &times; 10<sup>3</sup></td><td>0.004</td><td>0.1%</td></tr><tr><td>Total</td><td></td><td></td><td>7.214</td><td>100.0%</td></tr></tbody></table>

Figure 4 compares the convergence of the FV estimator with the MC estimator of the expected cost of the system as the number of arrival events $T$ increases, while keeping a fixed $N$ set to 500 and an absorption set $A$ with size $J = [1, 2]$ in each job class dimension. We observe:

* a consistently smaller variance of the FV estimator compared to the MC estimator.

* a higher accuracy of the FV estimator (although with

larger variance) when the start state of the FV simulation is chosen following the stationary distribution of states in $\bar{\partial}A^c$ compared to the uniformly random choice case. In the latter case (Figure 4(a)) the FV estimates tend to underestimate the expected cost, as made apparent by the median value represented by the horizontal line in the middle of each violin plot.

* a tendency of the MC estimator to overestimate the expected cost.

### 4.2.2 Learning the optimum blocking sizes using FVRL

The FVRL algorithm of the optimum blocking sizes $K^* \in \mathbb{N}^I$ of a loss network serving $I$ jobs classes learns the optimum $\theta_i$ of $I$ parameterised acceptance policies, each of the form (8). For each $\theta_i$ observed during learning, the deterministic blocking size $K_i$ of each policy is defined as $K_i = [\theta_i] + 1$, but, as in the $M/M/1/K$ case, the estimated optimum blocking size obtained at the end of the learning process is defined as $\hat{K}_i^* = \text{round}(\hat{\theta}_i^*) + 1$.

Upon arrival of a job of class $i$, the respective $\pi_{\theta_i}$ acceptance policy is applied, making the system's acceptance policy equal to $\pi_{\theta}(a = 1|x) = \sum_{i=1}^{I} \pi_{\theta_i}(a = 1|x_i)1_{L_i}$, where $L_i$ denotes the arrival event of a class-i job. Its derivative w.r.t. $\theta_i$ is non-zero only at states $x$ for which $x_i = K_i - 1$ as long as they satisfy the $R$-server constraint for a possible job acceptance, $x^T1 < R$. At those states the derivative is equal to +1 for action $a = 1$ (accept) and -1 for action $a = 0$ (reject).
we can write $\nabla_{\theta}g_x(\pi_{\theta})$ in (8) as

$$
\begin{align*}
\nabla_{\theta} g_x(\pi_{\theta}) &= \sum_{a \in \mathcal{U}} \mathbb{E}^{\pi_{\theta}} \left[ \sum_{n=1}^{\infty} (R_n - g(\pi_{\theta})) |X_0 = x, A_0 = a \right] \nabla_{\theta} \pi_{\theta}(a|x) \\
&= \sum_{a \in \mathcal{U}} \mathbb{E}^{\pi_{\theta}} \left[ \sum_{n=1}^{\tau} (R_n - g(\pi_{\theta})) |X_0 = x, A_0 = a \right] \nabla_{\theta} \pi_{\theta}(a|x) \\
&\quad + \mathbb{E}^{\pi_{\theta}} \left[ \sum_{n=\tau+1}^{\infty} (R_n - g(\pi_{\theta})) |X_0 = x \right] \sum_{a \in \mathcal{U}} \nabla_{\theta} \pi_{\theta}(a|x) \\
&= \sum_{a \in \mathcal{U}} \mathbb{E}^{\pi_{\theta}} \left[ \sum_{n=1}^{\tau} R_n |X_0 = x, A_0 = a \right] \nabla_{\theta} \pi_{\theta}(a|x) - \mathbb{E}^{\pi_{\theta}} [\tau] g(\pi_{\theta}) \sum_{a \in \mathcal{U}} \nabla_{\theta} \pi_{\theta}(a|x) \\
&= \sum_{a \in \mathcal{U}} \mathbb{E}^{\pi_{\theta}} \left[ \sum_{n=1}^{\tau} R_n |X_0 = x, A_0 = a \right] \nabla_{\theta} \pi_{\theta}(a|x),
\end{align*}
$$

where we have repeatedly used the sum-to-zero property of the policy gradient, namely $\sum_{a \in \mathcal{U}} \nabla_{\theta} \pi_{\theta}(a|x) = \nabla \sum_{a \in \mathcal{U}} \pi_{\theta}(a|x) = \nabla 1 = 0, \forall x \in S$. ■

Using Lemma 3, we propose an estimator of $\nabla_{\theta}g_x(\pi_{\theta})$ based on simulating the $|\mathcal{U}|$ copies of the process under policy $\pi_{\theta}$ as described by the lemma. When any two copies meet at the same state, we impose they continue evolving together forever. The simulation stops when either all copies meet or when a maximum number of steps taken by each copy is reached. Note that, due to Doeblin's theorem on coupling, [5, Section 4.3.1], the coupling time $\tau$ of all processes is almost surely finite and actually has an exponential moment when the state space is finite.

We remark that this estimation procedure can only improve over the traditional estimation procedure of independently simulating separate process copies for a predefined number of steps $S$. In fact, should the coupling time of all processes be smaller than $S$, the simulation will terminate earlier than the traditional method which only terminates after each copy has performed $S$ steps.

**The algorithm** The FVRL policy gradient algorithm is summarized in Algorithm 1.
<figure><img src="image_5.png" /><figcaption>(a) The start state of each FV particle is chosen uniformly at random in the boundary set $A^c$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) The start state of each FV particle is chosen following the stationary distribution of the states in the boundary set $A^c$</figcaption></figure>

Figure 4: Violin plots showing the convergence properties of the Fleming-Viot (left) and Monte-Carlo (right) estimators of the expected cost in a two-class loss network with $R = 6$ servers, $\lambda = [1, 5]$, $\rho = [0.3, 0.1]$, blocking at $K = [4, 6]$ with costs $C = [2.5 \times 10^3, 4.9 \times 10^6]$, as the number of arrival events $T$ increases. The absorption set sizes by job class are $J = [1, 2]$ and the number of particles is $N = 500$. The middle horizontal line in each violin plot is the median estimated value among 10 replications. In the corresponding MC experiments, the average number of observed return cycles to the absorption set is used on the horizontal axis. Finally, the top horizontal axes show the average number of events observed in the experiments run in each set, which by design coincide between each paired FV-MC execution, as mentioned in the text.

Hence, the partial derivative of $v^{\pi\theta}$ w.r.t. $\theta_i$ becomes:

$$ \frac{\partial v^{\pi\theta}}{\partial \theta_i} = \sum_{\mathbf{x} \in S_i} p^{\pi\theta}(\mathbf{x}) [Q^{\pi\theta}(\mathbf{x}, 1) - Q^{\pi\theta}(\mathbf{x}, 0)], \text{ for } i = 1, \dots, \text{ and } \mathbf{x}^T \mathbf{1} < R. $$

where $S_i$ is the set of states $\mathbf{x}$ where $x_i = K_i - 1$ and $\mathbf{x}^T \mathbf{1} < R$.

To illustrate the FVRL algorithm, we consider a two-job-class loss network with the following characteristics: $R = 6$ servers, job arrival rate by class $\lambda = [1, 5]$, service rate by class $\mu = [2.0, 16.67]$ (resulting in a load by class $\rho = [0.5, 0.3]$), and a blocking cost by class $C = [2 \times 10^3, 2 \times 10^4]$. For this choice of parameters, the optimum blocking sizes of the threshold policy are equal to $K^* = [4, 6]$ (i.e. $\theta^* = [3, 5]$), which will be estimated with FVRL.

The setup of the learning experiments is as follows: we set the initial blocking size guess to $K = [0.1, 0.1]$, and the following simulation parameters for the estimation of the stationary probabilities $p^{\pi\theta}(\mathbf{x})$ for each $\mathbf{x} \in S_i: J/K = 0.5$ for both classes, the number of FV particles $N = 100$, and the number of arrival events $T = 500$. The start state of each FV particle is chosen following the known stationary distribution of the states in $\partial A^{c(7)}$. The estimation of the gradient in (11) is completed by estimating the $Q$-difference at each $\mathbf{x}$ following the same coupling procedure used for the $M/M/1/K$ case. The learning rate $\alpha$ is set initially to 1.0 and is then decreased inversely proportional to the learning step.

The results of the above procedure run on 20 replications are shown in Figure 5 and described in the caption. We see that FVRL outperforms MC learning in terms of median learning

speed of parameter $\theta_1$ with non-trivial optimum, and presents a smaller path variability as well. On the other hand, the FVRL and MC average learning speeds are similar and are a little larger than the ideal learning scenario depicted in blue on the left, where the learning process uses the true stationary probabilities in the gradient expression, i.e. only the $Q$-differences are estimated for the computation of the gradient. The FVRL algorithm is closer than MC to such scenario.

We also observe that the non-trivial $\theta_1$ parameter goes through large variations in all learning scenarios –see Figure 5(a)– which tend to happen when the value of $\theta_1$ goes below its optimum of 3. This is due to the fact that the non-deterministic probability at $K_1 - 1$ of the acceptance policy suddenly changes from being close to 0 to being close to 1, which makes the learning agent suddenly receive a signal that it should reject incoming jobs at a larger threshold value.

<footnote n="7">⁷When the stationary distribution of the states in $\partial A^c$ is unknown, one should use the same procedure described in Subsection 4.2.1 for the selection of the start state.</footnote>
**Algorithm 1** FVRL algorithm for finding a local optimizer $\theta^* \in \mathbb{R}^d$ of $g(\pi_\theta)$, the long-run expected reward defined in (6), using a parameterised policy $\pi_\theta$ for the MDP introduced in Section 2.

**Data:**

A. Simulator or emulator of the MDP.

B. Absorption set $\mathcal{A}$ required by the FVEE algorithm presented in Section 3.2.

C. $N_0, M_0$: respectively, the initial number of FV particles and the initial number of return cycles to $\mathcal{A}$, used by the FVEE algorithm to estimate the stationary probability distribution, $p^{\pi_\theta}$. We let $N_k$ and $M_k$ denote the value of the same concepts at policy iteration $k$.

D. $V, S_0$: respectively, the number of replications and the initial maximum number of time steps allowed to estimate the function $\nabla_\theta g_x(\pi_\theta)$ defined in (8), with $\theta \in \mathbb{R}^d$. We let $S_k$ denote the value of the same concept at policy iteration $k$.

E. $\mathcal{K}$: Compact set defining the support for $\theta$, e.g. a hyper-rectangle.

F. $\theta_0$: initial guess of parameter $\theta \in \mathcal{K}$.

G. $\epsilon > 0$: exponent defining the rates of increase of $N_k, M_k$ and $S_k$ in terms of policy iteration $k$.

H. $\alpha_0 > 0$: initial learning rate or step size, used in the $\theta$ update performed by the algorithm.

**Result:** An estimate of a local optimizer $\theta^* \in \mathbb{R}^d$ for the long-run expected reward $g(\pi_\theta)$.

Initialize policy iteration, $k \leftarrow 0$

Let $S'$ be the set of states where the policy gradient is not zero for at least one action, i.e. $S' = \{x \in S : \exists \text{ action } a \in \mathcal{U} \text{ s.t. } \nabla_{\theta_k} \pi_{\theta_k}(a|x) \neq 0\}$.

repeat

a) Simulate or emulate $V$ times $|\mathcal{U}|$ parallel copies of the $\{X_n\}_{n \in \mathbb{N}}$ process, each starting at $x \in S'$ and at a different $a \in \mathcal{U}$, and collect the rewards $\{R_{v,n}^{(x,a)}\}_{n \ge 1}$ observed by each $a$-th copy for replication $v$. Strongly couple any two chains after they meet and run the simulation until all $|\mathcal{U}|$ simulated chains meet or until $S_k$ steps are taken on each parallel chain.

b) Estimate the contribution to the gradient by each $x \in S'$:

$$ \hat{\nabla}_{\theta} g_x(\pi_{\theta_k}) \leftarrow \frac{1}{V} \sum_{v=1}^{V} \sum_{a \in \mathcal{U}} \sum_{n=1}^{\min(\tau_v, S_k)} R_{v,n}^{(x,a)} \nabla_{\pi_{\theta_k}} \pi_{\theta_k}(a|x), $$
where $\tau_v$ is the coupling time of all $|\mathcal{U}|$ chain copies described above for replication $v$.

c) Compute $\hat{p}_{\text{FV}}^{\pi_{\theta_k}}(x)$, the FV estimator of the stationary probability $p^{\pi_{\theta_k}}$ using $N_k$ particles and $M_k$ return cycles to $\mathcal{A}$, as described in Section 3.2.

d) $\hat{\nabla}_{\theta_k} g(\pi_{\theta_k}) \leftarrow \sum_{x \in S'} \hat{p}_{\text{FV}}^{\pi_{\theta_k}}(x) \hat{\nabla}_{\theta} g_x(\pi_{\theta_k})$

e) $\theta_{k+1} \leftarrow \mathcal{P}_{\mathcal{K}}[\theta_k + \alpha_k \hat{\nabla}_{\theta_k} g(\pi_{\theta_k})]$, where $\mathcal{P}_{\mathcal{K}}$ is a projection operator on compact set $\mathcal{K}$.

f) $k \leftarrow k + 1$

g) Increase $N_k, M_k$ and $S_k$ as $N_k = \lceil N_0 k^\epsilon \rceil, M_k = \lceil M_0 k^\epsilon \rceil, S_k = \lceil S_0 \log k^\epsilon \rceil$, where $\lceil \cdot \rceil$ is the ceiling operator.

h) Reduce the learning step by $\alpha_k = \alpha_0/k$

until done;
<figure><img src="image_7.png" /><figcaption>Figure 5: Comparison of FVRL learning (middle) with Monte-Carlo learning (right) and learning using the true stationary probabilities (left), which gives the scenario closest to the true learning process, where only the difference in the *Q* values is estimated. The plots show the results over 20 experiments run over 30 learning steps on a loss network with *R* = 6 servers, *λ* = [1, 5], *ρ* = [0.5, 0.3], with blocking costs *C* = [2 × 10³, 2 × 10⁴]. At each learning step, the FVRL setup sets the absorption set size by job class at *J*ᵢ/*K*ᵢ = [0.3, 0.5], the number of particles at *N* = 300 and the number of arrival steps at *T* = 500. In the MC learning case, each experiment uses as many number of events per learning step as the average number of events observed in FVRL by learning step. In all cases the initial guess is *θ* = [0.1, 0.1] (*K* = [2, 2]), the learning parameter starts at *α* = 1.0 and is then decreased inversely proportional to the learning step. The top row shows each of the 20 realized paths. The bottom row shows the distribution of the visited *θ* values by learning step, where the darker band includes 50% of the paths around their median and the lighter band includes 100% of the paths. The interesting learning paths to compare are those associated to *θ*₁, as *θ*₂ is learned equally well as 4.9 by all methods (which is shown as dotted lines in (a)).</figcaption></figure>

# 5 Conclusions and future work

The Fleming-Viot particle system was presented as an efficient alternative to Monte-Carlo for the exploration of environments where rewards are sparse and their occurrence is rare. Its application to the estimation of the blocking probability in an *M/M/1/K* queueing system and of the expected rejection

cost in a loss network serving different job classes served as test benches. In the *M/M/1/K* system, the method proved to be much more efficient than Monte-Carlo for large capacities *K*, where the latter completely fails, and in the loss network system it was able to estimate the expected cost more accurately and more precisely than Monte-Carlo.

Results on optimal control of the above systems using pol-
### 3.3.2 Convergence guarantees of policy gradient learning with FVRL

FVRL is a stochastic gradient ascent algorithm with a biased gradient estimate. In this context, Theorem 2.1 of [14, Section 5.2] can be used to provide a convergence guarantee to a local optimizer of $g(\pi_\theta)$. The convergence property and the conditions under which it is valid are stated in the following proposition.

**Proposition 4 (Convergence of FVRL to a local optimizer of $g(\pi_\theta)$).** Assume $\pi_\theta, \theta \in K \subset \mathbb{R}^d$ compact, is continuously differentiable. Assume moreover that the hyperparameters defined in the FVRL policy gradient Algorithm 1, scale as follows:

$$M_k \sim k^e, N_k \sim k^e, S_k \sim \log k^e, \text{ and } \alpha_k \sim 1/k.$$

Then the FVRL policy gradient algorithm converges to a local optimizer in $K$ of the long-run expected reward $g(\pi_\theta)$ defined in (6).

The proof is given in Appendix C and makes use of the following lemma bounding the finite sample bias $\beta = \mathbb{E}^{\pi_\theta}[\hat{\nabla}_{\theta}g(\pi_\theta)] - \nabla_{\theta}g(\pi_\theta)$ of the FVRL estimator of the gradient.

**Lemma 5 (Bias of the estimated gradient).** The bias $\beta$ of the FVRL estimator of the gradient of the long-run expected reward is bounded by

$$|\beta| \le A \left( \frac{1}{\sqrt{M}} + \frac{1}{\sqrt{N}} \right) + Be^{-cS} + Ce^{-cS} \left( \frac{1}{\sqrt{M}} + \frac{1}{\sqrt{N}} \right),$$

for some constants $A, B, C, c > 0$, where $M$ and $N$ are defined in Theorem 2, and $S$ is the maximum number of time steps allowed in the simulation used to estimate the function $\nabla_{\theta}g_x(\pi_\theta)$ defined in (8).

The proof is given in Appendix C.

We note that under a general parameterised policy setting, global convergence is only ensured under convexity assumptions. See for example [19] for the finite sample convergence analysis of stochastic gradient descent under a convex objective, and [1] for a detailed convergence analysis of policy gradient methods under the discounted reward criterion.

## 4 Application to stochastic networks with blocking

In this section we apply the methodologies outlined in Sections 3.2 and 3.3 to two different systems:

1. an $M/M/1$ queue, which is a unidimensional state system that will help us illustrate the basic concepts;

2. a loss network system with $R$ servers receiving two classes of jobs. which we will use to illustrate the application to a multidimensional state system.

In both cases, we assume there exists a small subset $C$ of the state space where costs may be generated due to the rejection of an incoming job (blocking), and that there is no reward for job acceptance. This makes the reward landscape sparse. For each system we will consider two goals:

a) the accurate estimation of the long-run expected cost, $\mathbb{E}^\pi(\eta)$, defined in (1), where $\eta(x)$ is the cost of rejecting an incoming job when the system is at state $x$;

b) the learning of the blocking sizes that minimize $\mathbb{E}^\pi(\eta)$.
icy gradient were also presented. The proposed FVRL algo-
rithm tends to find the optimum parameters significantly faster
than Monte-Carlo learning. In this case, we note that the ac-
curacy of the stationary probability estimator is not as crucial
as in the estimation problem, because the algorithm is able to
learn as long as it receives a signal from the rare states.

In future work, we intend to extend the FV and the FVRL
algorithms to environments other than queues and networks,
to more traditional RL environments such as labyrinths or the
mountain car, where it will be crucial to define the absorption
set A adaptively, i.e. based on the discovery of the states that
give no rewards during exploration of the environment.

Acknowledgements

This work was partially supported by the French National Re-
search Agency through projects ANR-11-LABX-0040 of the
"Investments for the Future" program and ANR-22-CE25-0013-
02 (EPLER), and by project LAGOON of the Stic-amsud pro-
gram.

References

[1] S. Asmussen. *Applied Probability and Queues*. Applications of mathematics : stochastic modelling and applied probability. Springer, 2003.

[2] A. Asselah, P.A. Ferrari, and P. Groisman. Quasi-stationary distributions and Fleming-Viot processes in finite spaces. *J. Appl. Probab.*, 48(2):322–332, 2011.

[3] Thomas Bonald, Matthieu Jonckheere, and Alexandre Proutière. Insensitive load balancing. *ACM Sigmetrics Performance Evaluation Review*, 32(1):367–377, 2004.

[4] Léon Bottou, Frank E. Curtis, and Jorge Nocedal. Optimization methods for large-scale machine learning. *SIAM Review*, 60(2):223–311, jan 2018.

[5] Pierre Brémaud. Markov chains: Gibbs fields, Monte Carlo simulation, and queues, volume 31. Springer Science & Business Media, 2013.

[6] Yuri Burda, Harri Edwards, Deepak Pathak, Amos Storkey, Trevor Darrell, and Alexei A Efros. Large-scale study of curiosity-driven learning. arXiv preprint arXiv:1808.04355, 2018.

[7] K. Burdzy, R. Holyst, and P. March. A Fleming-Viot particle representation of the Dirichlet Laplacian. *Comm. Math. Phys.*, 214(3):679–703, 2000.

[8] Bertrand Cloez and Josué Corujo. Uniform in time propagation of chaos for a moran model. arXiv preprint arXiv:2107.10794, 2021.

[9] Bertrand Cloez and Marie-Noémie Thai. Quantitative results for the fleming-viot particle system and quasi-stationary distributions in discrete space. *Stochastic Processes and their Applications*, 126(3):680–702, 2016.

[10] J. N. Darroch and E. Seneta. On quasi-stationary distributions in absorbing continuous-time finite Markov chains. *J. Appl. Probability*, 4:192–196, 1967.

[11] Gabriel Dulac-Arnold, Daniel J. Mankowitz, and Todd Hester. Challenges of real-world reinforcement learning. CoRR, abs/1904.12901, 2019.

[12] P. Groisman and M. Jonckheere. Simulation of quasi-stationary distributions on countable spaces. Markov Process. Related Fields, 19(3):521–542, 2013.

[13] G. Yin H. J. Kushner. Stochastic Approximation and Recursive Algorithms and Applications. Springer-Verlag, 2003.

[14] Ger Koole. Structural results for the control of queueing systems using event-based dynamic programming. *Queueing systems*, 30(3):323–339, 1998.

[15] Ger Koole. Monotonicity in Markov reward and decision chains: Theory and applications, volume 1. Now Publishers Inc, 2007.

[16] Antonio Massaro, Francesco De Pellegrini, and Lorenzo Maggi. Optimal trunk-reservation by policy learning. In *IEEE INFOCOM 2019*, apr 2019.

[17] Maja J Mataric. Reward functions for accelerated learning. In *Machine learning proceedings 1994*, pages 181–189. Elsevier, 1994.

[18] Deepak Pathak, Pulkit Agrawal, Alexei A. Efros, and Trevor Darrell. Curiosity-driven exploration by self-supervised prediction. In Doina Precup and Yee Whye Teh, editors, *Proceedings of the 34th International Conference on Machine Learning*, volume 70 of *Proceedings of Machine Learning Research*, pages 2778–2787. PMLR, 06–11 Aug 2017.

[19] Alexey Piunovskiy and Yi Zhang. Continuous-time markov decision processes. Probability Theory and Stochastic Modelling, 2020.

[20] M. L. Puterman. Markov Decision Processes: Discrete Stochastic Dynamic Programming. John Wiley & Sons, 2005.

[21] Keith W. Ross. Multiservice Loss Models for Broadband Telecommunication Networks. Springer, 1995.

[22] Keith W. Ross and Danny H.K. Tsang. The stochastic knapsack problem. IEEE Transactions on Communications, 37(7):740–747, July 1989.

[23] David Silver, Julian Schrittwieser, Karen Simonyan, Ioannis Antonoglou, Aja Huang, Arthur Guez, Thomas Hubert, Lucas Baker, Matthew Lai, Adrian Bolton, et al. Mastering the game of go without human knowledge. *nature*, 550(7676):354–359, 2017.

[24] Richard S Sutton, David A McAllester, Satinder P Singh, and Yishay Mansour. Policy gradient methods for reinforcement learning with function approximation. In *Advances in neural information processing systems*, pages 1057–1063, 2000.

[25] Martin J Wainwright. *High-dimensional statistics: A non-asymptotic viewpoint*, volume 48. Cambridge University Press, 2019.
To simplify the exposition and to be able to do theoretical computations against which results are compared, we fix the parameters of the underlying Markov process to known values. Nevertheless, the method is designed to be applied in practice to systems with possibly unknown parameters, the only condition being that the system can be simulated or emulated, i.e. that its next state and the reward received are known once an action is applied at any given system state.

The algorithms for the FVEE estimation of the expected cost and of its minimisation using FVRL, applied to both the $M/M/1$ queue system and the loss network, are presented in Algorithms 2 and 3 in Appendix E, respectively, and implemented in the following repository: https://github.com/mastropi/RL, which includes a few details about the estimation and learning process. We note that to keep implementation simple the values of hyperparameters $M, N$, and $S$ are not scaled by the policy iteration step as required by Proposition 4.

In each application problem, the method's performance is compared with vanilla Monte-Carlo (MC) which is used as benchmark. To ensure a fair comparison, the MC estimation of expectations is based on the same number of events as the FVEE estimation, and the MC learning is based on the same number of events per policy iteration as the average number of events observed in FVRL across the whole policy learning process. The details of the estimation and learning processes by MC are further described in each subsection below.

## 4.1 M/M/1 system

We consider the MDP representing the dynamics of the well-known $M/M/1$ queue system under threshold-type admission control policies. Following the terminology defined in Section 2, for some $L \in \mathbb{N}$, the MDP can be described as: ($S = \{0, 1, \dots, L\}, U = \{0, 1\}, q = \{\lambda, \mu\}, \mathcal{R} = r(x, a) = \mathbb{1}_{\{x=K, a=0\}}, \bar{\Pi}_L$) where the Markov process $X_t^\pi$ is the length of the queue having job arrival rate $\lambda$ and service rate $\mu$, actions 0 and 1 respectively indicate "reject" and "accept" of an incoming job, and $\bar{\Pi}_L$ is the class of admission control policies that reject an incoming job at just one state, when $X_t^\pi = K$ for some $0 \le K \le L$.

All presented results use $\lambda = 0.7$ and $\mu = 1.0$.

### 4.1.1 Estimation of the expected cost using FVEE

We apply the methodology outlined in Section 3.2 to efficiently estimate the long-run expected cost (thereafter "expected cost") when blocking is a rare event. Since the rejection cost is 1, estimating the expected cost is tantamount to estimating the blocking probability, which can be quite small depending on the $M/M/1$ system parameters and admission control policy (see numerical examples in the Introduction).

In order to estimate the blocking probability using Fleming-Viot, we define the function $\eta$ introduced in Section 2 as $\eta^{\pi\theta}(x) = \mathbb{1}_{\{x=K\}}$, which is sparse and non-zero in the single-state set $C = \{K\}$. As a consequence and using the PASTA property of Poisson arrivals, the expectation $\mathbb{E}^\pi(\eta)$ in (1) becomes the blocking probability, which can be written as

$$p^{\pi}(K) = \frac{\int_{0}^{\infty} \phi_{t}^{J}(K) \mathbb{P}_{J}(T_{K} > t) dt}{\mathbb{E}_{J-1}T_{A}}, \qquad (10)$$

where we have used that any set $A = \{0, 1, \dots, J-1\}$ with $J \le K$ is a valid absorption set, making $\bar{\partial}A$ and $\bar{\partial}A^c$ two single-state entrance boundary sets equal to $\{J-1\}$ and $\{J\}$, respectively, which makes it possible to simplify (3) into (10).

According to Theorem 2, the FV estimator of $p^\pi(K)$ converges to its true value as both the number $N$ of particles of the FV system, and the number $M$ of return cycles to $A$ increase. While $N$ is a hyperparameter that is easy to define, as it is a fixed number, $M$ is not as straightforward, as it is a random number (it depends on the number of return cycles to $A$ observed during the Markov chain excursion). We will therefore parameterise the Markov chain simulation for the estimation of the denominator, $\mathbb{E}_{J-1}T_A$, by the number of arrival events $T$ that must be observed
# A Estimation of $\mathbb{E}^{\pi}(\eta)$ using Fleming-Viot particle systems

An estimator of $\mathbb{E}^{\pi}(\eta)$ using (3) is constructed from estimators $\hat{f}_{\eta}$, $\hat{g}$ of functions $f_{\eta}$, $g$ as

$$ \hat{\mathbb{E}}^{\pi}(\eta) = \int_{0}^{\infty} \hat{f}_{\eta}(t) \hat{g}(t) dt. \qquad (12) $$

Firstly, we explain the construction of the estimator of $g$, which is a ratio of the quantities $\mathbb{P}_{\vec{\partial}A^c}(T_K > t)$ and $\mathbb{E}_{\vec{\partial}A}(T_A)$, each of which is estimated from two separate simulations, as follows: (i) a simulation of the Markov process starting at an arbitrary state $x \in \vec{\partial}A$, and (ii) a simulation of a predefined number of $N$ independent copies of the Markov process starting at a randomly selected state $x \in \vec{\partial}A^c$ following the entrance state distribution into $A^c$ under stationarity, estimated from the first simulation. As explained below, simulations (i) and (ii) contribute to the estimation of $\mathbb{E}_{\vec{\partial}A}(T_A)$ while simulation (ii) contributes to the estimation of $\mathbb{P}_{\vec{\partial}A^c}(T_K > t)$.

The details are as follows: we set $\tau_{A,0} \doteq 0$ and define a sequence of stopping times $\tau_{A^c,i}$, $\tau_{A,i}$, associated to the events of entry, in the $i$-th cycle, into $A^c$ and into $A$, respectively, as follows:

$$ \tau_{A^c, i} = \inf_{t > \tau_{A, i-1}} \{X_t^\pi \in A^c\}, $$

$$ \tau_{A,i} = \inf_{t > \tau_{A^c,i}} \{X_t^\pi \in A\}, $$

for $i \ge 1$. We also define $T_{E,i} = \tau_{A^c,i} - \tau_{A,i-1}$, the entry time into $A^c$ within cycle $i$, and $T_{K,i} = \tau_{A,i} - \tau_{A^c,i}$, the killing time within cycle $i$. We note that $T_{A,i} = T_{E,i} + T_{K,i}$, which will be used below to construct the estimator of the denominator of $g$, $\mathbb{E}_{\vec{\partial}A}(T_A)$.

The first simulation consists of running the process $X_t^\pi$ until a predefined number $M_0 + M$ entry times $\{T_{E,i}\}_{i=1}^{M_0+M}$ are observed, where we consider the first $M_0$ observations to be burn-in in order to assume stationarity of the process thereafter. From this simulation we compute the empirical entrance state distribution into $A^c$, which can be considered as an estimator of $p_{\vec{\partial}A^c}^{\pi}$, as $\hat{p}_{\vec{\partial}A^c}^{\pi} = \frac{1}{M} \sum_{i=M_0+1}^{M_0+M} \mathbf{1}_{X_{\tau_{A^c,i}}^{\pi}}$.

The second simulation of the $N$ independent copies of $X_t^\pi$ is started at states $x_i \in \vec{\partial}A^c$ randomly chosen according to the estimated stationary distribution $\hat{p}_{\vec{\partial}A^c}^{\pi}$ (this is trivial when $\vec{\partial}A^c$ has a single state), yielding $N$ killing times $\{T_{K,i}\}_{i=1}^N$.

From the last $M$ entry times and the $N$ killing times, we define a Monte-Carlo estimator of $g$ as the ratio

$$ \hat{g}(t) = \frac{\frac{1}{N} \sum_{i=1}^{N} \mathbf{1}_{T_{K,i} > t}}{\frac{1}{M} \sum_{i=M_0+1}^{M_0+M} T_{E,i} + \frac{1}{N} \sum_{i=1}^{N} T_{K,i}}. \qquad (13) $$

Next, we explain the construction of the estimator of $f$, which is a $\eta$-weighted sum of $\phi_t^{\vec{\partial}A^c}$ over all states in $A^c$, using the Fleming-Viot $N$-particle system driven by $X_t^\pi$. The Fleming-Viot system, denoted by $(\xi_t^\nu)_{t \ge 0}$, is simulated as described in Section 3.2, using $A$ as the absorption set and starting each particle $i$ at $x_i$ a randomly chosen state according to the distribution $\nu(x) = p_{\vec{\partial}A^c}^{\pi}(x) \mathbf{1}_{x \in \vec{\partial}A^c}$. That is, all particles start at the boundary of $A^c$, as required by the quantity to estimate $\phi_t^{\vec{\partial}A^c}$. We let $m(\cdot, \xi) : A^c \to [0, 1]$ denote the empirical distribution of the $N$ particles with positions described by vector $\xi$, defined as the empirical mean $m(x, \xi) \doteq \frac{1}{N} \sum_{i=1}^N \mathbf{1}_{\xi(i)=x}$, $\forall x \in A^c$. Since $m(\cdot, \xi^\nu)$ is an estimator of $\phi_t^{\vec{\partial}A^c}$ (because $\nu$ is restricted to the boundary of $A^c$ where it is equal to the entrance state distribution into $A^c$ under stationarity), $f$ can be estimated by

$$ \hat{f}_{\eta}(t) = \sum_{x \in A^c} \eta(x) m(x, \xi_t^{\nu}). \qquad (14) $$

We note that, by construction, $\hat{g}(t) = 0$ for $t > T_{K,\max} = \max\{T_{K,i} : 1 \le i \le N\}$. Therefore, since we wish to compute $\hat{\mathbb{E}}^{\pi}(\eta) = \int_{0}^{\infty} \hat{f}_{\eta}(t) \hat{g}(t) dt$, we only need to simulate the Fleming-Viot process until time $T_{K,\max}$ is reached. Also, since both $\hat{f}_{\eta}(t)$ and $\hat{g}(t)$ are almost surely piecewise constant functions and $\hat{g}(t) = 0$ for $t > T_{K,\max}$, the integral $\int_{0}^{\infty} \hat{f}_{\eta}(t) \hat{g}(t) dt$ is a finite sum that can be easily computed.
before ending the simulation, which we will refer to as “the number of arrival events for estimating $E(T_A)$” thereafter
(we omit the subscript of the expectation for conciseness).

The three quantities in (10) contributing to the blocking probability are estimated following the steps described in Algorithm 2 in Appendix E, which implements the estimation methodology described in Appendix A.

**Remark 3.** There is a trade-off in the choice of $J$, the state that defines the size of the absorption set $\mathcal{A}$: for smaller $J$, the return times to $\mathcal{A}$ will be smaller, requiring a smaller number $T$ of arrival events for estimating $E(T_A)$ in (10); at the same time, however, visiting the rare blocking state $K$ becomes rarer, requiring a larger number of particles $N$ for an accurate estimation of the numerator in (10). The opposite is true for larger values of $J$. A detailed analysis on this trade-off, in terms of the impact on estimation accuracy of different $N$ and $T$ choices given the value of $J$, is presented in Appendix D.

We now study the convergence of FVEE as either $N$ or $T$ increases, and compare it with the benchmark MC estimator, obtained from a direct application of expression (2), i.e. as the fraction of the time spent at state $K$ and the total time of return cycles to the initial state $x = J - 1$ observed during the simulation. To guarantee a fair comparison between the two methods, we start the simulation at $x = J - 1$, so that both methods start at the same "distance" from the blocking state $K$, and let the simulation run until the same number of events observed by FVEE is reached.

Figure 2 shows violin plots of FVEE and of the MC estimator, as $N$ increases on the left column plots, and as $T$ increases on the right column plots. Their values were chosen as described in Appendix D. We considered the cases $K = 20$ (top row) and $K = 40$ (bottom row), which are regarded to represent moderate and large capacities based on their blocking probabilities at the considered value for $\rho = 0.7$ of order $10^{-4}$ and $10^{-7}$, respectively. The size $J$ of the absorption set $\mathcal{A}$ is held fixed at $J = 12$, which corresponds to choosing the states with stationary probability larger than 0.5%, for both $K = 20$ and $K = 40$.

We observe the following in terms of convergence of FVEE and of the MC estimators to the true blocking probability, which in the plots is represented by a horizontal gray line:

1. Both FVEE and MC converge to the true blocking probability when $K = 20$ (Figures 2(a,b)). MC presents a smaller variability than FVEE in the convergence analysis with $N$ (Figure 2(a)) but a larger variability than FV in the convergence analysis with $T$ (Figure 2(b)). This is due to the fact that Figure 2(a) is obtained from simulations whose events largely outnumber the events in Figure 2(b) by as much as 10 times ($\sim 20,000$ vs. $\sim 3,000$), which allows MC to observe the blocking event at $x = K$ frequently enough for an accurate estimation. On the other hand, Figure 2(b) tells us that a much smaller number of events ($\sim 3,000$) is sufficient for FVEE to estimate the blocking probability accurately enough, but is not sufficient for the MC estimator. This demonstrates the higher efficiency of FV than MC in discovering the rare event at $K$.

2. When $K = 40$ (Figures 2(c,d)) the MC estimator basically fails as it almost never observes the blocking state $K$.

3. Increasing $N$ has a larger impact on increasing FVEE's computational complexity than increasing $T$, but with no particular gain on the estimation accuracy or variability. This conclusion is obtained by observing that the number of events at the rightmost violin plot in Figure 2(a) ($\sim 20,000$ events), where $N \sim 4,000$ and $T \sim 4,000$, is about 10 times larger than the number of events at the rightmost violin plot in Figure 2(b) ($\sim 3,000$ events), where $N = 30$ and $T \sim 4,000$, while the distribution of the estimated probabilities is almost the same. That is, a 10-fold increase of the number of particles $N$ increased ten times the total number of observed events but did not increase the estimation accuracy nor decrease its variability. On the other hand, a mere 4-fold increase in $T$ does not impact the order of magnitude of the total number of observed events (see panels (b) and (d)) but decreases the variability of the blocking probability estimate (to about half in the case of panel (d)). This confirms the conclusion of the analysis presented in Appendix D that a larger value of $T$ is more important for a quality estimation of the blocking probability than a larger value of $N$ in the $M/M/1/K$ queue system.
B Proof of Theorem 2: Error of the Fleming-Viot estimator

Using the notation $f_\eta, g, \hat{f}_\eta, \hat{g}$ introduced in (3) and (14), (13), we have $\mathbb{E}^\pi(\eta) = \int_0^\infty f_\eta(t)g(t)dt$ and $\hat{\mathbb{E}}^\pi(\eta) = \int_0^\infty \hat{f}_\eta(t)\hat{g}(t)dt$.
We are thus interested in bounding:

$$ \mathbb{E} \left| \int_{0}^{\infty} \hat{f}_{\eta} \hat{g} dt - \int_{0}^{\infty} f_{\eta} g dt \right| $$

We start by decomposing the problem of upper bounding the above quantity into two subproblems in the following way:

$$ \mathbb{E} \left| \int_{0}^{\infty} \hat{f}_{\eta} \hat{g} dt - \int_{0}^{\infty} f_{\eta} g dt \right| \leq \mathbb{E} \left| \int_{0}^{\infty} (f_{\eta} - \hat{f}_{\eta}) g dt \right| + \mathbb{E} \left| \int_{0}^{\infty} \hat{f}_{\eta} (\hat{g} - g) dt \right| \quad (15) $$

We start by bounding the first term on the right hand side. For this purpose we will need the uniform propagation of chaos bound presented in (4), that is:

$$ \sup_{\|\phi\|_{\infty} \le 1} \sup_{t \ge 0} \mathbb{E} \left| [m(\cdot, \xi_t^{\nu})(\phi)] - \phi_t^{\nu}(\phi) \right| \le \frac{C_{\text{FV}}}{\sqrt{N}}, $$

from which it follows that:

$$ \sup_{t \ge 0} \mathbb{E} \left| \hat{f}_{\eta}(t) - f_{\eta}(t) \right| \le \frac{C_{\text{FV}}}{\sqrt{N}}. \quad (16) $$

As was mentioned in subsection 3.3, this bound follows directly from [8, Theorem 1.4]. The assumptions of [8, Theorem 1.4] have a very general form, but it is easy to check that they are trivially satisfied in our simple case. The assumption (I) on initialization is satisfied by our assumption that the FV particle system is started at the position of $N$ i.i.d samples from $p_{\bar{\partial}_A^c}^\pi$. Assumption (C1) has several parts: the uniform bound on selection rates (which are, in our case, the rates of jumps out of $\mathcal{A}^c$), follows from the fact that the state space is finite); the rest of assumption (C1) is trivially satisfied when we take $V_\mu^d(x)$ to be the rate of jump out of $\mathcal{A}^c$ from the state $x \in \mathcal{A}^c$ for any $\mu$, and set function $V_\mu^b(y)$, $V_\mu^s(x, y)$ equal to zero. Finally, assumption (C2) follows from the fact that we are working with an irreducible Markov chain on a finite state space. Therefore, using the triangle inequality and the inequality (16), we obtain:

$$ \begin{aligned} \mathbb{E} \left| \int_{0}^{\infty} (f_{\eta} - \hat{f}_{\eta}) g dt \right| &\leq \mathbb{E} \int_{0}^{\infty} \left| (f_{\eta} - \hat{f}_{\eta}) \right| g dt \\ &\leq \int_{0}^{\infty} \mathbb{E} \left| (f_{\eta} - \hat{f}_{\eta}) \right| g dt \\ &\leq \frac{C_{\text{FV}}}{\sqrt{N}} \int_{0}^{\infty} g dt = \frac{\mathbb{E}_{\bar{\partial}_A^c} T_K}{\mathbb{E}_{\bar{\partial}_A} (T_A)} \frac{C_{\text{FV}}}{\sqrt{N}}, \end{aligned} $$

where in the last line we also use the 'wedding cake decomposition', $\int_0^\infty \mathbb{P}_{\bar{\partial}_A^c}(T_K > t)dt = \mathbb{E}_{\bar{\partial}_A^c} T_K$.

Since $\eta(x)$ is bounded, without loss of generality we can consider that $\sup_{x \in \mathcal{A}^c} |\eta(x)| \le 1$, and use that $|\hat{f}_\eta(t)| \le 1$ for $t \ge 0$ to get:

$$ \left| \int_{0}^{\infty} \hat{f}_{\eta}(\hat{g} - g) dt \right| \le \int_{0}^{\infty} |\hat{g} - g| dt. $$

We thus wish to estimate $\mathbb{E} \int_0^\infty |\hat{g} - g| dt$. For convenience, we define a random variable $\bar{T}_A$ with the distribution of $T_A$ when $X_t^\pi$ is started with distribution $p_{\bar{\partial}_A}^\pi$, and a random variable $\bar{T}_K$ with the distribution of $T_K$ when $X_t^\pi$ is started with distribution $p_{\bar{\partial}_A^c}^\pi$. Since we start the simulation of $X_t^\pi$ for the purpose of estimating $g$ with distribution $p_{\bar{\partial}_A}^\pi$, we do not need any burn-in. We therefore take $M_0 = 0$. We also note that since we start the simulation at the distribution $p_{\bar{\partial}_A}^\pi$, it follows from renewal theory [1] that the inter-arrival times $\mathcal{T}_{A,i}$ used to construct the estimator $\hat{g}$ are i.i.d. with distribution $\bar{T}_A$.

We also introduce additional shorthand notation for the numerators and denominators of $g(t)$ and $\hat{g}(t)$. We denote $N_t = \mathbb{P}_{\bar{\partial}_A^c}(T_K > t)$ and $D_A = \mathbb{E}_{\bar{\partial}_A} T_A$, $D_K = \mathbb{E}_{\bar{\partial}_A^c} T_K$. We also denote by $\tilde{N}_t$, $\tilde{D}_A$ the estimators of $N_t$, $D_A$, that is $\tilde{N}_t = \frac{1}{N} \sum_{i=1}^N 1_{T_K, i} > t$ and $\tilde{D}_A = \tilde{D}_E + \tilde{D}_K$, with

$$ \hat{D}_E = \frac{1}{M} \sum_{i=1}^{M} T_{E,i}, $$

$$ \hat{D}_K = \frac{1}{N} \sum_{i=1}^{N} T_{K,i}, $$

We thus have $g(t) = \frac{N_t}{D}$ and $\hat{g}_t = \frac{\tilde{N}_t}{D_A}$.

We are interested in bounding:

$$ \mathbb{E} \int_{0}^{\infty} \left| \frac{\tilde{N}_{t}}{\tilde{D}_{A}} - \frac{N_{t}}{D_{A}} \right| $$
<figure><img src="image_4.png" /><figcaption>Figure 2: Violin plots showing the Fleming-Viot (left, green) (ours) and Monte-Carlo (right, red) estimations of the blocking probability for an M/M/1 queue system with capacity K and $\lambda = \rho = 0.7$, for two different capacities K = 20, 40, and different hyperparameter settings N and T of the Fleming-Viot expectation estimator (FVEE). In total, 8 replications were run for each hyperparameter setting, but not all of them yielded an estimate of the blocking probability by FVEE due to a requirement of observing at least 5 cycles in the estimation of $E(T_A)$. The left plots (a) and (c) analyze the convergence characteristics of the estimator as the number of particles N increases with T fixed at the indicated value. The right plots (b) and (d) analyze the convergence characteristics of the estimator as the number T of arrival events for estimating $E(T_A)$ increases with N fixed at the indicated value. The absorption set size is $J = 12$ in all cases, whose value was chosen as described in the text. The interior colored horizontal line in each violin plot represents the median. In the corresponding MC experiments, since neither N nor T are defined as hyperparameters of the simulation, they cannot be used as horizontal axis. In their place, the average of the number M of observed return cycles to the initial state $J-1$, across replications, is shown at the bottom axis. In both the FVEE and the MC plots, the top horizontal axis reports the average number of events observed across replications, which coincide in FVEE and MC by experiment design (see Section 4).</figcaption></figure>

### 4.1.2 Learning an optimal blocking size K with FVRL

We illustrate the FVRL algorithm for learning an optimal blocking size K by defining a cost function of rejecting an incoming job that is exponentially increasing with the queue size at the time of the job arrival. Note that the cost needs to be exponentially increasing with the queue size in order to obtain an optimisation problem with a non-trivial optimum (i.e. where the optimum K is finite). This is due to the fact that the stationary probability of the M/M/1 queue is exponentially decreasing with the queue size when $\rho < 1$.

More precisely, given the queue size x, the cost function is defined as $r(x, a) = B(1 + b^{x-x_{ref}})1_{\{a=0\}}$, where B, b
Using the triangle inequality $|\frac{\hat{N}_t}{\hat{D}_A} - \frac{N_t}{D_A}| \le |\frac{\hat{N}_t}{\hat{D}_A} - \frac{N_t}{D_A}| + |\frac{\hat{N}_t}{\hat{D}_A} - \frac{N_t}{D_A}|$ we get:

$$ \mathbb{E} \int_0^\infty \left| \frac{\hat{N}_t}{\hat{D}_A} - \frac{N_t}{D_A} \right| dt \le \mathbb{E} \int_0^\infty \hat{N}_t \left| \frac{1}{\hat{D}_A} - \frac{1}{D_A} \right| dt + \mathbb{E} \int_0^\infty \frac{1}{D_A} \left| \hat{N}_t - N_t \right| dt, $$

Using the formula $\int_0^\infty \hat{N}_t dt = \frac{1}{N} \sum_{i=1}^N T_{\kappa,i} = \hat{D}_\kappa$, the first term on the right hand side of the above bound is equal to $\mathbb{E} |\frac{\hat{D}_\kappa}{\hat{D}_A} - \frac{\hat{D}_\kappa}{D_A}|$. To bound this quantity, we introduce an event $B = \{\hat{D}_A < \frac{1}{2} D_A\}$. We use the decomposition

$$ \mathbb{E} \left| \frac{\hat{D}_\kappa}{\hat{D}_A} - \frac{\hat{D}_\kappa}{D_A} \right| = \mathbb{E} \mathbf{1}_B \left| \frac{\hat{D}_\kappa}{\hat{D}_A} - \frac{\hat{D}_\kappa}{D_A} \right| + \mathbb{E} \mathbf{1}_{B^c} \left| \frac{\hat{D}_\kappa}{\hat{D}_A} - \frac{\hat{D}_\kappa}{D_A} \right| \quad (17) $$

and bound each of the terms separately.

Since we always have $0 \le \frac{\hat{D}_\kappa}{\hat{D}_A} \le 1$ and on the set $B$ we have $0 \le \frac{\hat{D}_\kappa}{D_A} \le \frac{1}{2}$, we have:

$$ \mathbb{E} \mathbf{1}_B \left| \frac{\hat{D}_\kappa}{\hat{D}_A} - \frac{\hat{D}_\kappa}{D_A} \right| \le \mathbb{P}(B). $$

Since the Markov Process $X_t^\pi$ is irreducible and the state space $S$ is finite, it is geometrically ergodic [5]. It follows then that there exist constants $C, \lambda > 0$, such that $\mathbb{P}(\bar{T}_A > t) \le C \exp(-\lambda t)$. Therefore, by [25, Theorem 2.13], the random variable $\bar{T}_A$ is subexponential. Furthermore $\mathbb{P}(B) \le \mathbb{P}_{\bar{\delta}_A^c}(|\hat{D}_A - D_A| \ge \frac{1}{2} D_A)$. From the concentration bound for the standard estimator of the mean of subexponential variables [25][Equation 2.18], it follows that there exists $c > 0$:

$$ \mathbb{P}(B) \le e^{-c\sqrt{1/N+1/M}}. $$

To bound the second term in (17), we observe that the function $h(x) = 1/x$ is Lipschitz continuous on $[a, \infty)$ for any $a > 0$, with the Lipschitz constant $L_a = \sup_{x \in [a, \infty)} |h'(x)| = \frac{1}{a^2}$. Using this fact with $a = D_A/2$, we have:

$$ \mathbb{E} \mathbf{1}_{B^c} \hat{D}_\kappa \left| \frac{1}{\hat{D}_A} - \frac{1}{D_A} \right| \le \frac{4}{D_A^2} \mathbb{E} \hat{D}_\kappa \left| \hat{D}_A - D_A \right|. $$

Using the Cauchy-Schwartz inequality, we get:

$$ \begin{aligned} \mathbb{E} \hat{D}_\kappa |\hat{D}_A - D| &\le (\mathbb{E} \hat{D}_\kappa^2)^{1/2} (\mathbb{E} |\hat{D}_A - D_A|^2)^{1/2} \\ &\le \left( (\mathbb{E} \hat{D}_\kappa)^2 + \mathrm{Var}(\hat{D}_\kappa) \right)^{1/2} \left( \mathrm{Var}(\hat{D}_A) \right)^{1/2} \\ &\le \left[ \frac{(\mathrm{Var}_{\bar{\delta}_A} T_E)^{1/2}}{\sqrt{M}} + \frac{(\mathrm{Var}_{\bar{\delta}_A} T_K)^{1/2}}{\sqrt{N}} \right] \left( (D_\kappa)^2 + \frac{1}{M} \mathrm{Var}_{\bar{\delta}_A^c} (D_\kappa) \right)^{1/2} \end{aligned} $$

We therefore obtain

$$ \mathbb{E} \int_0^\infty \hat{N}_t \left| \frac{1}{\hat{D}_A} - \frac{1}{D_A} \right| dt \le \mathbb{E}_{\bar{\delta}_A^c} T_K \left[ \frac{(\mathrm{Var}_{\bar{\delta}_A} T_E)^{1/2}}{\sqrt{M}} + \frac{(\mathrm{Var}_{\bar{\delta}_A} T_K)^{1/2}}{\sqrt{N}} \right] + \mathcal{O}(\frac{1}{M}). $$

We are left with bounding

$$ \frac{1}{D_A} \mathbb{E} \int_0^\infty |\hat{N}_t - N_t| dt. $$

Since $\hat{N}_t$ is an average of $M$ Bernoulli random variables with mean $N_t$, we have:

$$ \begin{aligned} \frac{1}{\hat{D}_A} \mathbb{E} \int_0^\infty |\hat{N}_t - N_t| dt &= \frac{1}{\hat{D}_A} \int_0^\infty \mathbb{E} |\hat{N}_t - N_t| dt \\ &\le \frac{1}{\hat{D}_A} \int_0^\infty \sqrt{\mathbb{E} |\hat{N}_t - N_t|^2} dt \\ &= \frac{1}{\sqrt{N} D_A} \mathbb{E} \int_0^\infty \sqrt{N_t(1 - N_t)}, \end{aligned} $$
and $x_{\text{ref}}$ are positive constants. As described in Section 4.1, the MDP is driven by a threshold-type policy $\bar{\pi} \in \bar{\Pi}_L$, where $\bar{\pi}$ rejects an incoming job when $x = K$ for some $0 \le K \le L$. Under policy $\bar{\pi}$, the long-run expected cost to minimize becomes:

$$g(\bar{\pi}) = p^{\bar{\pi}}(K)r(K, 0) = \frac{(1-\rho)\rho^K}{1-\rho^{K+1}} B(1+b^{K-x_{\text{ref}}}). \quad (11)$$

Given the queue load $\rho < 1$, we choose $b > 1/\rho^{1+\rho}$ so that this objective function is a convex function of $K$ with a non-trivial (i.e. finite) minimum$^{(3)}$. The value of $x_{\text{ref}}$ is the reference queue size that is instrumental in defining the optimum threshold, $K^*$ (see Figure 3).

We now use the FVRL Algorithm 1 presented in Section 3.3.1 to learn the optimum threshold $K^*$. Following [17], we propose a parameterised acceptance policy $\pi(a = 1|x)$ that is a linear step function of the state $x$, which is deterministic for $x$ outside the interval $(\theta, \theta+1)$ and decreases linearly from 1 to 0 within the interval. That is, the acceptance policy parameterised by the positive real-valued $\theta$, is defined as:

$$\pi_{\theta}(a = 1|x) = \begin{cases} 1 & \text{if } x \le \theta, \\ 1 - (x - \theta) & \text{if } \theta < x < \theta + 1, \\ 0 & \text{if } x \ge \theta + 1. \end{cases} \quad (12)$$

Note that the policy is fully deterministic for integer-valued $\theta$, in which case the blocking size is $K = \theta + 1$. Otherwise, $K$ is defined as $\lceil \theta \rceil + 1$.

Using (7), the gradient of $g(\pi_{\theta})$ becomes

$$\frac{\partial g(\pi_{\theta})}{\partial \theta} = p^{\pi_{\theta}}(K-1)[Q^{\pi_{\theta}}(K-1, 1) - Q^{\pi_{\theta}}(K-1, 0)], \quad (13)$$

where $K-1$ is the smallest integer that is larger than or equal to $\theta$. We observe that this parameterisation leads, as expected, to gradients being 0 for $x \ne K-1$, that is, the gradient function is sparse in $x$. We note also that the gradient is discontinuous at $\theta$ and $\theta+1$, making the assumptions of Proposition 4 not fully satisfied. However, these two points have measure zero and therefore, with probability 1, no discontinuity is observed.

Under the piecewise-linear parameterised policy proposed in (12), the long-run expected cost becomes:

$$\begin{aligned} g(\pi_{\theta}) &= p^{\pi_{\theta}}(K-1)r(K-1, 0)(K-1-\theta) + p^{\pi_{\theta}}(K)r(K, 0) \\ &= p^{\pi_{\theta}}(\lceil \theta \rceil)r(\lceil \theta \rceil, 0)(\lceil \theta \rceil - \theta) + p^{\pi_{\theta}}(\lceil \theta \rceil + 1)r(\lceil \theta \rceil + 1, 0), \end{aligned} \quad (14)$$

where we have used that $K = \lceil \theta \rceil + 1$, with $\lceil \cdot \rceil$ the ceiling operator.

Details of the learning algorithm are given in Algorithm 3 in Appendix E.

To illustrate the FVRL algorithm, we consider an MDP with the following characteristics: the system load is $\rho = 0.7$, the blocking cost function $r(x, a)$ is defined with the parameters $b = 3$ ($> 1/\rho^{1+\rho} = 1.834$ as indicated above), $B = 5$, and $x_{\text{ref}} = 18$, giving $K^* = 17$. To help understand the optimisation problem, Figure 3 shows a plot of the long-run expected cost, $E^{\pi}(C) = \sum_{x \in S} p^{\pi}(x) \sum_{a \in A} r(x, a)\pi(a|x)$, both for the threshold-type deterministic policy $\pi = \bar{\pi}$ as a function of $K$, and for the parameterised stochastic piecewise-linear policy $\pi = \pi_{\theta}$ defined in (12).

The setup of the learning experiments is as follows: we choose the value $(x_{\text{ref}} + 10)$ for the initial blocking size guess, so that, already at the onset, blocking occurs rarely. Since the value of $K$ is no longer fixed (as was the case

<footnote>$^3$ The convexity of $g(\bar{\pi})$ when $b > 1/\rho^{1+\rho}$ is derived as follows: we regard $g(\bar{\pi})$ as a function $f(x)$ of a real-valued variable $x$ taking the place of the integer-valued $K$. We then write $f(x)$ as $Ah(x) + Bp^{\bar{\pi}}(x)$, where $A = Bb^{-x_{\text{ref}}} (1-\rho)$ and $h(x) = (\rho b)^x / (1 - \rho^{x+1})$. It is easy to see that $p^{\bar{\pi}}(x)$ is convex for all $x > 0$, so it suffices to show that $h(x)$ is convex when $b > 1/\rho^{1+\rho}$. We get: $h'(x) = \frac{h(x) \log(\rho b)}{1-\rho^{x+1}}$ and $h''(x) = \frac{h(x) \log^2(\rho b)}{(1-\rho^{x+1})^2} [1 + \rho^{x+1} \frac{\log b}{\log(\rho b)}]$. The root of the first derivative is $x^* = x_{\text{ref}} + \frac{\log(-\log \rho / \log(\rho b))}{\log(b)}$ as long as $b > 1/\rho$ (this guarantees that $\log(\rho b) > 0$, making the argument of the outer log positive when $\rho < 1$). Under this condition, the convexity of $h(x)$ is assured as long as $1 + \rho^{x+1} \log \rho / \log(\rho b) > 0$, or equivalently $b > 1/\rho^{1+\rho^{x+1}}, \forall x \ge 0$ (recall $\log(\rho b) > 0$). Since $1/\rho^{1+\rho^{x+1}} \le 1/\rho^{1+\rho}, \forall x \ge 0$ when $\rho < 1$, it suffices that $b > 1/\rho^{1+\rho}$, as claimed.</footnote>
where in the first inequality we use $EY \le \sqrt{EY^2}$ which follows from Cauchy-Schwartz inequality. We note, that $N_t = P_{\tilde{\delta}_{A^c}}(T_K > t) = 1 - P_{\tilde{\delta}_{A^c}}(T_K \le t) = 1 - F_K(t)$. Thus we have

$$ \frac{1}{D_A} E \int_0^\infty |\hat{N}_t - N_t| dt \le \frac{1}{\sqrt{N} E_{\tilde{\delta}_A} T_A} \int_0^\infty \sqrt{F_K(t)(1 - F_K(t))} dt. $$

Combining all of the above inequalities in an obvious manner, we obtain the bound from the thesis.

It follows from exponential ergodicity of $X_t^\tau$ that $\bar{T}_A, \bar{T}_K$ have exponential tails, that is, there exist constants $C_A, \lambda_A, C_K, \lambda_K$ such that $P_{\tilde{\delta}_A}(T_A > t) \le C_A \exp(-\lambda_A t)$ and $P_{\tilde{\delta}_{A^c}}(T_K > t) \le C_K \exp(-\lambda_K t)$. Therefore all the moments and the integral above are finite.
<figure><img src="image_5.png" /><figcaption>Figure 3: Curves showing the long-run expected cost (log scale) for the M/M/1 queue system when the blocking cost is defined as a function of the blocking size $x$, as: $5(1 + 3^{x-x_{\text{ref}}})$, with $x_{\text{ref}} = 18$. The blue dotted curve is given by (11), evaluated at integer $K$ values and represents the case of the threshold-type deterministic policy $\pi$ with rejection at $x = K$, as a function of integer-valued $K$. The violet solid curve is given by (14) on an equally-spaced grid of real $\theta$ values and represents the case of the parameterised stochastic piecewise-linear policy $\pi_{\theta}$ defined in (12). The minimizer of the blue curve is $K^* = 17$, a value that is very close to $x_{\text{ref}}$. The almost constant horizontal separation between the two curves is due to the relationship between $K$ and $\theta$ as $K = \lceil \theta \rceil + 1$. The queue load is $\rho = 0.7$.</figcaption></figure>

when estimating the blocking probability) but is now learned by the algorithm, it is not possible to choose a fixed value $J$ for the size of the absorption set $\mathcal{A}$. Instead, we consider a fixed $J/K$ fraction that adapts $J$ to each value of $K$ at the start of each learning step. In order to experiment with different sizes of the absorption set $\mathcal{A}$, we consider two different scenarios: $J = \lceil 0.3K \rceil$ and $J = \lceil 0.5K \rceil$, where $\cdot$ denotes the ceiling operator. Following the conclusions of the choice of $J$, $N$ and $T$ and their impact in the FV estimator accuracy outlined in Appendix D, given $J$ and using the approximate relative error formulas provided therein, we adjust $N$ and $T$ at each learning step to obtain approximate expected relative errors of $\epsilon_{\phi} \sim 100\%$ and $\epsilon_{ET} \sim 20\%$, respectively for the estimations of $\phi_t^J(K-1)$ and $E_{J-1}T_A$; i.e. we set the approximate expected error for $\hat{\phi}_t^J(K-1)$ much larger than the approximate expected error for $\hat{E}_{J-1}T_A$. For comparison purposes and to further check the conclusions of the error analysis in Appendix D, we also considered a scenario where the pre-defined relative errors are inverted, namely $\epsilon_{\phi} \sim 20\%$ and $\epsilon_{ET} \sim 100\%$.

For each setup we ran the FVRL policy learner on 100 learning steps. At each learning step, using (13) and as long as a positive FVEE estimate of the stationary probability $p^{\pi_{\theta}}(K-1)$ has been obtained, we estimate the sparse $\nabla_{\theta}g_x(\pi_{\theta})$ function defined in Section 3.3. To estimate the difference of the two $Q$ functions in (13), we use the coupling procedure described in Lemma 3, allowing up to 250 arrival events in each of 100 replications. The value of $\hat{\nabla}_{\theta}g_x(\pi_{\theta})$ is finally computed as the average $Q$-difference over these replications. Parameter $\theta$ is then updated by gradient descent using a constant learning rate of $\alpha = 10$, bounding any negative value of $\theta$ to 0.1$^4$. The estimated optimum threshold is set to $\hat{K}^* = \text{round}(\hat{\theta}^*) + 1$, where $\hat{\theta}^*$ is the value of the $\theta$ parameter obtained at the end of the learning process$^5$.

The benchmark MC learning algorithm uses the same policy gradient approach as FVRL with the only difference

<footnote>$^4$The value of $\theta$ is not allowed to go negative (although it could theoretically be as small as $-1$) because it would lead to a degenerate value of the $J$ parameter, namely $J=1$ which makes state $K-1$ be equal to $J-1$. Since $J-1$ is at the boundary of the absorption set $\mathcal{A}$, the stationary probability $p^{\pi_{\theta}}(K-1)$ would not be estimated by Fleming-Viot, thus precluding the estimation of the policy gradient. Finally, $\theta$ is lower bounded by 0.1 and not 0 in order to keep a real-valued $\theta$.</footnote>

<footnote>$^5$Note that the estimated optimum threshold is not set to $\lceil \hat{\theta}^* \rceil + 1$, as in the definition of the parameterised policy, so that $\hat{K}^*$ is more naturally chosen to be e.g. 5 when $\hat{\theta}^* = 4.01$, rather than 6.</footnote>
C Heuristics for the choice of the simulation parameters of the Fleming-Viot estimator

In this section we leverage theoretical results on the well-studied M/M/1/K queue system to provide insights about the
appropriate choice of parameter J –which in this unidimensional case fully defines the size of the absorption set A– and adapting
the values of the number of particles N and of the number of arrival events T (which directly impacts the number of cycles M
defined in Appendix B controlling, together with N, the convergence rate of the FV estimator) with the goal of obtaining an
accurate estimation of the blocking probability, i.e. of the stationary probability of state K. The conclusions of this analysis could
be used as an initial guideline for the choice of these parameters in a more general setting such as multidimensional problems.

We focus the analysis on the two quantities (out of the three involved in expression (7)) whose estimation is highly affected
by the trade-off of the choice of J described in remark 3, $\phi_t^J(K)$ and $\mathbb{E}_{J-1}(T_A)$, since the third quantity, $\mathbb{P}_J(T_K > t)$, does not
depend strongly on $J$ (see Appendix B)(8). The trade-off described in the remark states that, for constant values of $N$ and $T$,
the error in the estimation of $\phi_t^J(K)$ tends to decrease as $J$ increases, while the error in the estimation of $\mathbb{E}_{J-1}(T_A)$ tends to
increase. For the purpose of this impact analysis, we invert our reasoning: we fix $J$, and set the values of $N$ and $T$ required to
approximately satisfy predefined expected relative errors in the respective estimators. We then run simulations using each $N-T$
pair dictated by each relative error pair considered, and study statistics on the relative error obtained in the FV estimator of the
blocking probability as a function of the values of $N$ and $T$.

The values of $N$ and $T$ for each expected relative error pair are determined based on the following heuristics that approximately
express the two expected relative errors as a function of $J$, $N$, $T$, and the system's characteristics:

1. **Relative error of $\hat{\phi}_t^J(K)$**: Considering that in an $M/M/1/K$ queue system with $\rho < 1$, the conditioned blocking probability $\phi_t^{J=1}(K)$ converges as $t \to \infty$ towards $\sqrt{K}\rho^{K/2}$ [10], i.e. $\sim O(\rho^{K/2})$, similarly, the conditioned blocking probability for any absorption set size $J$, $\phi_t^J(K)$, increases to $\sim O(\rho^{(K-J+1)/2})$ as $t \to \infty$. As a very rough approximation and to get order of magnitude relations, we could think of $\phi_t^J(K)$ as the blocking probability $q$ of a single-server queue system with capacity equal to $\lceil \frac{(K-J+1)}{2} \rceil$, namely $q = \frac{(1-\rho)\rho^{K_A}}{1-\rho^{K_J+1}}$, where $K_J = \lceil (K-J+1)/2 \rceil$. Under this assumption, and ignoring the interdependence among FV particles, the relative error of estimating $\phi_t^J(K)$ with the empirical mean on $N$ samples is derived from the variance of a Binomial($q$) random variable as $\epsilon_\phi \sim \sqrt{(1-q)/N}q$. Thus, given $q$, the value of $N$ to approximately satisfy a desired expected relative error $\epsilon_\phi$ for $\hat{\phi}_t^J(K)$ is obtained as $N \sim \lceil \frac{1-q}{q\epsilon_\phi^2} \rceil \approx \lceil \frac{1}{q\epsilon_\phi^2} \rceil$, if $q \ll 1$ which is usually the case.

2. **Relative error of $\hat{\mathbb{E}}_{J-1}(T_A)$**: In an $M/M/1/K$ queue system with $\rho < 1$, the expected return time to the state $J-1$, when starting at $J-1$, is equal to $\mathbb{E}_{J-1}(T_{J-1}) = \frac{1}{\lambda_p(1+\rho^{-1})}$ [5], where $p$ is the stationary probability of the state $x=J-1$, i.e. $p = \frac{(1-\rho)\rho^{J-1}}{1-\rho^{K+1}}$. Given the absorption set $\mathcal{A} = \{0, 1, ..., J-1\}$, it can easily be shown that the expected return time to $\mathcal{A}$, when starting at $J-1$, is equal to $\mathbb{E}_{J-1}(T_A) = \mathbb{E}_{J-1}(T_{J-1})(1+\rho^{-1})$. If we want to observe $M$ return cycles to $\mathcal{A}$, we should simulate the queue system for as long as $t_S = M\mathbb{E}_{J-1}(T_A) = M\mathbb{E}_{J-1}(T_{J-1})(1+\rho^{-1})$, i.e. $t_S = \frac{M}{\lambda_S}$. Since $\lambda t_S$ is the expected number of arrival events $T$ observed in the time span $t_S$, we should simulate the system for as long as $T \sim M/p$ arrival events. On the other hand, given $M$ return cycles to $\mathcal{A}$, the standard error of the moment estimator of $\mathbb{E}_{J-1}(T_A)$ is given by $\sigma(T_A)/\sqrt{M}$. It is reasonable to assume (confirmed by experiments) that $\sigma(T_A) \sim \mathbb{E}_{J-1}(T_A)$, hence the relative error of the estimator of $\mathbb{E}_{J-1}(T_A)$, $\epsilon_{ET}$, is of the order of $1/\sqrt{M}$, which makes $M \sim 1/\epsilon_{ET}^2$. Thus, given $p$, the value of $T$ to approximately satisfy a desired expected relative error $\epsilon_{ET}$ for $\hat{\mathbb{E}}_{J-1}(T_A)$ is obtained as $T \sim \lceil \frac{1}{p\epsilon_{ET}^2} \rceil$.

From the above, the following common aspects are observed about the minimum values of $N$ and $T$ required to satisfy
predefined expected relative errors in the estimators of $\phi_t^J(K)$ and $\mathbb{E}_{J-1}(T_A)$:

1. $N$ affects the relative error of $\hat{\phi}_t^J(K)$ and $T$ affects the relative error of $\hat{\mathbb{E}}_{J-1}(T_A)$.

2. The relationship of $N$ and $T$ with their respective relative errors is of the same form, i.e. proportional to the inverse of a
stationary probability and to the inverse of the squared expected relative error.

On the other hand, the following difference is observed: for $\rho < 1$ and sufficiently large $K^9$, the value of $N$ as a function of the stationary probability $q$ is dominated by an increasing exponential function of $(K-J)/2$, i.e. $\sim O(\rho^{-(K-J)/2})$ whereas the value of $T$ as a function of the stationary probability $p$ is dominated by an increasing exponential function of $J$, i.e. $\sim O(\rho^{-J})$ which, importantly, does not depend on $K$. Thus, as mentioned in remark 3, the closer is $J$ to 0, the smaller the required $T$ and the larger the required $N$ for fixed expected relative errors, while the opposite is true when $J$ gets closer to $K$.

<footnote n="8">⁸Although the error of the estimator of $\mathbb{P}_J(T_K > t)$ depends on $N$, it is sensible to assume that its estimation error decreases as $N$ increases similarly or faster than the error of the estimator of $\phi_t^J(K)$, so we can limit our analysis to the error of the estimator of the latter.</footnote>

<footnote n="9">⁹$K$ is considered sufficiently large in this context when $\rho^{K-J}$ can be neglected w.r.t. 1.</footnote>
that it estimates the probability $p^{\pi\theta}(K - 1)$ in expression (13) using Monte-Carlo instead of Fleming-Viot, i.e. based on a single trajectory of $X_t^\pi$ as the ratio between the continuous time that the system spends at $K - 1$ and the total return cycle time to $J - 1$ (by renewal theory), where $J$ is defined as a function of $K$ as in FVRL. Each replication of the MC learner is started at $J - 1$ and is stopped when the average number of observed events over all learning steps in the respective FVRL replication is observed. These two conditions allow a fair comparison between the two methods.

The results of the above procedure run on 20 replications are shown in Figure 4. We see that FVRL clearly outperforms the benchmark MC learning in three out of the four scenarios considered, depicted in Figures 4(a), 4(b) and 4(d). In the scenario depicted in Figure 4(c), learning by the two methods is similar because the average number of events is about 10 times larger than the one observed in the other three scenarios (~ 40,000 vs. ~ 5,000), which allows MC to also observe the rare blocking state and thus learn almost as fast as FVRL. This is the same situation previously observed about the probability estimation results presented in Figure 2.

Of the two setups, ($\epsilon_\phi \sim 100\%$, $\epsilon_{ET} \sim 20\%$) used for Figures 4(a) and 4(b), and ($\epsilon_\phi \sim 20\%$, $\epsilon_{ET} \sim 100\%$) used for Figures 4(c) and 4(d), the former is the one that provides fastest learning (observe the steeper decrease of the curves in panels (a) and (b) compared to those in panels (c) and (d)), the smallest sample complexity, and the highest stability. This is consistent with the conclusions of the analysis presented in Appendix D of the choice of $J$, $N$ and $T$ in the FV estimation procedure, namely that, for an accurate estimation of the stationary probability, more importance should be given to achieve a smaller $\epsilon_{ET}$ error in the estimation of $E_{J-1}T_A$ than to achieve a smaller $\epsilon_\phi$ error in the estimation of $\phi_t^J(K)$.

Regarding the two values considered for the absorption set size, $J = 0.3K$ and $J = 0.5K$, results are very similar although using $J = 0.3K$ yields learning curves that tend to be closer to the optimum $\theta^*$ parameter.

Finally, in terms of stability we note that, for the least convenient setup of ($\epsilon_\phi \sim 20\%$, $\epsilon_{ET} \sim 100\%$) of Figure 4(d), in one replication the $\theta$ parameter suddenly increases from 28 to 140 (not shown but made apparent by the mean learning curve being significantly above the median learning curve). This overshoot impedes further learning because the blocking probability becomes extremely small at $K = 140$. This clearly illustrates the risks of allowing a too large error for the estimation of $E_{J-1}(T_A)$, which may yield a blocking probability largely overestimated (due to an underestimation of $E_{J-1}T_A$) generating such out-of-control excursions of the $\theta$ parameter.

Summing up, the best FVRL setup in terms of learning speed, sample complexity, and stability is to use an intermediate value of the $J/K$ factor, such as a value in the range [0.3K, 0.5K]. Once $J$ is defined, a large enough value of $T$ guaranteeing a small error in the estimation of $E_{J-1}T_A$ (say less than 50%) should be preferred over a large number of particles $N$ controlling the estimation error of $\phi_t^J(K - 1)$, which can be as large as 100%. Note that, although larger $J$ values suggest larger $T$ values, increasing $T$ doesn't affect as much the sample complexity as increasing $N$ does, as seen in the estimation problem in Section 4.1.1.
More importantly, experiments showed that the computational complexity –in terms of number of observed events– required to satisfy a given $\epsilon_{\phi}$ value (which impacts the number of observed events of the FV simulation) is much larger than the computational complexity required to satisfy the same value in $\epsilon_{ET}$ (which impacts the number of observed events in the single simulation of the $X_t^\pi$ process). Concrete values of the respective computational complexities can be derived from the results shown in Figure 2 and are complemented by the related observations presented within that Section 4.1.1. Thus, from the computational perspective alone, it is more convenient to favour a smaller error in $\hat{\mathbb{E}}_{J-1}(T_A)$ than a smaller error in $\hat{\phi}_t^J(K)$. Below we will see that this is also the case in terms of the estimation error of the FV estimator.

We completed the study of the appropriate choice of $N$ and $T$ by analyzing the impact of the errors of estimating $\hat{\phi}_t^J(K)$ and $\mathbb{E}_{J-1}(T_A)$ on the FV estimator of the blocking probability. To this end, we ran experiments on different combinations of expected relative errors $\epsilon_{\phi}$ and $\epsilon_{ET}$ and measured the accuracy in the estimation of the blocking probability of an $M/M/1/K$ system with $\rho = \lambda = 0.7$, $K = 20$, using a constant absorption set size of $J = 12$. The results of these experiments are shown in Figure 6 in terms of the estimation accuracy of the blocking probability $\hat{p}_{FV}(K)$.

From this heatmap, we conclude that it is more important to control the relative error in $\hat{\mathbb{E}}_{J-1}(T_A)$ than the relative error in $\hat{\phi}_t^J(K)$, as the contour lines are almost parallel to the axis where the relative error in $\hat{\phi}_t^J(K)$ is plotted, and their values indicate that the estimated blocking probability is larger than 1.5 times the true blocking probability when the relative error in $\hat{\phi}_t^J(K)$ is larger than 30% - 40%. Note that, for each experiment run, the FV estimator of the blocking probability was computed only when a minimum of 5 return cycles to $J-1$ were observed (these cycles are used to estimate the denominator $\mathbb{E}_{J-1}(T_A)$) after the 10 initial transitions of the system. This 10 transitions were used as a burn-in period to allow the system to get closer to the stationary regime. Therefore, if the number of arrival events $T$ is not large enough, the sample size on which the plotted median FV estimator is computed may be smaller than the 7 replications used for each $N-T$ combination.

In Figure 2 we used these heuristics to choose the different values of $N$ and $T$ on which the convergence properties of the FV estimator were analyzed: given $J = 12$, for the left plots (a) and (c), $T$ was fixed at the value associated to an approximate expected relative error in $\hat{\mathbb{E}}_{J-1}(T_A)$ equal to $\epsilon_{ET} = 20\%$, whereas the values of $N$ were chosen for approximate expected relative errors in $\hat{\phi}_t^J(K)$ equal to $\epsilon_{\phi} = 20\%, 10\%, 5\%$ for $K = 20$, and equal to $\epsilon_{\phi} = 80\%, 60\%, 40\%^{10}$; for the right plots (b) and (d), $N$ was fixed at the value associated to an approximate expected relative error in $\hat{\phi}_t^J(K)$ equal to $\epsilon_{\phi} = 60\%$, whereas the values of $T$ were chosen for approximate expected relative errors in $\hat{\mathbb{E}}_{J-1}(T_A)$ equal to $\epsilon_{ET} = 40\%, 20\%^{11}$.

<footnote n="10">¹⁰The larger relative errors chosen for $K = 40$ compared to $K = 20$ have to do with obtaining the same orders of magnitude for $N$ in each $K$ scenario.</footnote>

<footnote n="11">¹¹The expected relative error $\epsilon_{ET}$ depends on $K$ only through $1 - \rho^{K+1}$, and this term can be safely approximated by 1 for large enough values of $K$ such as 20 and 40 when $\rho = 0.7$. Thus, doing this approximation, the value of $T$ satisfying a given $\epsilon_{ET}$ only depends on $J$ and is thus the same for both $K$ scenarios.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 4: Learning curves for 20 replications of FVRL (green, ours) vs. Monte-Carlo learning (red) of the optimum parameter $\theta^* = 16$ (horizontal dashed gray line) of the parameterised acceptance policy $\pi_\theta$ of incoming jobs to the $M/M/1$ queue system with $\lambda = \rho = 0.7$. For each learning method (color), the dashed line indicates the mean learning curve, the marked solid line the median learning curve, the dark and light shades indicate, respectively, the location of 50% and 80% of learning curves closer to the median. Each plot caption indicates the following learning settings for the FVRL approach which impact the estimations of the quantities in (10) at each learning step: the approximate expected relative error wished for the estimation of $\phi_t^J(K)$ (which defines $N$, the number of particles at each learning step), the approximate expected relative error wished for the estimation of $\mathbb{E}_{J-1}T_A$ (which defines $T$, the number of arrival events for estimating $\mathbb{E}(T_A)$, at each learning step), and the fraction $J/K$ defining the size of the absorption set. The range of $N$ and $T$ values used, as well as the average number of events over learning steps and across replications are also reported. In all cases the learning parameter is kept constant at $\alpha = 10$, and the initial $\theta$ guess is 28.1, which defines a "large" blocking size of $K = 30$ whose small stationary probability ($\sim 10^{-6}$) makes blocking a rare event already at the onset of learning.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: Heatmap showing the impact of different combinations of the expected relative errors for $\hat{\phi}_t^J(K)$ and $\hat{E}_{J-1}(T_A)$ (indicated respectively on the vertical and horizontal axis) on the accuracy of the FV estimator of the blocking probability, $\hat{p}_{FV}(K)$, in an $M/M/1/K$ queue system over up to 7 valid experiments carried out for each combination. The system characteristics are $K = 20$, $\lambda = 0.7$, $\mu = 1$. The accuracy is shown as the ratio between the median $\hat{p}_{FV}(K)$ value and the true blocking probability $p(K)$. Thus, a ratio of 1 (light blue) implies 100% median accuracy, a ratio smaller than 1 (dark blue) implies underestimation, and a ratio larger than 1 (from cyan to red) implies overestimation. Selected contour levels are overlaid. The gray diagonal represents the line of equal expected relative errors in the two analyzed dimensions, and the gray points indicate the $5 \times 6 = 30$ error combinations on which experiments were run, whose size is proportional to the number of experiments (7 at the largest points, down to 2 at the smallest points). Smaller points are associated to a larger expected relative error in the estimation of $\hat{E}_{J-1}(T_A)$ which tend to preclude its estimation due to an insufficient number of observed return cycles to the absorption set $A$, as described in the text. The color scale is chosen as the transformation $\log_2(1+z)$ where $z$ is the ratio between the estimated and the true blocking probability; note that this transformation maps 0 to 0 and 1 to 1.</figcaption></figure>
## 4.2 Loss network system

Consider now a loss network that processes jobs of $I$ different classes which, as in the $M/M/1$ system, arrive following independent Poisson processes and are served with independently exponentially distributed times by one of $R$ available servers. We are interested in analyzing the number of jobs of each class being served by the system at a given time $t$, $X_t^\pi = (X_{t,1}, ..., X_{t,I})$, whose dynamics is governed by the system characteristics and the policy $\pi$ being applied. As before, we also consider a policy $\pi$ that belongs to the family of deterministic threshold-type admission control policies that reject an incoming job of a given class $i$ when the system is already serving a predefined number of jobs $K_i$ of that class, in which case a class-dependent rejection cost is accrued.

Let us denote by $\lambda_i, \mu_i, C_i, K_i$ the job arrival rate, the service rate, the rejection cost, and the blocking size of each job class $i = 1, ..., I$. Thus, the system can be represented by an MDP ($S = \{X \in \mathbb{N}^I : 0 \le X_i \le R, i = 1, ..., I \text{ s.t. } X^T 1 \le R, \}$, $U = \{0, 1\}$, $q = \{\lambda_1, ..., \lambda_I; \mu_1, ..., \mu_I\}$, $\mathcal{R} = r(\mathbf{x}, a) = \sum_{i=1}^I C_i \mathbf{1}_{\{x \in S_i, \mathcal{L}_{i,a=0}\}}$, $\bar{\Pi}_R$), where $S_i$ is the set of states where blocking occurs following an $i$-class job arrival denoted by $\mathcal{L}_i$, and $\bar{\Pi}_R$ is the class of threshold-type admission control policies for a system with $R$ servers.

To facilitate illustration, we consider in our experiments the smallest multi-class loss network, namely one that serves just two classes of job.

### 4.2.1 Estimation of the expected cost using FVEE

Contrary to the $M/M/1$ case, estimating the expected cost in the loss network system is not equivalent to estimating the blocking probability, because blocking occurs at more than one state whose cost in general depends on the class of the rejected job.

The expression of the long-run expected cost $\mathbb{E}^\pi(C)$ (thereafter "expected cost") under the threshold policy $\pi \in \bar{\Pi}_R$ and costs $\mathcal{R}$ is obtained from the total probability theorem applied on all possible job arrival classes and all possible accept/reject actions on those arrivals, as:

$$ \mathbb{E}^\pi(C) = \sum_{x \in S} p^\pi(\mathbf{x})c(\mathbf{x}), \qquad (15) $$

where $p^\pi(\mathbf{x}) \doteq \mathrm{Pr}(\mathbf{X}_n^\pi = \mathbf{x})$ is the stationary probability that the SMDP $\mathbf{X}_n^\pi$ (defined in Section 4.1.2) is at state $\mathbf{x}$, and $c(\mathbf{x}) = \sum_{i=1}^I C_i \lambda_i / \Lambda$ $[\mathbf{1}_{\{\sum_{j=1}^I x_j = R\}} + \mathbf{1}_{\{\sum_{j=1}^I x_j < R\}} \mathbf{1}_{\{x_i = K_i\}}]$ is the expected cost of rejection at $\mathbf{x}$ over all possible arriving job classes, where $\Lambda \doteq \sum_{i=1}^I \lambda_i$ is the total job arrival rate, and the indicator terms reflect, respectively, that rejection can happen either when the network is operating at full capacity or when it is serving the maximum allowed number of jobs of the arriving class. To better isolate the impact of the estimation quality of the stationary probability $p^\pi(\mathbf{x})$ on the expected cost estimation, we use the true values of the job arrival rates $\lambda_i$ to compute the expected cost $c(\mathbf{x})$, although in practice these values are normally unknown and would have to be estimated from the simulation. Finally, the blocking probability is given by the same expression (15) when setting $C_i = 1$ in $c(\mathbf{x})$, for all $i = 1, ..., I$.

A loss network with the following characteristics was considered for the estimation of the expected cost: capacity $R = 6$ servers, $\lambda = [1, 5]$, $\mu = [3.33, 50.0]$ (hence $\rho = [0.3, 0.1]$), $C = [2 \times 10^3, 2 \times 10^5]$, $K = [4, 6]$. The choice of the rejection costs allows us to illustrate the benefits of the FV approach over MC as it makes a few of the cost-generating states with smaller probability as important as those with larger probability in terms of their contribution to the system's expected cost. These contributions are listed in Table 1 where states are listed in decreasing order of their respective true stationary probabilities, computed following the product form of the stationary distribution of the stochastic knapsack described in [24, chapter 4]. We observe that the top three states in terms of probability contribute to about 72% of the expected cost, while the bottom four, with probability smaller than $10^{-6}$, contribute to as much as 28%, therefore obtaining accurate estimates of the probability of a few of those smaller probability states is important for an accurate estimation of the expected cost.
D FV and FVRL algorithms

This section presents the two algorithms used throughout this work to apply the Fleming-Viot methodology to network systems in order to (i) estimate the expected rejection cost (FV algorithm) and (ii) learn the optimum blocking sizes (FVRL algorithm).
<table><caption>Table 1: Contribution to the expected cost $E^π(C)$ by each possible blocking state in the loss network system according to expression (15), with $R = 6$, $K = [4, 6]$, $\lambda = [1, 5]$, $\rho = [0.3, 0.1]$, $C = [2 \times 10^3, 2 \times 10^5]$, sorted by decreasing probability $p^π(x)$. The bottom four states with occurrence probability smaller than $10^{-6}$ contribute to ~28% of the expected cost.</caption><thead><tr><th>State x</th><th>Prob. p<sup>π</sup>(x)</th><th>Expected cost c(x)</th><th>p<sup>π</sup>(x)c(x)</th><th>% Exp. cost E<sup>π</sup>(C)</th></tr></thead><tbody><tr><td>[4, 0]</td><td>2.2 × 10<sup>-4</sup></td><td>0.333 × 10<sup>3</sup></td><td>0.1889</td><td>20.0%</td></tr><tr><td>[4, 1]</td><td>2.3 × 10<sup>-5</sup></td><td>0.333 × 10<sup>3</sup></td><td>0.0075</td><td>2.0%</td></tr><tr><td>[4, 2]</td><td>1.1 × 10<sup>-6</sup></td><td>167 × 10<sup>3</sup></td><td>0.0754</td><td>49.7%</td></tr><tr><td>[3, 3]</td><td>5.0 × 10<sup>-7</sup></td><td>167 × 10<sup>3</sup></td><td>0.0839</td><td>22.0%</td></tr><tr><td>[2, 4]</td><td>1.3 × 10<sup>-7</sup></td><td>167 × 10<sup>3</sup></td><td>0.0210</td><td>5.5%</td></tr><tr><td>[1, 5]</td><td>1.7 × 10<sup>-8</sup></td><td>167 × 10<sup>3</sup></td><td>0.0028</td><td>0.74%</td></tr><tr><td>[0, 6]</td><td>9.3 × 10<sup>-10</sup></td><td>167 × 10<sup>3</sup></td><td>0.0002</td><td>0.04%</td></tr><tr><td>Total</td><td></td><td></td><td>7.214</td><td>100.0%</td></tr></tbody></table>

The true blocking probability and true expected cost are computed from (15), against which the estimations by the analyzed methods are compared. The quality of their estimations is visualized using the violin plots presented in Figure 5 obtained from 10 replications of the estimation process. For each panel, the green violin plots on the left column of the figure analyze the convergence of FVEE, the FV expectation estimator of the quantity indicated on the vertical axis, in terms of increasing $T^6$ while keeping the number of particles $N$ fixed at 400. The absorption set $A$ is chosen as all the states whose server occupation is strictly smaller than 2 and 3, respectively for job classes 1 and 2, i.e. $A = \{[0, 0], [0, 1], [0, 2], [1, 0], [1, 1], [1, 2]\}$. On the other hand, the red violin plots in each panel on the right column of Figure 5 present the results of the benchmark estimator, which is based on an MC simulation of the loss network. In order to make a fair comparison, the start state for MC is chosen uniformly at random in $\tilde{\partial}A = \{[0, 2], [1, 2], [1, 0], [1, 1]\}$ –so that both FVEE and MC methods start at a similar "distance" from the set of blocking states– and the MC simulation is let run until the number of events observed in the FVEE simulation is reached, whose average over replications is indicated at the top horizontal axes of the plots.

From Figure 5 we observe that both FVEE and MC successfully estimate the blocking probability and the expected cost in terms of obtaining a distribution of estimates that include their respective true values. However, FVEE shows a consistently smaller error than MC –as indicated by the median values in each experiment pair– as well as a consistently smaller variability.

### 4.2.2 Learning blocking sizes with FVRL

The FVRL algorithm for learning the blocking sizes $K^* \in \mathbb{N}^I$ of a loss network serving $I$ jobs classes learns the values of $\theta_i$ of $I$ parameterised acceptance policies, each of the form (12). For each $\theta_i$ observed during learning, the deterministic blocking size $K_i$ of the $i$-th policy is defined as $K_i = \lceil \theta_i \rceil + 1$, but, as in the $M/M/1$ case, the estimated blocking size obtained at the end of the learning process is defined as $\hat{K}_i^* = \text{round}(\hat{\theta}_i^*) + 1$.

Upon arrival of a job of class $i$, the respective $\pi_{\theta_i}$ acceptance policy is applied, making the system’s acceptance policy equal to $\pi_{\theta_i}(a = 1|x) = \sum_{i=1}^I \pi_{\theta_i}(a = 1|x_i)1_{L_i}$, where $L_i$ denotes the arrival event of a class-$i$ job. Then, the policy derivative w.r.t. $\theta_i$ is non-zero only at states $x$ for which $x_i = K_i - 1$ as long as they satisfy the $R$-server constraint for a possible job acceptance, $x^T1 < R$. At those states the derivative is equal to +1 for action $a = 1$ (accept) and -1 for action $a = 0$ (reject).

<footnote n="6">⁶We chose the number $T$ of arrival events for estimating $E(T_A)$ for the convergence analysis –instead of choosing the number of particles $N$– based on the conclusion obtained in Appendix D about the higher impact of $T$ over $N$ in the estimation accuracy of the stationary probability in the $M/M/1/K$ queue system.</footnote>
## FV algorithm

### Data:

A. System characteristics: a loss network with $R$ servers serving jobs of $I$ different classes whose state is represented by the number of jobs of each class in the system, $\mathbf{x} = (x_1, ..., x_I)$.

B. System dynamics: the job arrival rate $\lambda_i$ and the service rate $\mu_i$ for each job class $i = 1, ..., I$ are known.

C. Job acceptance policy: an incoming job of class $i$ is either accepted ($a = 1$) or rejected ($a = 0$). It is accepted whenever the system is not operating at full capacity $R$ and when the number of jobs of the arriving class being served by the system is less than a constant $K_i$, $i = 1, ..., I$. Otherwise, it is rejected, in which case a cost $C_i$ is accrued. Thus, the job acceptance policy, when the system is at state $\mathbf{x}$ and a job of class $i$ arrives ($\mathcal{L}_i$), is $\pi(a = 1|\mathbf{x}, \mathcal{L}_i) = \mathbf{1}_{\{\sum_{j=1}^I x_j < R\}} \mathbf{1}_{\{x_i < K_i\}}$, and the set of blocking states is $\mathcal{C} = \{\mathbf{x} : \sum_{j=1}^I x_j = R \text{ or } \sum_{j=1}^I x_j < R, x_i = K_i \text{ for some } i = 1, ..., I\}$.

D. $J_i, i = 1, ..., I, 0 \le J_i \le K_i$: size of the absorption set $\mathcal{A}$ in dimension $i$, that is whenever the state of an FV particle visits a state $\mathbf{x}$ having $x_i = J_i - 1$ for some $i$, the particle is considered absorbed.

E. $N$: number of FV particles used to estimate $\phi_t^{\tilde{\mathcal{A}}^c}(K)$ and $\mathbb{P}_{\tilde{\mathcal{A}}^c}(T_K > t)$ in (3).

F. $T$: number of arrival events (incoming jobs), which has a direct impact on the number of cycles $M$ used to estimate $\mathbb{E}_{\tilde{\mathcal{A}}^c}(T_A)$ in (3).

G. $B$: number of burn-in time steps to assure stationarity, typically 10-20.

H. $M_0$: minimum number of return cycles to $\mathcal{A}$ to have a reliable estimate of $\mathbb{E}_{\tilde{\mathcal{A}}^c}(T_A)$, typically 5-10.

**Result:** An estimate of the expected rejection cost under stationarity.

**Steps:** The algorithm is divided into the following three steps:

1. Estimation of $\mathbb{E}_{\tilde{\mathcal{A}}^c}(T_A)$:

(a) Simulate the continuous-time Markov process $\mathbf{X}_t^\pi$ starting at $\mathbf{X}_0^\pi$ chosen uniformly at random from the set $\tilde{\mathcal{A}} = \{\mathbf{x} : x_i = J_i - 1 \text{ for some } i = 1, ..., I\}$.

(b) Record the observed values of the cycle time $T_A$ defined in Section 3.1, measured as return times to the set $\mathcal{A}$ after the first visit to a state in $\tilde{\mathcal{A}}$ following the burn-in period given by parameter $B$.

(c) Record the states at which the process enters $\mathcal{A}^c$.

(d) Stop the simulation when the number of arrival events is equal to $T$, and record the number $M$ of observed complete return cycles to $\mathcal{A}$ after the burn-in period.

(e) If $M > M_0$, estimate $\mathbb{E}_{\tilde{\mathcal{A}}^c}(T_A)$ as $\frac{1}{M} \sum_{k=1}^M T_{A,k}$, o.w. conclude that the expected rejection cost cannot be reliably estimated, and end the process here.

(f) Estimate the entrance distribution to $\mathcal{A}^c$, $\hat{p}_{\tilde{\mathcal{A}}^c}^\pi(\mathbf{x})$ as $\frac{1}{M} \sum_{k=1}^M \mathbf{1}_{\mathbf{x}(k)=\mathbf{x}}$, $\forall \mathbf{x} \in \tilde{\mathcal{A}}^c$, where $\overline{M}$ is the number of observed entrances to $\mathcal{A}^c$ at states $\mathbf{x}(k)$.

2. Estimation of $\mathbb{P}_{\tilde{\mathcal{A}}^c}(T_K > t)$ and $\phi_t^{\tilde{\mathcal{A}}^c}(\mathbf{x}_K)$, where $\mathbf{x}_K$ is any state in the set of blocking states $\mathcal{C}$:

(a) Simultaneously simulate $N$ trajectories (FV particles) that independently follow the law of the Markov process $\mathbf{X}_t^\pi$, starting at $\mathbf{X}_0^\pi$ chosen in $\tilde{\mathcal{A}}^c$ following $\hat{p}_{\tilde{\mathcal{A}}^c}^\pi$. Record the observed values of the first killing times $T_K$ at which each of the $N$ particles enters $\mathcal{A}$.

(b) Every time one of the particles is killed at $\mathcal{A}$, reinitialize it to the state of one of the other $N - 1$ particles, selected uniformly at random.

(c) Stop when all $N$ particles have been killed at least once, as this is the maximum time $t$ that will contribute to the integral in expression (3).

(d) For each $t$ in the set $\{T_{K,k}\}_{k=1...N}$, estimate $\mathbb{P}_{\tilde{\mathcal{A}}^c}(T_K > t)$ as $\frac{1}{N} \sum_{k=1}^N \mathbf{1}_{T_{K,k} > t}$.

(e) For each time $t$ at which one of the particles changes state, estimate $\phi_t^{\tilde{\mathcal{A}}^c}(\mathbf{x}_K)$ as the proportion of particles that are at state $\mathbf{X}_t^\pi = \mathbf{x}_K$.

3. Estimation of the expected rejection cost:

Estimate the integral in (3) where $\eta(\mathbf{x})$ is the expected cost accrued in state $\mathbf{x} \in \mathcal{C}$, namely

$$ \eta(\mathbf{x}) = \sum_{i=1}^{I} C_i \lambda_i / \Lambda \left[ \mathbf{1}_{\{\sum_{j=1}^{I} x_j = R\}} + \mathbf{1}_{\{\sum_{j=1}^{I} x_j < R\}} \mathbf{1}_{\{x_i = K_i\}} \right], \text{ where } \Lambda = \sum_{i=1}^{I} \lambda_i \text{ is the total job arrival rate. The integral can} \\ \text{easily be computed as a finite sum of the piecewise constant function of time resulting from the product } \hat{\mathbb{P}}_{\tilde{\mathcal{A}}^c}(T_K > t) \hat{\phi}_t^{\tilde{\mathcal{A}}^c}(\cdot), \text{ which} \\ \text{is 0 for } t > \max(\{T_{K,k}\}_{k=1...N}). $$

**Algorithm 1:** FV algorithm for the estimation of the expected rejection cost in an $M/M/I/R$ loss network serving $I$ different job classes with $R$ servers, of which the $M/M/1/K$ queue system is a particular case.
<figure><img src="image_7.png" /><figcaption>Figure 5: Violin plots showing the convergence properties of the Fleming-Viot (left, green) (ours) and Monte-Carlo (right, red) estimators of the blocking probability (top) and of the expected cost (bottom) in a two-job-class loss network with $R = 6$ servers, $\lambda = [1, 5]$, $\rho = [0.3, 0.1]$, blocking at $K = [4, 6]$ with costs $C = [2 \times 10^3, 2 \times 10^5]$, as the number $T$ of arrival events for estimating $\mathbb{E}(T_A)$ increases. The absorption set sizes by job class are $J = [2, 3]$ and the number of particles is $N = 400$. The horizontal tick inside each violin plot is the median estimated value over 10 replications. In the corresponding MC experiments, the average number of return cycles to the absorption set is used on the horizontal axis, which is a measure of sample size. Finally, the top horizontal axes show the average number of events observed in the experiments run in each set, which by design coincide between each paired FV-MC execution, as mentioned in the text.</figcaption></figure>

Hence, the partial derivative of $g(\pi_\theta)$ w.r.t. $\theta_i$ in (13) becomes:

$$ \frac{\partial g(\pi_\theta)}{\partial \theta_i} = \sum_{\mathbf{x} \in S_i} p^{\pi_\theta}(\mathbf{x}) [Q^{\pi_\theta}(\mathbf{x}, 1) - Q^{\pi_\theta}(\mathbf{x}, 0)], \quad \text{for } i = 1, \dots, I, \qquad (16) $$

where $S_i$ is the set of states $\mathbf{x}$ where $x_i = K_i - 1$ and $\mathbf{x}^T \mathbf{1} < R$.

To illustrate the FVRL algorithm, we consider the same loss network used for the estimation problem (see previous Section 4.2.1). For this choice of parameters, the optimum blocking sizes of the threshold policy are equal to $K^* = [4, 6]$ (i.e. $\theta^* = [3, 5]$), found by evaluating the true expected cost on all possible threshold combinations. Note that $\lambda_i$ and $\mu_i$ do not need to be known to apply the algorithm.

The setup of the learning experiments is as follows: we set the initial blocking size guess to $K = [0.1, 0.1]$, and the following simulation parameters for the estimation of the stationary probabilities $p^{\pi_\theta}(\mathbf{x})$ for each $\mathbf{x} \in S_i$: $J/K = 0.5$ for both classes, the number of FV particles $N = 100$ (four times smaller than the value of $N$ used in the estimation problem), and the number of arrival events for estimating $\mathbb{E}(T_A)$, $T = 500$ (very close to the smallest $T$ value considered in the estimation problem). The start state of each FV particle is chosen following the estimated entrance distribution of the states in $\vec{\partial}A^c$. The estimation of the gradient in (16) is completed by estimating the $Q$-difference at each $\mathbf{x}$ following the same coupling setup of Lemma 3 used for the $M/M/1$ case. The learning rate $\alpha$ is set initially to 1.0 and then decreased inversely proportional to the learning step. Finally, at each learning step, the estimated $\theta_i$ values are bounded to the interval $[0.1, K_i - 1 - 0.1]$.