**FVRL algorithm**

**Data:**

A. Characteristics of the loss network to optimise: arrival rate $\lambda_i$ and service rate $\mu_i$ of each job class $i = 1, \dots, I$.

B. Cost of rejecting an arriving job class, $C_i$, $i = 1, \dots, I$.

C. $\pi_{\theta_i}$ ("accept" $|x_i|$): the job acceptance policy for an arriving job of class $i$ parameterised by the positive real-valued $\theta_i$, as defined in expression (8).

D. $\theta = (\theta_i)_{i=1,\dots,I}$: positive non-integral initial values of the parameter to optimise, from where $K_i$ can be obtained as $K_i = \text{ceiling}(\theta_i + 1)$.

E. $F_i$: the size of the absorption set $\mathcal{A}$ in dimension $i$ as a fraction of $K_i$.

F. $L$: number of learning steps, i.e. the number of times an update of $\theta$ will be computed by the gradient-based algorithm.

G. $N, T$: respectively, number of FV particles and number of arrival events used to estimate the blocking probability with the FV estimator described in Algorithm 1.

H. $U, S$: respectively, number of replications and maximum number of time steps allowed to estimate the function $\eta(x)$ defined in Section 3.4, whose expectation is the gradient of the average state value. Ex: $R = 100, S = 250$.

**Result:** An estimate of the optimum blocking sizes $\hat{K}_i^*, i = 1, \dots, I$.

**Steps:**

1. Compute the deterministic blocking sizes $K_i = \text{ceiling}(\theta_i + 1), i = 1, \dots, I$. Set $J_i$, the size of the absorption set $\mathcal{A}$ in dimension $i$ as $J_i = \lceil F_i K_i \rceil$.

2. Estimate the stationary probabilities $p^{\pi_{\theta}}(\mathbf{x})$ for each $\mathbf{x}$ such that $x_i = K_i - 1$ and $\sum_{j=1}^I x_j < R$ using the FV algorithm described in Algorithm 1.

3. Simulate $U$ times two coupled systems, as described in Section 3.4, for as long as $S$ time steps to estimate the $Q$ differences $Q^{\pi_{\theta}}(\mathbf{x}, 1) - Q^{\pi_{\theta}}(\mathbf{x}, 0)$, for each $\mathbf{x}$ considered in the previous step, as the average of the difference observed on the replications where mixing of the two systems occurs.

4. Estimate the gradient of the average state value $\frac{\partial v^{\pi_{\theta}}}{\partial \theta_i}$ given in expression (11), using the estimates of steps (2) and (3).

5. Update $\theta$ following the classical gradient descent algorithm.

6. Repeat steps (1)-(5) until the number of learning steps $L$ is reached.

7. Use the final value $\theta^L$ to compute the optimum blocking sizes estimated by the algorithm, $\hat{K}_i^*$, as $\text{round}(\theta_i^L) + 1, i = 1, \dots, I$.

**Algorithm 2:** FVRL algorithm for finding the optimum blocking sizes $K_i$ 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.
The results of the above procedure, run on 20 replications, are shown in Figure 6 and described in the caption. Besides showing the FVRL and MC learning curves, the figure includes a reference learning path depicted in blue on the left, where the learning process uses the *true* stationary probabilities (presented in Table 1) in the gradient expression (16), so that only the *Q*-differences need to be estimated for the computation of the gradient. The aim is to have a reference for comparison that corresponds to the most ideal available learning scenario.

We observe that FVRL clearly outperforms MC learning as, in median, parameter $\theta_1$ is correctly learned by FVRL but not by MC, where it practically remains stuck near 5. FVRL also presents smaller variability in the learning paths than MC. As seen in the center panel in Figure 6(a)), however, 30 learning steps were not enough for all 20 replications to reach the optimal $\theta_1$ value, although most of them started learning.
<figure><img src="image_8.png" /><figcaption>Figure 6: Two comparisons of the FVRL method with Monte-Carlo learning on 20 replications. The top row shows each observed learning trajectory and the bottom row shows their distribution. On each row, the left panel represents the scenario closest to the true learning process as it corresponds to learning using the true stationary probabilities, where only the difference in the *Q* values is estimated. The center panel corresponds to FVRL and the right panel to Monte-Carlo learning. The plots show the results over 20 experiments run over 30 learning steps on a loss network with *R* = 6 servers, *λ* = [1, 5], *ρ* = [0.3, 0.1], 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.5, 0.5], the number of particles at *N* = 100 and the number of arrival steps for estimating E(*T*<sub>A</sub>) 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 80% 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 Conclusion

We presented the Fleming-Viot particle system (FV) as an alternative to Monte-Carlo for a more efficient exploration of reinforcement learning environments by leveraging prior knowledge about a set $\mathcal{A}$ of frequently visited states with no rewards. The method helps to accelerate the discovery of rarely observed states that may yield large rewards affecting optimisation objectives, thus speeding up policy learning. It is proposed particularly for situations that preclude the use of offline exploration techniques, such as importance sampling, and to avoid the need for expert knowledge to guide learning.

The use of FV is two-fold: to estimate a long-run optimisation objective, and to learn a policy for control problems. Its use was illustrated on the design of threshold-type admission control policies to minimize the long-run expectation cost in two stochastic network environments: an $M/M/1$ queue system and a loss network accepting two classes of job.

For the estimation problem, the FV estimation approach proved to be much more efficient than Monte-Carlo in the $M/M/1$ system for large capacities $K$, where Monte-Carlo completely fails using the same event budget. In the loss network system, FV was able to estimate the expected cost more accurately and precisely than Monte-Carlo.

For the control problem, the presented FVRL method –which uses FV to estimate sparse gradients of a policy gradient algorithm– is able to find the optimal thresholds of admission control policies considerably faster than Monte-Carlo, in both the $M/M/1$ and the loss network cases. We note that in the control problem the estimation accuracy of the gradient 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. We also note that the FVRL algorithm appears to converge to a global optimum in these problems, even though the theoretical convergence properties of Proposition 4 only guarantee reaching a local optimum, and even if not all convergence requirements are satisfied (e.g. increasing number of cycles $M$ and number of particles $N$).

The presented experimental results showed the validity of the FV methodology in rather small problems, namely on finite state and action spaces whose cardinality does not go over a couple of hundreds. Proving its feasibility on larger environments is left for future work.

We finally note that the presented FV methodology is applicable insofar as a “black-box” simulator or emulator of the system under study is available, which must be able to provide the next state and reward observed when the agent takes any possible action at any possible system’s state. This requirement also means that FV is currently not applicable to online or buffer replay contexts, i.e. when trajectories are collected from an ongoing or previously observed agent-environment interaction.

In future work we will generalize the algorithms so that they include a step where the set $\mathcal{A}$ of frequently visited states is defined from an initial exploration of the environment, as well as allow for the possibility of non-zero rewards received in $\mathcal{A}$.

We also intend to extend the algorithms to more classical RL environments, such as labyrinths, games, mountain car, etc., which are episodic in nature, and to optimisation problems with discount.
## References

[1] Alekh Agarwal, Sham M Kakade, Jason D Lee, and Gaurav Mahajan. On the theory of policy gradient methods: Optimality, approximation, and distribution shift. Journal of Machine Learning Research, 22(98):1–76, 2021.

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

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

[4] Dimitris Bertsimas and Cheol Woo Kim. Optimal control of multiclass fluid queueing networks: A machine learning approach. CoRR, abs/2307.12405, 2023.

[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] Ernesto Garcia, Paola Bermolen, Matthieu Jonckheere, and Seva Shneer. Probabilistic insights for efficient exploration strategies in reinforcement learning, 2025.

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

[13] Isaac Grosof, Siva Theja Maguluri, and R Srikant. Convergence for natural policy gradient on infinite-state queueing mdps. arXiv preprint arXiv:2402.05274, 2024.

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

[15] Mark Gluzman J. G. Dai. Queueing network controls via deep reinforcement learning. Stochastic Systems, 12(1):30–67, 2022.

[16] Chi Jin, Akshay Krishnamurthy, Max Simchowitz, and Tiancheng Yu. Reward-free exploration for reinforcement learning. In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 4870–4879. PMLR, 13–18 Jul 2020.

[17] Antonio Massaro, Francesco De Pellegrini, and Lorenzo Maggi. Optimal trunk-reservation by policy learning. In IEEE INFOCOM 2019, apr 2019.
[18] Maja J Mataric. Reward functions for accelerated learning. In *Machine learning proceedings 1994*, pages 181–189. Elsevier, 1994.

[19] Eric Moulines and Francis Bach. Non-asymptotic analysis of stochastic approximation algorithms for machine learning. In J. Shawe-Taylor, R. Zemel, P. Bartlett, F. Pereira, and K.Q. Weinberger, editors, *Advances in Neural Information Processing Systems*, volume 24. Curran Associates, Inc., 2011.

[20] 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.

[21] Brahma S. Pavse, Matthew Zurek, Yudong Chen, Qiaomin Xie, and Josiah P. Hanna. Learning to stabilize online reinforcement learning in unbounded state spaces. In *Proceedings of the 41st International Conference on Machine Learning*, ICML’24. JMLR.org, 2024.

[22] M. L. Puterman. *Markov Decision Processes: Discrete Stochastic Dynamic Programming*. John Wiley & Sons, 1994.

[23] Herbert E. Robbins. A stochastic approximation method. *Annals of Mathematical Statistics*, 22:400–407, 1951.

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

[25] Richard Sutton and Andrew Barto. *Reinforcement Learning, an introduction*. MIT Press, 2018.

[26] 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.

[27] Andrea Tirinzoni, Andrea Sessa, Matteo Pirotta, and Marcello Restelli. Importance weighted transfer of samples in reinforcement learning. In *International Conference on Machine Learning*, pages 4936–4945. PMLR, 2018.

[28] Martin J Wainwright. *High-dimensional statistics: A non-asymptotic viewpoint*, volume 48. Cambridge University Press, 2019.
A Estimation of the stationary probability distribution $p^\pi(x)$

using Fleming-Viot particle systems Based on expression (3) of Proposition 1, the FV estimator of the stationary probability $p^\pi(x)$ for states $x$ outside the absorption set $\mathcal{A}$ is constructed as follows: we first define functions $f(x, t)$ and $g(t)$ as:

$$ 
\begin{aligned}
f(x,t) &\doteq \phi_t^{\vec{\partial} \mathcal{A}^c}(x), \\
g(t) &\doteq \frac{\mathbb{P}_{\vec{\partial} \mathcal{A}^c}(T_K > t)}{\mathbb{E}_{\vec{\partial} \mathcal{A}} T_A},
\end{aligned}
$$

where $\phi_t^{\vec{\partial} \mathcal{A}^c}(x)$ is defined in Proposition 1 and the other concepts are defined in the Definitions Section 3.1. These functions, together with their estimators $\hat{f}$ and $\hat{g}$, will be instrumental in Appendix B in proving the consistency of the FV estimator of the stationary probability provided by Theorem 2.

Using the estimators $\hat{f}(x, t)$ and $\hat{g}(t)$ described below, the FV estimator of the stationary probability is computed as:

$$ \hat{p}_{FV}^{\pi}(x) = \int_{0}^{\infty} \hat{f}(x, t) \hat{g}(t) dt, \quad \forall x \in \mathcal{A}^{c}. \qquad (17) $$

Firstly, we explain the construction of the estimator $\hat{g}(t)$, which is a ratio of the quantities $\mathbb{P}_{\vec{\partial} \mathcal{A}^c}(T_K > t)$ and $\mathbb{E}_{\vec{\partial} \mathcal{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} \mathcal{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} \mathcal{A}^c$ following the entrance state distribution into $\mathcal{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} \mathcal{A}} T_A$ while simulation (ii) contributes to the estimation of $\mathbb{P}_{\vec{\partial} \mathcal{A}^c}(T_K > t)$.

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

$$ 
\begin{aligned}
\tau_{\mathcal{A}^c, i} &= \inf_{t > \tau_{\mathcal{A}, i-1}} \{X_t^\pi \in \mathcal{A}^c\}, \\
\tau_{\mathcal{A}, i} &= \inf_{t > \tau_{\mathcal{A}^c, i}} \{X_t^\pi \in \mathcal{A}\},
\end{aligned}
$$

for $i \ge 1$. We also define $T_{E,i} = \tau_{\mathcal{A}^c,i} - \tau_{\mathcal{A},i-1}$, the entry time into $\mathcal{A}^c$ within cycle $i$, and $T_{K,i} = \tau_{\mathcal{A},i} - \tau_{\mathcal{A}^c,i}$, the killing time within cycle $i$. We note that $T_{\mathcal{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} \mathcal{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 $\mathcal{A}^c$, which can be considered as an estimator of $\hat{p}_{\vec{\partial} \mathcal{A}^c}^\pi$, as $\hat{p}_{\vec{\partial} \mathcal{A}^c}^\pi = \frac{1}{M} \sum_{i=M_0+1}^{M_0+M} \mathbf{1}_{X_{\tau_{\mathcal{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} \mathcal{A}^c$ randomly chosen according to the estimated stationary distribution $\hat{p}_{\vec{\partial} \mathcal{A}^c}^\pi$ (this is trivial when $\vec{\partial} \mathcal{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 (18) $$

Next, we explain the construction of the estimator $\hat{f}(x, t)$, which is the conditional occupation probability $\phi_t^{\vec{\partial} \mathcal{A}^c}(x)$, defined in Proposition 1, at each $x \in \mathcal{A}^c$ for all times $t$, 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 $\mathcal{A}$ as the absorption set and starting each particle $i$ at $x_i$ a randomly chosen state according to the distribution $\nu(y) = p_{\tilde{\partial}\mathcal{A}^c}(y)\mathbf{1}_{y\in\tilde{\partial}\mathcal{A}^c}$. That is, all particles start at the boundary of $\mathcal{A}^c$, as required by the quantity to estimate, $\phi_t^{\tilde{\partial}\mathcal{A}^c}(x)$. 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$. Since $m(x, \xi_t^\nu)$ is an estimator of $\phi_t^{\tilde{\partial}\mathcal{A}^c}(x)$ (because $\nu$ is restricted to the boundary of $\mathcal{A}^c$ where it is equal to the entrance state distribution into $\mathcal{A}^c$ under stationarity), $f$ can be estimated by

$$ \hat{f}(x, t) = m(x, \xi_t^\nu). \qquad (19) $$

We note that, by construction, $\hat{g}(t) = 0$ for $t > T_{\mathcal{K},\max} = \max\{T_{\mathcal{K},i} : 1 \le i \le N\}$. Therefore, since we wish to compute $\hat{p}^\pi(x) = \int_0^\infty \hat{f}(x,t)\hat{g}(t)dt$, we only need to simulate the Fleming-Viot process until time $T_{\mathcal{K},\max}$ is reached. Also, since both $\hat{f}(x,t)$ and $\hat{g}(t)$ are almost surely piecewise constant functions and $\hat{g}(t) = 0$ for $t > T_{\mathcal{K},\max}$, the integral $\int_0^\infty \hat{f}(x,t)\hat{g}(t)dt$ is a finite sum that can be easily computed.

## B Proof of Theorem 2: error of the Fleming-Viot estimator of the stationary probability, $\hat{p}_{\mathrm{FV}}^\pi$

Using the notation $f, g$ introduced in Appendix A and estimators $\hat{f}, \hat{g}$ given in (18) and (19) respectively, we have $p^\pi(x) = \int_0^\infty f(x,t)g(t)dt$ and $\hat{p}_{\mathrm{FV}}^\pi(x) = \int_0^\infty \hat{f}(x,t)\hat{g}(t)dt$, for states $x$ outside the absorption set $\mathcal{A}$. We are thus interested in bounding:

$$ \mathbb{E} \left| \int_0^\infty \hat{f}\hat{g}dt - \int_0^\infty fgdt \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}\hat{g}dt - \int_{0}^{\infty} fgdt \right| \le \mathbb{E} \left| \int_{0}^{\infty} (f - \hat{f})gdt \right| + \mathbb{E} \left| \int_{0}^{\infty} \hat{f}(\hat{g} - g)dt \right| \quad (20) $$

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 (5) which bounds the estimation error of $m(\cdot, \xi_t^\nu)$ defined in Appendix A, 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_{\mathrm{FV}}}{\sqrt{N}}, $$

for some $C_{\mathrm{FV}} > 0$. It follows that:

$$ \sup_{t \ge 0} \mathbb{E} \left| \hat{f}(x,t) - f(x,t) \right| \le \frac{C_{\mathrm{FV}}}{\sqrt{N}}, \quad \forall x \in \mathcal{A}^c. \qquad (21) $$

As was mentioned in subsection 3.2, 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_{\tilde{\partial}\mathcal{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 (21), we
obtain:

$$
\begin{align*}
\mathbb{E} \left| \int_0^\infty (f - \hat{f}) g dt \right| &\le \mathbb{E} \int_0^\infty |f - \hat{f}| g dt \\
&\le \int_0^\infty \mathbb{E} |f - \hat{f}| g dt \\
&\le \frac{C_{\text{FV}}}{\sqrt{N}} \int_0^\infty g dt = \frac{\mathbb{E}_{\vec{\partial}_A^c} T_K}{\mathbb{E}_{\vec{\partial}_A} T_A} \frac{C_{\text{FV}}}{\sqrt{N}},
\end{align*}
$$

where in the last line we also used the 'wedding cake decomposition', $\int_0^\infty \mathbb{P}_{\vec{\partial}_A^c}(T_K > t)dt = \mathbb{E}_{\vec{\partial}_A^c}T_K$.
Since $\hat{f}(x, t)$ is a probability, we get:

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

We thus wish to bound $\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_{\vec{\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_{\vec{\partial}_A^c}^\pi$. Since we start the simulation of $X_t^\pi$ for the purpose of estimating $g$ with distribution $p_{\vec{\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_{\vec{\partial}_A}^\pi$, it follows from renewal theory [2] 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}_{\vec{\partial}_A^c}(T_K > t)$ and $D_A = \mathbb{E}_{\vec{\partial}_A} T_A$, $D_K = \mathbb{E}_{\vec{\partial}_A^c} T_K$. We also denote by $\hat{N}_t$, $\hat{D}_A$ the estimators of $N_t$, $D_A$, that is $\hat{N}_t = \frac{1}{N} \sum_{i=1}^N \mathbf{1}_{T_{K,i}>t}$ and $\hat{D}_A = \hat{D}_E + \hat{D}_K$, with

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

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

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

We are interested in bounding:

$$
\mathbb{E} \int_{0}^{\infty} \left| \frac{\hat{N}_{t}}{\hat{D}_{A}} - \frac{N_{t}}{D_{A}} \right|
$$

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{\hat{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_{K,i} = \hat{D}_K$, the first term on the right hand side of the above bound is equal to $\mathbb{E} \left| \frac{\hat{D}_K}{\hat{D}_A} - \frac{\hat{D}_K}{D_A} \right|$. 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}_K}{\hat{D}_A} - \frac{\hat{D}_K}{D_A} \right| = \mathbb{E} \mathbf{1}_B \left| \frac{\hat{D}_K}{\hat{D}_A} - \frac{\hat{D}_K}{D_A} \right| + \mathbb{E} \mathbf{1}_{B^c} \left| \frac{\hat{D}_K}{\hat{D}_A} - \frac{\hat{D}_K}{D_A} \right| \quad (22)
$$
and bound each of the terms separately.

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

$$ \mathbb{E} \mathbf{1}_B \left| \frac{\hat{D}_K}{\hat{D}_A} - \frac{\hat{D}_K}{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 [28, Theorem 2.13], the random variable $\bar{T}_A$ is subexponential. Furthermore $\mathbb{P}(B) \le \mathbb{P}_{\bar{\partial}_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 [28][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 (22), 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}_K \left| \frac{1}{\hat{D}_A} - \frac{1}{D_A} \right| \le \frac{4}{D_A^2} \mathbb{E} \hat{D}_K \left| \hat{D}_A - D_A \right|. $$

Using the Cauchy-Schwartz inequality, we get:

$$ \begin{aligned} \mathbb{E} \hat{D}_K |\hat{D}_A - D| & \le (\mathbb{E} \hat{D}_K^2)^{1/2} (\mathbb{E} |\hat{D}_A - D_A|^2)^{1/2} \\ & \le \left( (\mathbb{E} \hat{D}_K)^2 + \mathrm{Var}(\hat{D}_K) \right)^{1/2} (\mathrm{Var}(\hat{D}_A))^{1/2} \\ & \le \left[ \frac{(\mathrm{Var}_{\bar{\partial}_A} T_E)^{1/2}}{\sqrt{M}} + \frac{(\mathrm{Var}_{\bar{\partial}_A} T_K)^{1/2}}{\sqrt{N}} \right] \left( (D_K)^2 + \frac{1}{M} \mathrm{Var}_{\bar{\partial}_A^c} (D_K) \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 \leq \mathbb{E}_{\bar{\partial}_{A}^{c}} T_{K} \left[ \frac{(\mathrm{Var}_{\bar{\partial}_{A}} T_{E})^{1/2}}{\sqrt{M}} + \frac{(\mathrm{Var}_{\bar{\partial}_{A}} T_{K})^{1/2}}{\sqrt{N}} \right] + \mathcal{O}\left(\frac{1}{M}\right). $$

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}{D_A} \mathbb{E} \int_0^\infty |\hat{N}_t - N_t| dt &= \frac{1}{D_A} \int_0^\infty \mathbb{E} |\hat{N}_t - N_t| dt \\ &\le \frac{1}{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} $$
where in the first inequality we use $EY \le \sqrt{EY^2}$ which follows from Cauchy-Schwartz inequality. We note, that $N_t = \mathbb{P}_{\bar{\partial}A^c}(T_K > t) = 1 - \mathbb{P}_{\bar{\partial}A^c}(T_K \le t) = 1 - F_K(t)$. Thus we have

$$ \frac{1}{D_A} \mathbb{E} \int_0^\infty |\hat{N}_t - N_t| dt \le \frac{1}{\sqrt{N} \mathbb{E}_{\bar{\partial}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^\pi$ 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 $\mathbb{P}_{\bar{\partial}A}(T_A > t) \le C_A \exp(-\lambda_A t)$ and $\mathbb{P}_{\bar{\partial}A^c}(T_K > t) \le C_K \exp(-\lambda_K t)$. Therefore all the moments and the integral above are finite.

# C Proof of the convergence of the FVRL algorithm to a local optimizer

## C.1 Proof of Proposition 4

Our optimisation control problem falls within the scope of Theorem 2.1 in [14, Section 5.2]. This theorem provides a convergence guarantee to a limit set of a projected ordinary differential equation (ODE), which includes the FVRL algorithm as a special case. To see this, we note that the $\theta$ update rule given by FVRL can be thought of as the discretisation of the projected ODE on the time variable $t$, $\dot{\theta}(t) = \nabla_{\theta}g(\pi_{\theta}) + z(t)$, where $g(\pi_{\theta})$ is defined in (6) and $z(t)$ is the projection term representing the minimum force needed to keep $\theta(t)$ in a pre-defined compact set $K \subset \mathbb{R}^d$ (for more details see [14, Section 4.3, "Projected ODE"]). This projected ODE is a special case of the more general case considered in Theorem 2.1, where the driving process defining the dynamics of $\theta(t)$ is a gradient function of $\theta(t)$ itself. Under this context, the limit set to which the trajectories of the projected ODE converge is actually called a stationary set. In turn, the stationary points in this set are easily seen to be local optimizers of the $g(\pi_{\theta})$ function. In fact, they are either stationary points of $g(\pi_{\theta})$ when inside the compact set $K$ where $z(t)$ is zero, or located at the boundary of the compact set, where $z(t) = -\nabla_{\theta}g(\pi_{\theta})$, making $\dot{\theta} = 0$.

We now prove that all the assumptions (1.1), (A2.1)-(A2.5) for Theorem 2.1 in [14] are satisfied by the FVRL algorithm, with which a biased estimate of the unknown gradient $\nabla_{\theta}g(\pi_{\theta})$ is obtained. For ease of notation we use thereafter $\nabla$ to mean $\nabla_{\theta}$.

**Assumptions (1.1) and (A2.4)** require that the learning rate satisfies the Robbins-Monro conditions for the convergence of stochastic approximation algorithms [23], $\sum_k \alpha_k = \infty$ and $\sum_k \alpha_k^2 < \infty$. This is directly satisfied by our choice of the learning rate scheduling as $\alpha_k = \alpha_0/k$ for some $\alpha_0 > 0$, as per item h) of the FVRL Algorithm 1.

**Assumption (A2.1)** requires finite second moments of the estimated gradient for all iterates $k$. This follows from:

$$ \sup_k \mathbb{E}^{\pi_{\theta_k}} \left\| \hat{\nabla} g(\pi_{\theta_k}) \right\|^2 = \sup_k \mathbb{E}^{\pi_{\theta_k}} \left[ \sum_{x \in S} \sum_{y \in S} \hat{p}_{FV}^{\pi_{\theta_k}}(x) \hat{p}_{FV}^{\pi_{\theta_k}}(y) \times \sum_{a \in \mathcal{U}} \sum_{b \in \mathcal{U}} \langle \hat{Q}^{\pi_{\theta_k}}(x, a) \nabla_{\pi_{\theta_k}}(a|x), \hat{Q}^{\pi_{\theta_k}}(y, b) \nabla_{\pi_{\theta_k}}(b|y) \rangle \right] \\ \le |S|^2 |\mathcal{U}|^2 \sup_k \mathbb{E}^{\pi_{\theta_k}} \left\{ \max_{z \in S, c \in \mathcal{U}} \left[ \hat{Q}^{\pi_{\theta_k}}(z, c) \|\nabla_{\pi_{\theta_k}}(c|z)\| \right]^2 \right\}. $$

Because rewards are by assumption bounded, $\hat{Q}^{\pi_{\theta_k}}(z, c)$ is bounded for all states $z$ and actions $c$, and since policy $\pi_{\theta}$ is assumed to be a continuously differentiable function of $\theta$ in a compact set $K \in \mathbb{R}^d$, $\|\nabla_{\pi_{\theta_k}}(c|z)\|$ is also bounded for all states $z$ and actions $c$, and this holds for all $k$. Hence the supremum over $k$ of the second moment of the estimated long-run expected reward gradient is finite.

**Assumption (A2.2)** requires that the expectation of the $\theta$-update driving process conditional to the past be decomposed as a measurable function of $\theta$ plus a random process $\beta_k$. In our case, this is automatically satisfied because
the $\theta$-update driving process, $\hat{\nabla}g(\pi_{\theta_k})$, is a function of $\theta_k$ and thus independent of its past values, $\hat{\nabla}g(\pi_{\theta_l}), l < k$. Hence, $\mathbb{E}[\hat{\nabla}g(\pi_{\theta_k})|\theta_0, \hat{\nabla}g(\pi_{\theta_l}), l < k] = \mathbb{E}[\hat{\nabla}g(\pi_{\theta_k})]$. The assumption then follows from the fact that $\mathbb{E}[\hat{\nabla}g(\pi_{\theta_k})] = \nabla g(\pi_{\theta_k}) + \beta_k$, where $\beta_k$ is a deterministic bias term.

**Assumption (A2.3)** requires that $\nabla g(\pi_{\theta})$ is continuous. This immediately follows from the expression for $\nabla g(\pi_{\theta})$ in (7), using the proposition's assumption that $\pi_{\theta}$ is a continuous differentiable function of $\theta$ in a compact set $\mathcal{K}$, that the state and action spaces are finite, and that the action value function $Q^{\pi_{\theta}}(x, a)$ is bounded because rewards are bounded.

**Assumption (A2.5)** requires that the bias term $\{\beta_k\}_{k \ge 0}$ of the gradient goes to zero fast enough, namely such that $\sum_k \alpha_k|\beta_k| < \infty$. If hyperparameters $M_k, N_k$ scale at least as fast as $k^\epsilon$ and hyperparameter $S_k$ scales at least as fast as $\log k^\epsilon$ for some $\epsilon > 0$, namely, $M_k = \lceil M_0 k^\epsilon \rceil, N_k = \lceil N_0 k^\epsilon \rceil, S_k = \lceil S_0 \log k^\epsilon \rceil$, then using Lemma 5, there exist constants $A, B, C, c > 0$ such that the bias $\beta_k$ is bounded by

$$ \begin{aligned} |\beta_k| & \le A \left( \frac{1}{\sqrt{M_0 k^\epsilon}} + \frac{1}{\sqrt{N_0 k^\epsilon}} \right) + B e^{-c S_0 \log k^\epsilon} + C e^{-c S_0 \log k^\epsilon} \left( \frac{1}{\sqrt{M_0 k^\epsilon}} + \frac{1}{\sqrt{N_0 k^\epsilon}} \right) \\ & = \left( \frac{A'}{k^{\epsilon/2}} + \frac{B}{k^{\epsilon c S_0}} + \frac{C'}{k^{\epsilon(1/2+cS_0)}} \right), \end{aligned} $$

where $A' = A(\frac{1}{\sqrt{M_0}} + \frac{1}{\sqrt{N_0}})$ and $C' = C(\frac{1}{\sqrt{M_0}} + \frac{1}{\sqrt{N_0}})$.

This yields the convergence of the $\sum_k \alpha_k|\beta_k|$ series, since by assumption $\alpha_k = \alpha_0/k$, and

$$ \sum_k \alpha_k |\beta_k| \le \alpha_0 \left( A' \sum_k \frac{1}{k^{1+\epsilon/2}} + B \sum_k \frac{1}{k^{1+\epsilon c S_0}} + C' \sum_k \frac{1}{k^{1+\epsilon(1/2+cS_0)}} \right) < \infty, $$

as all the exponents of $k$ are strictly larger than 1 for all $\epsilon, c > 0$.

**Remark 4.** If a more general learning rate scheduling of the FVRL algorithm were used, such as $\alpha_k = \alpha_0/k^\gamma$ for some $\frac{1}{2} < \gamma \le 1$ (still satisfying assumptions (1.1) and (A2.4)), Assumption (A2.5) would no longer be satisfied for all $\epsilon > 0$. The reason is that the exponents of $k$ in the two last series of the last expression in the proof would be equal to $\gamma + \epsilon c S_0$ and $\gamma + \epsilon(1/2 + cS_0)$, respectively. This would require knowledge of the constant $c$ in order to determine the minimum value of $\epsilon$ guaranteeing that those exponents are greater than 1, as required for the series in (A2.5) to converge.

## C.2 Proof of Lemma 5

Let $\eta^{\pi_{\theta}}(x) \doteq \nabla_{\theta} g_x(\pi_{\theta})$ and its estimator be denoted by $\hat{\eta}^{\pi_{\theta}}(x)$. Then, the bias $\beta$ of the gradient estimated by item d) of the FVRL Algorithm 1 is given by:

$$ \begin{aligned} \beta &= \mathbb{E}^{\pi_{\theta}} [\hat{\nabla}_{\theta} g(\pi_{\theta})] - \nabla_{\theta} g(\pi_{\theta}) = \sum_{x \in S'} \mathbb{E}^{\pi_{\theta}} [\hat{p}_{FV}^{\pi_{\theta}}(x)] \mathbb{E}^{\pi_{\theta}} [\hat{\eta}^{\pi_{\theta}}(x)] - \sum_{x \in S} p_{FV}^{\pi_{\theta}}(x) \eta^{\pi_{\theta}}(x) \\ &= \sum_{x \in S} (p^{\pi_{\theta}}(x) + b_p(x))(\eta^{\pi_{\theta}}(x) + b_{\eta}(x)) - \sum_{x \in S} p_{FV}^{\pi_{\theta}}(x) \eta^{\pi_{\theta}}(x) \\ &= \sum_{x \in S} b_p(x) \eta^{\pi_{\theta}}(x) + \sum_{x \in S} b_{\eta}(x) p^{\pi_{\theta}}(x) + \sum_{x \in S} b_p(x) b_{\eta}(x), \end{aligned} \quad (23) $$

where we have used the fact that the estimators $\hat{p}_{FV}^{\pi_{\theta}}(x)$ and $\hat{\eta}^{\pi_{\theta}}(x)$ are mutually independent (since they are estimated from independent simulations), and we have extended the sum over $S'$ to the sum over $S$, given that the gradient of the policy (intervening in $\hat{\eta}^{\pi_{\theta}}$) is zero outside $S'$ (by definition of $S'$ as the set of states where the policy gradient is not zero for at least one action –see initialization step of Algorithm 1).
We now bound each of the three bias terms in (23).

**Bound for the bias** $b_p(x)$ of $\hat{p}_{FV}^{\pi_{\theta}}(x)$: Given the number of particles $N$ and the number $M$ of return cycles to $\mathcal{A}$, by Theorem 2 the bias $b_p(x)$ of the stationary probability estimator $\hat{p}_{FV}^{\pi_{\theta}}(x)$ is upper bounded by:

$$|b_p(x)| = |\mathbb{E}[\hat{p}_{FV}^{\pi_{\theta}}(x)] - p^{\pi_{\theta}}(x)| \le \mathbb{E}|\hat{p}_{FV}^{\pi_{\theta}}(x) - p^{\pi_{\theta}}(x)| \le C_p \left( \frac{1}{\sqrt{M}} + \frac{1}{\sqrt{N}} \right),$$

for some $C_p > 0$.

**Bound for the bias** $b_{\eta}(x)$ of $\hat{\eta}^{\pi_{\theta}}(x)$: Using Lemma 3, the state contribution to the gradient, $\eta^{\pi_{\theta}}(x)$, can be estimated following the procedure described by items a) and b) for the FVRL Algorithm 1 on at most $S$ time steps, as

$$\hat{\eta}^{\pi_{\theta}}(x) = \frac{1}{V} \sum_{v=1}^{V} \sum_{a \in \mathcal{A}} \sum_{n=1}^{\min(\tau_v, S)} R_{v,n}^{(x,a)} \nabla \pi_{\theta}(a|x),$$

where $V$ is the number of replications run and $\tau_v$ is the value of the coupling time defined in Lemma 3 observed at replication $v$.

Using that all rewards $\{R_{v,n}^{(x,a)}\}_n$ for each replication $v$ are identically distributed, the expectation of $\hat{\eta}^{\pi_{\theta}}(x)$ can be written as:

$$\mathbb{E}^{\pi_{\theta}}[\hat{\eta}^{\pi_{\theta}}(x)] = \sum_{a \in \mathcal{A}} \mathbb{E}^{\pi_{\theta}} \left[ \sum_{n=1}^{\tau} R_n \mathbf{1}_{\{\tau \le S\}} + \sum_{n=1}^{S} R_n \mathbf{1}_{\{\tau > S\}} \middle| X_0 = x, A_0 = a \right] \nabla \pi_{\theta}(a|x), \quad (24)$$

where $\{R_n\}_n$ is the reward process and $\tau$ is the coupling time defined in Lemma 3.

Hence, using (24), the expression of $\eta^{\pi_{\theta}}(x)$ given by Lemma 3, and the assumption of bounded rewards, there exists a positive constant $C_r$ such that $|R_n| \le C_r$ and the bias of $\hat{\eta}^{\pi_{\theta}}(x)$ is bounded by:

$$\begin{aligned}
|b_{\eta}(x)| &= |\mathbb{E}^{\pi_{\theta}}[\hat{\eta}^{\pi_{\theta}}(x)] - \eta^{\pi_{\theta}}(x)| \\
&= \left| \sum_{a \in \mathcal{A}} \mathbb{E}^{\pi_{\theta}} \left[ \sum_{n=\tau+1}^{S} R_n \mathbf{1}_{\{\tau>S\}} \middle| X_0 = x, A_0 = a \right] \nabla \pi_{\theta}(a|x) \right| \\
&\le \sum_{a \in \mathcal{A}} \mathbb{E}^{\pi_{\theta}} \left[ \sum_{n=1}^{S} |R_n| \mathbf{1}_{\{\tau>S\}} \middle| X_0 = x, A_0 = a \right] \| \nabla \pi_{\theta}(a|x) \| \\
&\le C_r S \sum_{a \in \mathcal{A}} \mathbb{P}(\tau > S | X_0 = x) \| \nabla \pi_{\theta}(a|x) \|,
\end{aligned}$$

Since the state space is finite and the MDP is assumed to be ergodic for all policies in the parameterised policy class, $\mathbb{P}(\tau > S | X_0 = x) \le e^{-c_0 S}$ for some positive constant $c_0$ (Doeblin's coupling theorem for ergodic Markov chains, [5, Section 4.3.1]). Finally, because $\theta$ takes values in a compact set and the policy is assumed to be a continuous differentiable parameterisation by $\theta$, the norm of the policy gradient is bounded by some $C_{\pi} > 0$, and this yields the following bound for the bias:

$$|b_{\eta}(x)| \le C_r C_{\pi} |\mathcal{U}| S e^{-c_0 S} \le C_{\eta} e^{-cS} \text{ for some } C_{\eta} > 0, c > 0,$$

i.e. the bias of $\hat{\eta}^{\pi_{\theta}}(x)$ tends to zero exponentially fast with an increasing $S$.

Finally, using that the state and action spaces are finite, that $\pi_{\theta}$ is a continuous differentiable function of $\theta$ in a compact set $\mathcal{K}$, and that the action value function $Q^{\pi_{\theta}}(x, a)$ is bounded because rewards are bounded by assumption,
$\eta^{\pi\theta}(x)$ is upper bounded by some $B_{\eta} > 0$. Thus, the bias term $\beta$ can be bounded as:

$$ \begin{aligned} |\beta| &\le \sum_{x \in S} |b_p(x)\eta^{\pi\theta}(x)| + \sum_{x \in S} |b_{\eta}(x)|p^{\pi\theta}(x) + \sum_{x \in S} |b_p(x)||b_{\eta}(x)| \\ &\le |S| \left( C_p B_{\eta} \left( \frac{1}{\sqrt{M}} + \frac{1}{\sqrt{N}} \right) + C_{\eta}e^{-cS} + C_p C_{\eta}e^{-cS} \left( \frac{1}{\sqrt{M}} + \frac{1}{\sqrt{N}} \right) \right). \end{aligned} $$

The lemma then follows by defining the positive constants $A = |S|C_pB_{\eta}$, $B = |S|C_{\eta}$, $C = |S|C_pC_{\eta}$.

## D Heuristics for the choice of the simulation parameters for FVEE

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 hyperparameters $J$, $N$ and $T$ of the Fleming-Viot estimation procedure of the blocking probability described in Section 4.1.1, so that they can be used as an initial guideline for the choice of these parameters in a more general setting such as multidimensional problems.

According to Theorem 2, the error of the Fleming-Viot estimator is controlled by two independent quantities: $N$, the number of Fleming-Viot particles, and $M$, the number of return cycles to $A$ observed in the simulation used to estimate the denominator, $\mathbb{E}(T_A)$ (where the expectation subscript has been dropped for conciseness). While $N$ is a hyperparameter of the Fleming-Viot simulation, $M$ is a random variable that depends on both $J$ and $T$. As stated in Section 4.1.1, $J$ is problem dependent (e.g. it should be chosen as a function of the state visit frequency), which in practice makes $M$ a function of hyperparameter $T$. Clearly, for a fixed value of $J$, a larger number of cycles $M$ is expected for a larger number of arrivals $T$ allowed in the $\mathbb{E}(T_A)$ simulation, since more time is allowed for the simulation.

Therefore, in what follows our analysis focuses on how the error of the Fleming-Viot estimator varies as $N$ and $T$ vary, while keeping $J$ fixed to a predefined quantity. The value of $J$ is chosen so that a few meaningful $(N, T)$ value pairs can be considered for the estimation of the blocking probability which involve simulations that take a reasonable amount of time to run. The criterion that $J$ defines the number of states in the $M/M/1/K$ system having a stationary probability larger than 0.5% proved reasonable. For a blocking size $K = 20$ and system load $\rho = 0.7$, this criterion yields $J = 12$.

We focus the estimator error analysis on the two quantities (out of the three involved in expression (10) of Section 4.1.1) 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$; the remaining quantity, $\mathbb{P}_J(T_K > t)$, does not strongly depend on $J$ (see Appendix B)$^{(7)}$. 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. As mentioned above, we now invert our reasoning: we fix $J$, and set the values of $N$ and $T$ required to approximately satisfy predefined expected relative errors in these estimators. We then run simulations using each $(N, T)$ value pair dictated by the relative error pairs considered, and study statistics on the relative error of the Fleming-Viot estimator of the blocking probability to understand how it varies as a function 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)$: For fixed $J$, the expected relative error of $\hat{\phi}_t^J(K)$ is approximately inversely proportional to $\sqrt{N}$.

*Derivation:* 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

<footnote>$^7$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>
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_J$ of a single-server queue system with capacity equal to $\lceil \frac{K-J+1}{2} \rceil$ ($\lceil \cdot \rceil$ is the ceiling operator), namely $q_J = \frac{(1-\rho)\rho^{K_J}}{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_J$) random variable as $\epsilon_\phi \approx \sqrt{(1-q_J)/Nq_J}$. Thus, given $q_J$, the value of $N$ to approximately satisfy a desired expected relative error $\epsilon_\phi$ for $\hat{\phi}_t^J(K)$ is obtained as $N \approx \lceil \frac{1-q_J}{q_J\epsilon_\phi^2} \rceil$, which in turn can be approximated as $N \approx \lceil \frac{1}{q_J\epsilon_\phi^2} \rceil$, if $q_J \ll 1$, i.e. when $J$ is sufficiently small compared to $K$. We finally solve for $\epsilon_\phi$ to get its inverse relationship with $\sqrt{N}$ stated above.

2. **Relative error of $\hat{E}_{J-1}T_A$:** For fixed $J$, the expected relative error of $\hat{E}_{J-1}T_A$ is approximately inversely proportional to $\sqrt{T}$.

*Derivation:* 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 $E_{J-1}(T_{J-1}) = \frac{1}{\lambda_{p_J}(1+\rho^{-1})}$ [5], where $p_J$ is the stationary probability of the state $x=J-1$, i.e. $p_J = \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 $E_{J-1}T_A = 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 = ME_{J-1}T_A = ME_{J-1}(T_{J-1})(1+\rho^{-1})$, i.e. $t_S = \frac{M}{\lambda_{p_J}}$. 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 \approx M/p_J$ arrival events. On the other hand, given $M$ return cycles to $\mathcal{A}$, the standard error of the moment estimator of $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) \approx E_{J-1}T_A$, hence the relative error of the estimator of $E_{J-1}T_A$, $\epsilon_{ET}$, is of the order of $1/\sqrt{M}$, which makes $M \approx 1/\epsilon_{ET}^2$. Thus, given $p_J$, the value of $T$ to approximately satisfy a desired expected relative error $\epsilon_{ET}$ for $\hat{E}_{J-1}T_A$ is obtained as $T \approx \lceil \frac{1}{p_J\epsilon_{ET}^2} \rceil$. We finally solve for $\epsilon_{ET}$ to get its inverse relationship with $\sqrt{T}$ stated above.

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 $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{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^8$, the value of $N$ as a function of the stationary probability $q_J$ 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_J$ is dominated by an increasing exponential function of $J$, i.e. $\sim O(\rho^{-J})$ which, remarkably, does not depend on $K$. Thus, as mentioned in remark 3, the closer $J$ is 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$.

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 of Section 4.1.1. Thus, from the computational

<footnote n="8">⁸$K$ is considered sufficiently large in this context when $\rho^{K-J}$ can be neglected w.r.t. 1.</footnote>
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 error of the FV estimator of the blocking probability.

We completed the study of the appropriate choice of $N$ and $T$ by analyzing the impact of the errors of estimating $\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$, which corresponds to choosing the states with stationary probability larger than 0.5%. The results of these experiments are shown in the heatmap of Figure 7 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 $T$ of arrival events for estimating $\mathbb{E}(T_A)$ 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\%^9$; 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\%^{10}$.

<footnote n="9">⁹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="10">¹⁰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_9.png" /><figcaption>Figure 7: Heatmap showing the impact of different combinations of the expected relative errors of $\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 $\mathcal{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>

## E FVEE and FVRL algorithms for threshold-type policies

This section presents the two algorithms used throughout this work to apply the Fleming-Viot methodology to network systems under threshold-type admission control policies, in order to (i) estimate the rejection probability and the expected rejection cost using the Fleming-Viot expectation estimator (FVEE Algorithm 2), and (ii) learn locally optimal blocking sizes of threshold-type admission control policies using the Fleming-Viot reinforcement learning method (FVRL Algorithm 3). Note that the latter is essentially the FVRL Algorithm 1 presented in the paper's body, adapted to piecewise-linear admission control policies defined by (12).
**Algorithm 2** FVEE 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* queue system is a particular case (with *I* = 1, *R* = 1).

**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 assumed known for the estimation problem.

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^{\vec{\partial}_{\mathcal{A}^c}}(K)$ and $\mathbb{P}_{\vec{\partial}_{\mathcal{A}^c}}(T_K > t)$ in (3).

F. *T*: number of arrival events (incoming jobs) for estimating $\mathbb{E}(T_A)$, which has a direct impact on the number *M* of return cycles to *A* used to estimate $\mathbb{E}_{\vec{\partial}_{\mathcal{A}}} 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}_{\vec{\partial}_{\mathcal{A}}} 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}_{\vec{\partial}A} 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 $\vec{\partial}A = \{\mathbf{x} : x_i = J_i - 1 \text{ for some } i = 1, \dots, I\}$.

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

(c) Record the states at which the process enters $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 $A$ after the burn-in period.

(e) If $M > M_0$, estimate $\mathbb{E}^{\vec{\partial}A} 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 $A^c$, $\hat{p}_{\vec{\partial}A^c}^\pi(\mathbf{x})$ as $\frac{1}{M} \sum_{k=1}^{\overline{M}} \mathbf{1}_{\mathbf{x}(k)=\mathbf{x}}$, $\forall \mathbf{x} \in \vec{\partial}A^c$, where $\overline{M}$ is the number of observed entrances to $A^c$ at states $\mathbf{x}(k)$.

2. Estimation of $\mathbb{P}_{\vec{\partial}A^c}(T_K > t)$ and $\phi_t^{\vec{\partial}A^c}(\mathbf{x}_C)$, where $\mathbf{x}_C$ is any state in the set of blocking states $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 $\vec{\partial}A^c$ following $\hat{p}_{\vec{\partial}A^c}^\pi$. Record the observed values of the first killing times $T_K$ at which each of the $N$ particles enters $A$.

(b) Every time one of the particles is killed at $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}_{\vec{\partial}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^{\vec{\partial}A^c}(\mathbf{x}_C)$ as the proportion of particles that are at state $\mathbf{X}_t^\pi = \mathbf{x}_C$.

3. Estimation of the expected rejection cost:

Let $\eta(\mathbf{x})$ be the expected cost accrued at rejection over all possible arriving job classes when the system is in state $\mathbf{x} \in C$, namely: $\eta(\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\}}$], where $\Lambda = \sum_{i=1}^I \lambda_i$ is the total job arrival rate. Estimate the expected rejection cost, $\mathbb{E}(\eta)$, using (1), with the stationary probability estimated using (4). The integral giving $\hat{p}_{FV}^\pi(\mathbf{x})$ is easily computed as a finite sum of the piecewise constant function of time, resulting from the product $\hat{\mathbb{P}}_{\vec{\partial}A^c}(T_K > t) \hat{\phi}_t^{\vec{\partial}A^c}(\mathbf{x})$, which is 0 for $t > \max(\{T_{K,k}\}_{k=1...N})$.
**Algorithm 3** FVRL algorithm for piecewise-linear admission control policies to find blocking sizes $K_i$ that locally minimize

the expected cost defined in (15).

The system is an $M/M/I/R$ loss network serving $I$ different job classes with $R$ servers, of which the $M/M/1$ queue system is a particular case (with $I = 1, R = 1$).

**Data:**

A. A simulator or emulator of the loss network or queue system to control (in our case, a simulator was obtained by defining the arrival rate $\lambda_i$ and the service rate $\mu_i$ of each job class $i = 1, \dots, I$).

B. Cost of rejecting an arriving job class, $C_i$, $i = 1, \dots, I$.

C. $\pi_{\theta_i}$ ("accept" $|x_i|$): the job acceptance policy for an arriving job of class $i$ parameterised by the positive real-valued $\theta_i$, as defined in expression (12).

D. $\theta = (\theta_i)_{i=1,...,I}$: positive non-integral initial values of the parameter to optimise, from where $K_i$ can be obtained as $K_i = \text{ceiling}(\theta_i + 1)$.

E. $F_i$: the size of the absorption set $\mathcal{A}$ in dimension $i$ as a fraction of $K_i$.

F. $L$: number of learning steps, i.e. the number of times an update of $\theta$ will be computed by the gradient-based algorithm.

G. $N, T$: respectively, the number of FV particles and the number of arrival events for estimating $\mathbb{E}(T_A)$ used to estimate the blocking probability with the FVEE Algorithm 2.

H. $V, S$: respectively, the number of replications and the maximum number of time steps allowed to estimate the function $\nabla_{\theta}g_x(\pi_{\theta})$ defined in (8).

I. $\alpha_0 > 0$: initial learning rate used in the $\theta$ update (e.g. $\alpha_0 = 1$).

**Result:** An estimate of blocking sizes $\hat{K}_i^*$, $i = 1, \dots, I$, that locally minimize the expected rejection cost defined in (15).

**Steps:**

1. Compute the deterministic blocking sizes $K_i = \text{ceiling}(\theta_i + 1)$, $i = 1, \dots, I$. Set $J_i$, the size of the absorption set $\mathcal{A}$ in dimension $i$ as $J_i = [F_i K_i]$.

2. Estimate the stationary probabilities $p^{\pi_{\theta}}(\mathbf{x})$ for each $\mathbf{x}$ such that $x_i = K_i - 1$ and $\sum_{j=1}^{I} x_j < R$ using the FVEE Algorithm 2.

3. For each $\mathbf{x}$ considered in the previous step simulate $V$ times two independent $\{X_n\}_{n \in \mathbb{N}}$ processes in parallel under the current policy $\pi_{\theta}$, one with initial action $a = 0$ (reject first incoming job) and the other with initial action $a = 1$ (accept first incoming job). Collect the rewards observed by each copy, respectively $\{R_{v,n}^{(x,1)}\}_{n \ge 1}$ and $\{R_{v,n}^{(x,0)}\}_{n \ge 1}$ and stop the simulation when the two processes meet for the first time, or when each takes $S$ steps. Use Lemma 3 to compute the moment estimate of the $\nabla_{\theta}g_x(\pi_{\theta})$ function defined in (8), which in this piecewise-linear policy case, as per (16), reduces to: $\hat{\nabla}_{\theta}g_x(\pi_{\theta}) = \frac{1}{V} \sum_{v=1}^{V} \sum_{n=1}^{\min(\tau_v, S)} (R_{v,n}^{(x,1)} - R_{v,n}^{(x,0)})$, where $\tau_v$ is the first time the two processes meet at replication $v$.

4. Use the estimates of steps 2 and 3 to estimate the gradient of the long-run expected cost, $\hat{\nabla}g(\pi_{\theta}) = \{\frac{\partial g(\pi_{\theta})}{\partial \theta_i}\}_{i=1,...,I}$, using (16).

5. Update $\theta$ using gradient descent: $\theta \leftarrow \theta - \frac{\alpha_0}{k} \hat{\nabla}g(\pi_{\theta})$, where $k$ is the current learning step. Bound $\theta_i$ ($i = 1, \dots, I$) to the valid range, if applicable.

6. Repeat steps (1)-(5) until the number of learning steps $L$ is reached.

7. Use the final value $\theta^L$ to compute the final blocking sizes estimated by the algorithm, $\hat{K}_i^*$, as $\text{round}(\theta_i^L) + 1$, $i = 1, \dots, I$.