2. Related literature

Here we only review papers that are most closely related to our study and refer interested readers to the classic text by Gittins et al. (Reference Gittins, Glazebrook and Weber2011) for a systematic and comprehensive treatment of MABPs and to the recent paper by Verloop (Reference Verloop2016) for other related references.

Existing literature on the MABP and RMABP often assumes a stationary activation budget (i.e. the maximum number of bandits that can be activated per period is independent of time) and only two possible actions per bandit (Verloop (Reference Verloop2016)). To the best of our knowledge, a nonstationary activation budget has only been considered in Cohen et al. (Reference Cohen, Zhao and Scaglione2014). The contrast between Cohen et al. (Reference Cohen, Zhao and Scaglione2014) and our work is that we allow for an arbitrary finite number of actions and states, and we focus on developing an asymptotically optimal policy instead of deriving sufficient conditions for optimality of a certain policy form.

RMABPs with multiple actions are often referred to as superprocesses and were first considered for the MABP in the setting where only one bandit can be activated per period (Whittle (Reference Whittle1980)). For this setting, it was shown that, under the condition that each state has a dominant action, there is an optimal policy that is indexable. A less strict condition is given in Gittins et al. (Reference Gittins, Glazebrook and Weber2011), but still has the same purpose of ensuring that there is an optimal policy that is indexable. Multiple actions were also considered in Verloop (Reference Verloop2016) and generalize the concept of dominant action to their setting. Our proposed RAC makes no such restriction to a class of dominant action policies.

Verloop (Reference Verloop2016) analyzed multi-class restless bandits with a long-run average cost criteria. A class of priority policies that do not require indexability are introduced, and their asymptotic optimality is analyzed for a fixed and dynamic population of bandits that can arrive and depart from the system. Conditions for asymptotic optimality require that the differential equation that describes the fluid model associated with the RMABP have a globally attracting equilibrium point, similar to the conditions needed in Weber and Weiss (Reference Weber and Weiss1990). This condition guarantees that the equilibrium point induces a priority policy and that the process under this priority policy converges to the equilibrium point independent of the initial conditions. To guarantee the condition holds, the family of processes that scales to the fluid model must each have a unique invariant probability distribution with finite first moment, and the collection of these unique invariant probability distributions must be tight and uniformly integrable. For a fixed population of bandits, these conditions are satisfied when the generated Markov process is unichain (i.e. the state space for the Markov chain has a single recurrent class C and there is no other absorbing set disjoint from C), so that the resulting Markov chain has a unique equilibrium distribution. For a dynamic population of bandits, they are satisfied when the generated Markov process is irreducible and state 0 (i.e. the ‘empty’ state) is positive recurrent for any bandit that is never activated, i.e. inactivated bandits eventually leave the system.

Our dynamic population framework differs from that in Verloop (Reference Verloop2016) in four ways. First, we allow time-varying constraints on the number of active bandits per period (again, our results can also be extended to the setting with nonstationary transition matrices and nonstationary costs). Second, we consider a finite-horizon setting, so we do not require a global attractor condition. Third, we make no assumptions about the transition probability matrices or the underlying dynamics generated by the policies of consideration, such as being unichain. Fourth, we propose a new type of policy that is neither a priority policy nor an index policy, but is still asymptotically optimal under certain relatively mild conditions.

Kelly (Reference Kelly1981) considered an expected long-run average reward criterion for a finite population bandit. Each arm generates a sequence of Bernoulli random variables whose parameters are themselves random variables, and are independent with common distribution satisfying a regularity condition. He showed that, as the discount factor approaches 1, the optimal policy converges to that which pulls the arm(s) that has(ve) incurred the least number of ‘failures’. Moreover, as the discount factor approaches 1, the performance of the optimal discounted expected reward policy as judged with respect to the expected average reward criterion approaches the best possible.

Our work is also related to the literature that uses LP relaxation to approximate RMABPs. Bertsimas and Nino-Mora (Reference Bertsimas and Niño-Mora2000) were the first to consider a sequence of LP relaxations to obtain a primal-dual index policy for the infinite-horizon discounted-cost RMABP. Le Ny et al. (Reference Le Ny, Dahleh and Feron2008) extended the work of Bertsimas and Nino-Mora (Reference Bertsimas and Niño-Mora2000) to include switching times/costs between activating bandits. Both Bertsimas and Nino-Mora (Reference Bertsimas and Niño-Mora2000) and Le Ny et al. (Reference Le Ny, Dahleh and Feron2008) considered LP relaxations of the dynamic program formulation of infinite-horizon RMABPs. An alternative LP relaxation can be derived by considering a fluid model of the RMABP that only takes into account mean drifts of bandit dynamics (Weber and Weiss (Reference Weber and Weiss1990) and Verloop (Reference Verloop2016)). This fluid approach is also called the certainty equivalent approach in the operations research literature and is closest to the LP formulation in this paper.

Finite-horizon MABPs have been studied in Robbins (Reference Robbins1952) and Bradt et al. (Reference Bradt, Johnson and Karlin1956), who focused on the case where engaging a project corresponds to sampling from a Bernoulli population with unknown success probability and the objective is to maximize the expected number of successes over a finite number of plays. We refer interested readers to the monograph by Berry and Fristedt (Reference Berry and Fristedt1985) for additional references on finite-horizon MABPs. More recently, Caro and Gallien (Reference Caro and Gallien2007), considered the setting motivated by a dynamic assortment problem in the fashion retail industry. Nino-Mora (2011) (see also the references therein) considered a class of finite-horizon discrete-state bandit problems whose optimal policy is known to be of index type (i.e. the counterpart of the Gittins index for a finite-horizon discrete-state bandit), and proposed an efficient and exact algorithm to compute the index.

Finite horizon RMABPs have been less studied. To the best of our knowledge, the work by Hu and Frazier (Reference Hu and Frazier2017) is the closest to the current study. Hu and Frazier (Reference Hu and Frazier2017) considered the RMABP with a finite number of independent and identically distributed bandits which can be either activated or not, except over a finite horizon and with a fixed and stationary number M of multiple pulls per period. They proposed an index policy for this model and showed that it is asymptotically optimal. Their asymptotically optimal index policy was obtained in three main steps. First, they constructed a Lagrangian relaxation of the original RMABP, which provides an upper bound to the optimal value of the RMABP, and computed Lagrange multipliers that minimize the value of the RMABP relaxation using subgradient descent. Second, they decomposed the relaxation into smaller Markov decision processes (MDPs), one for each bandit, and each of which can be optimally solved using backward induction. For each smaller MDP, they obtained an ‘index’ for each time period. They then constructed an index policy for the original RMABP by activating the M largest indices corresponding to the smaller MDPs in each time period. One potential issue is that there can be ties between indices so that it may not be possible to choose exactly M indices. As a result, the third step is a tie-breaking rule that activates the remaining subprocesses between tied states according to a probability distribution induced by an optimal policy that is obtained by solving a linear program with respect to occupation measures and that depends on the Lagrange multipliers.

In addition to our population and multi-action setting, there are three other components of our approach that distinguish it from that of Hu and Frazier (Reference Hu and Frazier2017) and which make their approach not directly applicable to ours. First is the objective. We introduce a set of undesirable states, which are the states in which a bandit will incur a high penalty cost of at the end of the horizon. In particular, a cost is incurred for each excess bandit that ends up in an undesirable state immediately after the end of the horizon and we allow for at most a fixed number of bandits to be in any of the undesirable states at this time without being penalized. The dependency between the bandits introduced by this terminal cost implies that the Lagrangian relaxation and the subsequent decomposition into sub-MDPs in Hu and Frazier (Reference Hu and Frazier2017) is not directly applicable. The second difference is that we allow for the number of bandits to be activated in each period to be upper bounded by a time-dependent constant. Because the latter constraint is nonbinding, the Lagrangian relaxation of Hu and Frazier is again not directly applicable. Finally, the index policy in Hu and Frazier (Reference Hu and Frazier2017) is based on aggregating ‘indices’ from smaller MDPs based on the Lagrangian relaxation of the original RMABP and on a tie-breaking rule, akin to our formulation of our stochastic control problem, based on an LP formulation for state occupation measures and the Lagrange multipliers. By contrast, ours is based on an LP relaxation where the random variables in the original RMABP are replaced by their means.

Another potential distinction from Hu and Frazier (Reference Hu and Frazier2017) is that our analysis can be extended to include nonstationary transition probability matrices and costs and state-dependent terminal costs. It can also be generalized to the setting where we have several constraints for several subsets of bandits or one constraint on the total number of activated bandits over the entire finite-time horizon. The proposed approach would also work for the multi-action setting where each action contributes a different weight to the budget. To the best of our knowledge, we are the first to analyze finite-horizon RMABPs with a nonstationary activation budget and asymptotically optimal policies that are nonindexable.

3. Fixed population model

We now describe the setting of a fixed population model where the number of bandits in the system remains constant throughout the horizon (i.e. no new bandits arrive and no existing bandits depart); this setting will be assumed throughout Sections 3 to 5. We consider a finite-horizon, discrete-time model, where time $t \in \{1, \dots, T + 1 \}$ with $T \lt\infty$. Let $\mathbb{J} = \{j\colon 1 \le j \le J\}$ with $J \lt \infty$ denoting the set of feasible states, and let $\mathbb{A} = \{a\colon 0 \le a \le A\}$ with $A \lt \infty$ denoting the set of feasible actions. Our results and analysis will still apply when each state has its own set of feasible actions; so, without loss of generality, we will simply assume that all states share the same set of feasible actions $\mathbb{A}$. We also introduce a set of states $\mathbb{U} \subseteq \mathbb{J}$, called the undesirable states, which are the states in which a bandit will incur a high penalty cost at the end of the horizon. Undesirable states can represent machine failure in maintenance and replacement models, customers that did not complete service before the end of the horizon or abandoned the system before completing service in queueing scheduling applications, or deteriorated/critical health conditions in healthcare applications. We refer to action $a = 0$ as no action (or no treatment) and any action $a \gt 0$ as a proper action (or proper treatment). Moreover, using the standard terminology from the RMABP literature, we will call a bandit active when a proper action is applied to it and passive otherwise. We assume that the state of each bandit transitions from state i to state j under action a according to a probability $p^a_{ij}$ and that the maximum number of active bandits at time t is $b_t \gt 0$, which we call the activation budget. Two types of costs can be incurred by the system. First, a cost $c^a_j$ is incurred each time action a is applied to a bandit in state j. Second, a cost $\phi$ is incurred for each excess bandit that ends up in an undesirable state at time $T+1$; we allow for at most $m \ge 0$ bandits to be in any of the undesirable states at time $T+1$ without being penalized. Finally, we use $n_{{\rm{tot}}}$ to denote the total number of bandits in the system. For ease of exposition, we have used a time-independent transition matrix and cost function, but remark that our analysis can be applied to the setting with a time-dependent transition matrix and cost function.

The decision-making scenario is as follows. At the beginning of time t, the decision maker views the state of each bandit in the system and then decides the action/treatment to be applied to each bandit. After receiving the treatment, each bandit incurs a cost and its state transitions to a potentially new state at the beginning of time $t+1$. The objective of the decision maker is to minimize the expected total treatment and penalty costs, or simply, the expected total costs.

Let $\Pi$ denote the set of all nonanticipating controls, and let $\pi$ denote a feasible control in $\Pi$. For a given control $\pi$, since all bandits share identical transition dynamics and costs, we do need to keep track of the individual bandits; instead, it is sufficient to simply keep track of the number of bandits in state j that receive treatment a at time t, which we denote by $X^{\pi, a}_j(t)$, as well as the number of bandits in state i that receive treatment a at time t whose states transition to state j at time $t+1$, which we denote by $\smash{Y^{\pi,a}_{ij}(t)}$. Let $n_{j}$ denote the number of bandits in state j at time $t=1,$ and let $V^{S}$ denote the expected total costs under an optimal control $\pi^*$. Then

(1a)\begin{equation}V^{S} = \min_{\pi \in \Pi} \mathbb{E}^{\pi} \bigg[\sum_{t=1}^T\sum_{a=0}^A \sum_{j=1}^J c^a_j \cdot X^{\pi, a}_j(t) + \phi \cdot \bigg(\sum_{a=0}^A \sum_{j \in \mathbb{U}} X^{\pi,a}_j(T+1) - m \bigg)^+\bigg]\end{equation}

such that

(1b)\begin{equation}\sum_{a=0}^A X^{\pi, a}_j(t) = \sum_{a=0}^A\sum_{i=1}^J Y^{\pi, a}_{ij}(t-1) \quad\text{for all } j\ge 1, \, t \ge2, \end{equation}

(1c)\begin{equation}\sum_{a=0}^A X^{\pi, a}_j(1) = n_{j} \quad\text{for all } j \ge 1,\end{equation}

(1d)\begin{equation}\sum_{a=1}^A \sum_{j=1}^J X^{\pi, a}_j(t) \le b_t \quad\text{for all } t\ge 1 ,\end{equation}

where the constraints hold almost surely (i.e. with probability 1).

We remark that the seemingly simple stochastic control model (1) is in fact surprisingly general. As alluded to in Section 1, it can be used for machine maintenance and replacement problems or capacity management and/or resource allocation in healthcare, among others. Additionally, this model can be modified to include dynamic populations (i.e. ‘birth’ and ‘death’ RMABP) by including nonstationary random arrival processes, nonstationary random service capacity (e.g. nonstationary budget), random service completions, and random abandonments. Lastly, we remark that the formulation (1) includes the budget constraints explicitly, which are typically included in the definition of the class $\Pi$ of admissible controls. In what follows, and for ease of notation, whenever it is clear from the context which policy is being used, we will suppress the notational dependency on $\pi$.

In theory, the stochastic control problem (1) can be solved using the standard DP algorithm. Unfortunately, this approach suffers from the well-known curse of dimensionality. Thus, rather than solving (1) exactly, we instead construct and analyze a provably near-optimal heuristic control. To describe our heuristic control in Section 4, we need to introduce the following family of linear programs indexed by $\boldsymbol{\epsilon} \,:\!= ({\epsilon_1, \epsilon_2,\dots, \epsilon_T) \geq \mathbf{0} = (0, 0, \dots, 0)$:

(2a)\begin{equation}V^D(\boldsymbol{\epsilon}) = \min_{x,z} \sum_{t = 1}^T \sum_{a=0}^A \sum_{j =1}^{J} c_{j}^a \cdot x_{j}^a(t, \boldsymbol{\epsilon}) + \phi \cdot z(\boldsymbol{\epsilon})\end{equation}

such that

(2b)\begin{equation}\sum_{a=0}^{A} x_{j}^a(t, \boldsymbol{\epsilon}) =\sum_{a=0}^A \sum_{i=1}^{J} x_{i}^a(t-1, \boldsymbol{\epsilon}) \cdot p^a_{ij}\quad\text{for all } j\ge 1 , \, t \ge 2, \end{equation}

(2c)\begin{equation}\sum_{a = 0}^A x_{j}^{a} (1, \boldsymbol{\epsilon}) = n_{j} \quad\text{for all }j \ge 1, \end{equation}

(2d)\begin{equation}\sum_{a=1}^A \sum_{j=1}^J x^a_j(t, \boldsymbol{\epsilon}) \le b_t - \epsilon_t\quad\text{for all } t \ge 1, \end{equation}

(2e)\begin{equation}z(\boldsymbol{\epsilon}) \ge \sum_{j \in \mathbb{U}} \sum_{a=0}^A \sum_{i=1}^J x^a_i(T,\boldsymbol{\epsilon}) \cdot p^a_{ij} - m, \end{equation}

(2f)\begin{equation}z(\boldsymbol{\epsilon}), \, x^a_j(t, \boldsymbol{\epsilon}) \ge 0 \quad\text{for all } a\ge 0, \, j \ge 1, \, t \ge 1. \end{equation}

One can arrive at linear program (2) by replacing all the random variables in (1) with their expected values, replacing the original activation budget $b_t$ in (1) with $b_t -\epsilon_t$, and introducing a new variable $z(\boldsymbol{\epsilon})$ to capture the number of excess bandits in undesirable states at the end of the time horizon. In the special case when $\boldsymbol{\epsilon} = \mathbf{0}$, linear program (2) is simply the deterministic relaxation of the stochastic control problem (1). Otherwise, when $\boldsymbol{\epsilon}> \mathbf{0}$, we interpret vector $\boldsymbol{\epsilon}$ as a buffer for making sure that the activation budgets are not exceeded when a heuristic control constructed using the solution of (2) is implemented in the original stochastic system.

We end this section with the following result.

Lemma 1. It holds that $V^D(\mathbf{0}) \le V^S$.

Proof. Let $\pi^*$ denote an optimal policy for our original stochastic control problem (1), and let $\smash{x^a_j(t,\mathbf{0}) =\mathbb{E}[X^{\pi^* a}_j(t)]}$ for all a, j, t. By Jensen’s inequality,

\begin{equation}\mathbb{E}\bigg[\bigg(\sum_{a=0}^A \sum_{j \in U} X^{\pi^* a}_j(T+1) -m\bigg)^+\bigg] \ge \sum_{a=0}^A \sum_{j\in U} x^a_j(T+1,\mathbf{0}) - m.\end{equation}

Thus, let

\begin{align*}z(\mathbf{0}) =\begin{cases}\sum_{a=0}^A \sum_{j\in U} x^a_j(T+1,\mathbf{0}) - m & \text{if$\sum_{a=0}^A \sum_{j\in U} x^a_j(T+1,\mathbf{0}) \geq m$},\\0 & \text{otherwise.}\end{cases}\end{align*}

It follows that $\smash{x^a_j(t,\mathbf{0})}$ for all a, j, t and $z(\mathbf{0})$ provide a feasible solution to (2) with $\boldsymbol{\epsilon} \,:\!= (\epsilon_1, \epsilon_2, \dots, \epsilon_T) =\mathbf{0} = (0, 0, \dots, 0)$ having total costs in (2) smaller than the total expected costs in (1) under $\pi^*$. This yields $V^D(\mathbf{0}) \le V^S$.□

Lemma 1 allows us to use $V^D(\mathbf{0})$ as a proxy for $V^S$ and study the performance of any feasible control $\pi$ by analyzing the difference between $V^{\pi}$ and $V^D(\mathbf{0})$, as claimed.

4. Randomized activation control

We now discuss our heuristic control, which we call randomized activation control (RAC). To do this, we first need to introduce some notation. For a given $\boldsymbol{\epsilon} \geq \mathbf{0}$, let $x^{*,a}_j(t,\boldsymbol{\epsilon})$ and $z^*(\boldsymbol{\epsilon})$ denote an optimal solution of linear program (2). Define $n^*_j(t, \boldsymbol{\epsilon})$ and $q^{*,a}_j(t, \boldsymbol{\epsilon})$ as

\begin{align*}n^*_j(t, \boldsymbol{\epsilon}) = \sum_{a=0}^A x^{*,a}_j(t, \boldsymbol{\epsilon}), q^{*,a}_j(t, \boldsymbol{\epsilon})=\begin{cases}\frac{x^{*,a}_j(t, \boldsymbol{\epsilon})}{ n^*_j(t, \boldsymbol{\epsilon}) } &\text{if } n^*_j(t, \boldsymbol{\epsilon}) \gt 0,\\1\{a=0\} & \text{if } n^*_j(t, \boldsymbol{\epsilon}) = 0.\end{cases}\end{align*}

Note that, by definition, we have $\smash{n^*_j(1, \boldsymbol{\epsilon})} =n_{j}$ for all j. Note that $\smash{n^*_j(t,\boldsymbol{\epsilon})}$ is the total number of bandits in state j at time t (in the deterministic system) and $\smash{q^{*,a}_j(t,\boldsymbol{\epsilon})}$ is the fraction of bandits in state j at time t that receive action a. Let $Z_{jl}(t,\boldsymbol{\epsilon})$ be a categorical random variable on the set $\mathbb{A}$ with the property that

\begin{equation}\mathbb{P}(Z_{jl}(t, \boldsymbol{\epsilon}) = a) = q^{*,a}_j(t, \boldsymbol{\epsilon})\quad\text{for all } a \in \mathbb{A},\end{equation}

where $\ell$ denotes a bandit. Also, let $\boldsymbol{\varphi} = (\varphi_1,\varphi_2, \dots, \varphi_{n_{{\rm tot}}})$ denote a permutation (or ordering) of all of the bandits in the system, where $n_{{\rm tot}} =\smash{\sum_{j=1}^J n_{j}},$ and let $j(t, \ell)$ denote the state of bandit $\ell$ at time t before any action is applied. The description of RAC is given below.

The idea behind RAC is as follows. At each time t, we check the bandits in the order prescribed by $\varphi$. Suppose that we are currently checking bandit $\ell$. If its current state is j then we apply action a to $\ell$ with probability $\smash{q^{*,a}_j(t,\mathbf{\epsilon})}$ as long as we have not exceeded the allowed budget in period t; otherwise, we simply do not apply a proper treatment to $\ell$. It turns out that the order in which we assess the bandits in RAC (i.e. the permutation $\boldsymbol{\varphi}$) does not affect our analysis and results. Thus, for simplicity, we will simply use $\varphi(k) = k$ for all k. The choice of $\boldsymbol{\epsilon}$, on the other hand, significantly affects the performance of RAC. To see this, note that we can decompose the loss of RAC as

\begin{equation}0 \leq V^{{\rm RAC}} - V^{D}(\mathbf{0}) = [V^{{\rm RAC}} -V^D(\boldsymbol{\epsilon})] + [V^D(\boldsymbol{\epsilon}) - V^D(\mathbf{0})].\end{equation}

First note that the first inequality follows from Lemma 1 since $\smash{V^{{\rm RAC}} \geq V^S}$. Next, observe that if all components of $\boldsymbol{\epsilon}$ are close to 0, then $V^D(\boldsymbol{\epsilon})- V^D(\mathbf{0})$ is also close to 0. However, $V^{{\rm RAC}} -V^D(\boldsymbol{\epsilon})$ could be large. This is so because the deviations of $X^a_j(t)$ from $x^{*,a}_j(t,\boldsymbol{\epsilon})$ (due to the stochastic nature of the system) are likely to lead RAC to exhaust the activation budget at time t before those bandits that should have been activated are even considered for activation. As a result, many bandits may end up in undesirable states at the beginning of period $T+1$. On the other hand, if the components of $\boldsymbol{\epsilon}$ are large, $V^{{\rm RAC}} -V^D(\boldsymbol{\epsilon})$ will be close to 0. This is because there is a negligible probability that we will exhaust the activation budget in each period and we will be able to activate all bandits that need to be activated. However, the term $V^D(\boldsymbol{\epsilon}) - V^D(\mathbf{0})$ can be large (see Lemma 2 below). It follows that care must be taken when choosing $\boldsymbol{\epsilon}$.

The following lemma provides an upper bound for $V^D(\boldsymbol{\epsilon}) -V^D(\mathbf{0})$.

Lemma 2. There exists a constant $M \gt 0$, independent of T, and a vector $\boldsymbol{\epsilon} \ge \boldsymbol{0}$ satisfying $\epsilon_t \le b_t$ for all t such that

\begin{equation}V^D(\boldsymbol{\epsilon}) - V^D(\mathbf{0}) \le M \cdot \bigg[\sum_{t=1}^T\epsilon_t\bigg].\end{equation}

Lemma 2 says that $V^D(\boldsymbol{\epsilon}) - V^D(\mathbf{0})$ is roughly proportional to $\smash{\sum_{t=1}^T \epsilon_t}$. The proof of this lemma can be found in Appendix A and utilizes a standard duality argument for LP sensitivity analysis.

Let $C^{{\rm RAC}}$ denote a realization of total costs under RAC. By definition, we have $\mathbb{E}[C^{{\rm RAC}}] = V^{{\rm RAC}}$. The following lemma tells us that $C^{{\rm RAC}} - V^D(\epsilon)$ too is roughly proportional to $\smash{\sum_{t=1}^T \epsilon_t}$, with a positive probability.

Lemma 3. For any $\boldsymbol{\epsilon} \ge \boldsymbol{0}$, we have

(3)\begin{equation}C^{{\rm RAC}} - V^D(\boldsymbol{\epsilon}) \le [\phi + \max_{a,j} c^a_j ] \cdot \bigg[\sum_{t=1}^T \epsilon_t\bigg]\label{eqn11}\end{equation}

with probability at least

(4)\begin{align}1 - 3 &\cdot A \cdot J^2 \cdot \sum_{t=1}^T \exp\bigg\{-\frac{\epsilon^2_t}{12 \cdot A^2 \cdot J^4 \cdot [\sum_{\ell = 1}^Jn_{\ell}] }\bigg\} \label{eqn12}\end{align}

\begin{align}&\quad- 3 \cdot A \cdot J^2 \cdot \sum_{t=1}^T \exp\bigg\{- \frac{\epsilon_t}{6\cdot A \cdot J^2 }\bigg\}. \nonumber\end{align}

If the components of $\boldsymbol{\epsilon}$ are uniformly small then the bound in (3) is small and holds with a small probability. If, on the other hand, the components of $\boldsymbol{\epsilon}$ are uniformly large, the bound in (3) is large and holds with a large probability. Ideally, we would like to choose $\boldsymbol{\epsilon}$ that yields a small bound in (3) that holds with a large probability.

To better characterize the impact of the choice of $\boldsymbol{\epsilon}$ on the performance of RAC, we now turn to an asymptotic setting where $n_{j}$ (for all j), $b_t$ (for all t), and m are uniformly scaled by a factor of $\theta \gt 0$. This is the standard asymptotic setting considered in the RMABP literature (e.g. Weber and Weiss (Reference Weber and Weiss1990) and Verloop (Reference Verloop2016)). Let $\smash{V^S_{\theta}}$ and $\smash{V^{\pi}_{\theta}}$ denote the expected total costs under the optimal control and a feasible control $\pi$, respectively, when we use the scaling factor $\theta$. Also, let $\smash{V^D_{\theta}(\boldsymbol{\epsilon})}$ denote the optimal value of the corresponding linear program (2). It is not difficult to see that the optimal solution of linear program (2) for $\boldsymbol{\epsilon} =\mathbf{0}$ and $\theta \gt 0$ is given by

\begin{equation}x^{*,a}_{\theta, j}(t, \mathbf{0}) = \theta \cdot x^{*,a}_j(t,\mathbf{0}) \quad\text{and}\quad z^{*}_{\theta}(\mathbf{0}) = \theta\cdot z^{*}(\mathbf{0}),\end{equation}

which implies that $\smash{V^D_{\theta}(\mathbf{0}) = \theta \cdot V^D(\mathbf{0})}$. We also define the terms $\smash{n^*_{\theta,j}(t,\boldsymbol{\epsilon})}$ and $\smash{q^{*,a}_{\theta,j}}(t, \boldsymbol{\epsilon})$ analogously to how $\smash{n^*_j(t, \boldsymbol{\epsilon})}$ and $\smash{q^{*,a}_j(t, \boldsymbol{\epsilon})}$ are defined, with $\smash{x^{*,a}_j(t, \mathbf{0})}$ being replaced with $\smash{x^{*,a}_{\theta, j}(t, \mathbf{0})}$. Let $n_{{\rm tot}} \,:\!=\smash{\sum_{i=1}^J n_{i}},$ and define $S(t) \,:\!= \smash{\{(a,j)\colon x^{*,a}_j(t, \mathbf{0}) \gt 0 \mbox{ and } c^a_j> 0\}}$. We state our first theorem below.

Theorem 1. Let d be a positive number. Suppose that $S(t) \not= \emptyset$ for all t, and let

\begin{equation} \epsilon_t = 6 \cdot A \cdot J^2 \cdot \sqrt{d \cdot n_{{\rm tot}} \cdot \theta \cdot \ln \theta }\end{equation}

for all time periods t. For all $\theta \gt \theta^*$ where

(5)\begin{equation}\frac{\theta^*}{\ln \theta^*} = \frac{36 \cdot A^2 \cdot J^4 \cdot d\cdot n_{{\rm tot}}}{\min_t b_t},\label{eqn13}\end{equation}

there exists a constant $M \gt 0$ independent of T, d, and $\theta$ such that

\begin{equation} \frac{V^{{\rm RAC}}_{\theta} - V^S_{\theta}}{V^S_{\theta}} \le M \cdot \bigg[\frac{T}{\theta^{d}} + \sqrt{\frac{d \cdot \ln\theta}{\theta}}\bigg].\end{equation}

Proof. Fix $d \gt 0$. It is not difficult to see that condition (5) is equivalent to $\theta \cdot b_t \gt \epsilon_t$ for all t. Its role is to guarantee that the right-hand side of budget constraints in linear program (2) is positive and, therefore, the linear program has a feasible solution. Now, note that $\smash{\sqrt{d \cdot n_{{\rm tot}}\cdot \theta \cdot \ln \theta}} \ge d \cdot \ln \theta$ for all large $\theta$. Thus, by Lemma 3, under the choice of $\boldsymbol{\epsilon}$ in Theorem 1,

\begin{equation}C^{{\rm RAC}}_{\theta} - V^D_{\theta}(\boldsymbol{\epsilon}) \, = \,O(T \cdot \sqrt{d \cdot \theta \cdot \ln \theta})\end{equation}

with probability at least

\begin{equation}1 - \Theta\bigg(\frac{T}{\theta^{3d}}\bigg) -\Theta\bigg(\frac{T}{\theta^d}\bigg) = 1 -\Theta\bigg(\frac{T}{\theta^d}\bigg).\end{equation}

Since there is a total of $\theta \cdot n_{{\rm tot}}$ bandits in each period, $\smash{C^{{\rm RAC}}_{\theta} - V^D_{\theta}(\boldsymbol{\epsilon})}$ is at most $\Theta(T \cdot \theta \cdot n_{{\rm tot}})$ with probability at most $\Theta({T}/{\theta^d})$. We conclude that

\begin{equation}V^{{\rm RAC}}_{\theta} - V^D_{\theta}(\boldsymbol{\epsilon}) = O\bigg(T \cdot \sqrt{d \cdot \theta \cdot \ln \theta} + \frac{T^2}{\theta^{d-1}}\bigg).\end{equation}

By Lemma 2, under the choice of $\boldsymbol{\epsilon}$ in Theorem 1, we have

\begin{equation}V^D_{\theta}(\boldsymbol{\epsilon}) - V^D_{\theta}(\mathbf{0}) = O(T \cdot \sqrt{d \cdot \theta \cdot \ln \theta}).\end{equation}

Putting the bounds for $V^{{\rm RAC}}_{\theta} -V^D_{\theta}(\boldsymbol{\epsilon})$ and $V^D_{\theta}(\boldsymbol{\epsilon}) -V^D_{\theta}(\mathbf{0})$ together yields

\begin{equation}V^{{\rm RAC}}_{\theta} - V^D_{\theta}(\mathbf{0}) = O\bigg(T \cdot \sqrt{d \cdot \theta \cdot \ln \theta} + \frac{T^2}{\theta^{d-1}}\bigg).\end{equation}

The proof is complete by noting that $\smash{{(V^{{\rm RAC}}_{\theta} -V^S_{\theta})}/{V^S_{\theta}}} \le \smash{{(V^{{\rm RAC}}_{\theta} -V^D_{\theta}(\mathbf{0}))}/{V^D_{\theta}(\mathbf{0})}}$ and $\smash{V^D_{\theta}(\mathbf{0})}$ is $\Theta(T \cdot \theta)$ (because $S(t) \not= \emptyset$ for all t, which implies that in each period we always incur a total cost that is at least proportional to $\theta$).□

The bound in Theorem 1 is for the expected total costs criteria, but our analysis can also be applied to the expected total discounted cost criteria. This is detailed in Proposition 1 below, which is the direct analogue of Theorem 1.

Proposition 1. Let d be a positive number. Suppose that we multiply the cost at time t with $\delta^{t-1}$ for some discount factor $\delta \in (0,1)$. Suppose that $S(t) \not= 0$ for all t. Let

\begin{equation}\epsilon_t = 6 \cdot A \cdot J^2 \cdot \sqrt{d \cdot n_{{\rm tot}} \cdot \ln(t+e-1) \cdot \theta \cdot \ln \theta \,}\end{equation}

for all t (note that $\epsilon_t$ now depends on t). For all $\theta> \theta^*$ where

\begin{equation}\frac{\theta^*}{\ln \theta^*} = \frac{36 \cdot A^2 \cdot J^4 \cdot d\cdot n_{{\rm tot}} \cdot \log (T + {\rm e} - 1)}{\min_t b_t},\end{equation}

there exists a constant $M \gt 0$ independent of T and $\theta$ such that

\begin{equation}\frac{V^{{\rm RAC}}_{\theta} - V^S_{\theta}}{V^S_{\theta}} \le\frac{V^{{\rm RAC}}_{\theta} -V^D_{\theta}(\mathbf{0})}{V^D_{\theta}(\mathbf{0})} \le M \cdot \bigg[\frac{1}{\theta^d} + \sqrt{\frac{d \cdot \ln\theta}{\theta}}\bigg].\end{equation}

Proof. The proof is akin to the proof of Theorem 1. Here, we will only provide its outline. Note that, when we multiply the cost at time t with $\delta^{t-1}$ for some discount factor $\delta \in (0,1)$, the bound in Lemma 2 becomes

\begin{equation}V^D(\boldsymbol{\epsilon}) - V^D(\mathbf{0}) \le M \cdot \bigg[\sum_{t=1}^T\delta^{t-1} \cdot \epsilon_t\bigg]\end{equation}

for some $M \gt 0$ independent of T and $\boldsymbol{\epsilon} \ge \mathbf{0}$. The bound can be shown using exactly the same arguments as used in the proof of Lemma 2. See Remark 4 in Appendix A.1. Similarly, the bound in (3) becomes

\begin{equation}C^{{\rm RAC}} - V^D(\boldsymbol{\epsilon}) \le \Big[\phi + \max_{a,j} \, c^a_j\Big] \cdot \bigg[\sum_{t=1}^T \delta^{t-1} \cdot \epsilon_t\bigg],\end{equation}

which can be shown using exactly the same arguments as used in the proof of Lemma 3. See Appendix A.2. Next, observe that, akin to the proof of Theorem 1, our choice of $\epsilon_t$ implies that

\begin{equation}C^{{\rm RAC}}_{\theta} -V^D_{\theta}(\boldsymbol{\epsilon}) =O\bigg(\sum_{t=1}^T\delta^{t-1} \cdot \sqrt{d \cdot \ln(t+{\rm e}-1) \cdot \theta \cdot \ln \theta}\bigg) = O(\sqrt{d \cdot \theta \cdot \ln\theta})\end{equation}

with probability at least

\begin{equation}1 - \Theta\bigg(\sum_{t=1}^T \frac{1}{\theta^{3d \cdot \ln(t + {\rm e} -1)}}\bigg) - \Theta\bigg(\sum_{t=1}^T \frac{1}{\theta^{d \cdot \ln(t+{\rm e}-1)}}\bigg) = 1 - \Theta\bigg(\sum_{t=1}^T \frac{1}{(t+{\rm e}-1)^{d\cdot \ln \theta}}\bigg).\end{equation}

Additionally, $\smash{C^{{\rm RAC}}_{\theta} -V^D_{\theta}(\boldsymbol{\epsilon})}$ is at most

\begin{equation}\Theta \bigg(\sum_{t=1}^T \delta^{t-1} \cdot \theta \cdot n_{{\rm tot}}\bigg) = \Theta(\theta \cdot n_{{\rm tot}})\end{equation}

with probability at most $\smash{\Theta(\sum_{t=1}^T {1}/{(t+{\rm e}-1)^{d\cdot \ln \theta}})}.$ We conclude that

\begin{equation}V^{{\rm RAC}}_{\theta} - V^D_{\theta}(\mathbf{0}) = O\bigg(\sqrt{d \cdot \theta \cdot \ln \theta} + \sum_{t=1}^T \frac{\theta}{(t+{\rm e}-1)^{d \cdot \ln \theta}}\bigg).\end{equation}

Since $\smash{V^D_{\theta}(\mathbf{0})}$ is $\smash{\Theta(\sum_{t=1}^T\delta^{t-1} \cdot \theta)} = \Theta(\theta)$ (by our assumption that $S(t) \not= \emptyset$ for all t) and

\begin{equation}\sum_{t=1}^T \frac{1}{(t+{\rm e}-1)^{d \cdot \ln \theta}} \le \frac{1}{{\rm e}^{d\cdot \ln \theta}} + \int_1^{\infty} \frac{{\rm d}x}{(x+{\rm e}-1)^{d \cdot \ln\theta}} \le \frac{2}{{\rm e}^{d \cdot \ln \theta}} = \frac{2}{\theta^{d}}\end{equation}

for large $\theta$, we have the desired bound.□

5. Numerical experiments

5.1 RAC performance

In this section, we test the performance of RAC using two experiments. In the first experiment, we consider an instance of the RMABP with two states and two actions; in the second experiment we consider an instance of the RMABP with five states and five actions. In each experiment, we use $\epsilon_t = \smash{\sqrt{\theta \cdot \log \theta}}$ for all t and a wide range of $\theta$. For both experiments, we randomly generated the cost coefficients according to a uniformly distributed random variable U(0,10). We also generated the transition probabilities by first generating a uniformly distributed random variable U(0,1), and then normalizing the generated rows to obtain transition probability matrix P. We used the same budget $b_t=1$ for all t for both experiments. The remaining details regarding all other parameters can be found in the online supplement Zayas-Cabán et al. (2019). In what follows, we report the percentage losses of RAC (i.e. $\smash{{(V^{{\rm{RAC}}}_{\theta} - V^D_{\theta}(\mathbf{0}))}/{V^D_{\theta}(\mathbf{0})}}$) and their corresponding confidence intervals out of 1000 Monte Carlo runs for the two experiments in Tables 1 and 2, respectively. Additional simulation scenarios and their corresponding results can be found in the online supplement Zayas-Cabán et al. (2019).

Table 1: Percentage loss for two states and two actions.

Table 2: Percentage loss for five states and five actions.

As predicted by Theorems 1 and 2, RAC performs better as $\theta$ grows large. This confirms the asymptotic optimality of RAC. The performance of RAC for small $\theta$ depends on the problem parameters (see the online supplement Zayas-Cabán et al. (2019)). However, in most cases, the percentage loss of RAC is only a single digit number starting from $\theta = 10$, which suggests the robustness of RAC even for problems with relatively small $\theta$. Finally, we remark that while the expected percentage loss should be positive, it is possible to have negative percentage loss in simulations. This is due to sampling error in estimating the average costs using Monte Carlo simulations, and most of the negative percentage losses are very close to 0.

5.2 Comparing RAC with an index policy

In this section, we compare the performance of the RAC with a strict priority rule (or index policy). In the two examples below, we consider two states and two actions (i.e. treatment/active or no treatment/inactive) and assume that there is a budget of 1 (i.e. a maximum of one active bandit) in each time period. Our priority rule works as follows: we use our heuristic to calculate the probability of activating an arm in each state and then activate the arm with the highest probability directly. Note that, since we are activating an arm in each period, we will use the entire budget in each period. We then compare the costs for this priority rule with our heuristic. The results are summarized in Tables 3 and 4. Note that RAC outperforms the specific priority rule. Moreover, it does so without always using the full budget in every period. (We report in Table 5 the percentage budget used in RAC. That is, if the budget in each period is b and the number of periods under consideration is T, then the total budget allowed is $b\times T$. If the total budget actually used during the T time periods is B, then the percentage budget is calculated as ${B}/{(b\times T)}$.) This is perhaps not too surprising.

Table 3: Percentage loss for two states and one action example using RAC.

Table 4: Percentage loss for two states and one action using priority rule.

Table 5: Percentage budget used for two states and one action.

Since a treatment could be costly, it may not be optimal to always use up all the budgets; it is also important to take into account the trade-off between the cost of a treatment and its benefit.

6. Dynamic population model

In this section, we allow new bandits to arrive in each time period, and we also consider random departures (e.g. random service completions or random abandonments). The same notation as before will be used, with an addition that we also use $\Lambda_{jt}$ to denote the (random) arrival of bandits in state j at time t, where $\Lambda_{jt}$ is a Poisson random variable with mean $\lambda_{jt} \gt 0$. The stochastic control formulation of our problem can now be written as

\begin{gather} V^S = \min_{\pi \in \Pi} \mathbb{E}\bigg[\sum_{t=1}^T \sum_{a=0}^A\sum_{j=1}^J c^a_j \cdot X^{\pi, a}_j(t) + \phi \cdot \bigg(\sum_{a=0}^A\sum_{j \in \mathbb{U}} X^{\pi, a}_j(T+1) - m \bigg)^+\bigg]\\[2pt]\text{such that} \sum_{a=0}^A X^{\pi, a}_j(t) = \sum_{a=0}^A\sum_{i=1}^J Y^{\pi, a}_{ij}(t-1) + \Lambda_{jt} \quad\text{for all }j\ge 1, \, t \ge 2,\\[2pt]\sum_{a=0}^A X^{\pi, a}_j(1) = n_{j} + \Lambda_{j1} \quad\text{for all }j \ge 1,\\[2pt]\sum_{a=1}^A \sum_{j=1}^J X^{\pi, a}_j(t) \le b_t \quad\text{for all } t\ge 1.\end{gather}

The corresponding deterministic formulation is given by

(6a)\begin{equation}V^D(\mathbf{\epsilon}) = \min_{x, z} \sum_{t = 1}^T \sum_{a=0}^A \sum_{j= 1}^{J} c_{j}^a \cdot x_{j}^a(t, \boldsymbol{\epsilon}) + \phi \cdot z(\boldsymbol{\epsilon})\end{equation}

such that

(6b)\begin{equation}\sum_{a=0}^{A} x_{j}^a(t, \boldsymbol{\epsilon}) =\sum_{a=0}^A \sum_{i=1}^{J} x_{i}^a(t-1, \boldsymbol{\epsilon}) \cdot p^a_{ij} \,+ \, \lambda_{jt} \quad\text{for all } j\ge 1 , \, t \ge 2,\end{equation}

(6c)\begin{equation}\sum_{a = 0}^A x_{j}^{a} (1, \boldsymbol{\epsilon}) = n_{j} + \lambda_{j1}\quad\text{for all } j \ge 1,\end{equation}

(6d)\begin{equation}\sum_{a=1}^A \sum_{j=1}^J x^a_j(t, \boldsymbol{\epsilon}) \le b_t - \epsilon_t\quad\text{for all } t \ge 1,\end{equation}

(6e)\begin{equation}z(\boldsymbol{\epsilon}) \, \ge \, \sum_{j \in \mathbb{U}} \sum_{a=0}^A \sum_{i=1}^Jx^a_i(T, \boldsymbol{\epsilon}) \cdot p^a_{ij} - m,\end{equation}

(6f)\begin{equation}z(\boldsymbol{\epsilon}), \, x^a_j(t, \boldsymbol{\epsilon}) \ge 0 \quad\text{for all }\, a \ge 0, \, j \ge 1, \, t \ge 1.\end{equation}

Both Lemmas 1 and 2 still hold in the new setting (we omit the details), and the definition of RAC is still the same as that presented in Section 4, with one difference: since the number of bandits in the system potentially changes at every time period, the permutation $\boldsymbol{\varphi}$, which provides the order in which we assess the state of existing bandits in step 3 in the original definition of RAC, must now be selected at the beginning of every time period. That said, as noted in Section 4, the exact choice of $\boldsymbol{\varphi}$ does not affect our asymptotic analysis and thus we can choose an arbitrary permutation $\boldsymbol{\varphi}$ that randomizes the order of existing bandits at each time period. The following result is the analogue of Lemma 3 and the proof can be found in Appendix B.

Lemma 4. For any $\boldsymbol{\epsilon} \ge \boldsymbol{0}$, we have

(7)\begin{equation}C^{{\rm RAC}} - V^D(\mathbf{\epsilon}) \le \Big[\phi + \max_{a,j} c^a_j\Big] \cdot \bigg[\sum_{t=1}^T \epsilon_t\bigg]\label{eqn20}\end{equation}

with probability at least

(8)\begin{align}&1 - 3 \cdot A \cdot J^2 \cdot \sum_{t=1}^T \exp\bigg\{\!-\frac{\epsilon^2_t}{48 \cdot A^2 \cdot J^4 \cdot [\sum_{\ell = 1}^Jn_{\ell}] }\bigg\}- 3 \cdot A \cdot J^2 \cdot \sum_{t=1}^T \exp\bigg\{\!-\frac{\epsilon_t}{12 \cdot A \cdot J^2 }\bigg\} \nonumber\\&- 3 \cdot A \cdot J^2 \cdot \sum_{t=1}^T\exp\bigg\{-\frac{\epsilon_t^2}{64 \cdot A^2 \cdot J^4 \cdot[\sum_{s=1}^t \sum_{i=1}^J \lambda_{is}]} \bigg\}.\label{eqn21}\end{align}

The last summation in (8) is due to the arrivals of new bandits at each time period. If $\lambda_{tj} = 0$ for all t and j then the last summation equals 0 and the bound in (8) is almost identical to the bound in (4), with the exception that the numbers 12 and 6 that appear in the denominator in the original bound have now been replaced with 48 and 12 in the denominator above. Let $\lambda_{{\rm tot}} \,:\!= \max_t \smash{\sum_{j=1}^J\lambda_{jt}}$. We consider the same asymptotic setting as in Section 4, where we also scale the arrival rate $\lambda_{jt}$ by $\theta$. Let $n_{{\rm tot}}$ and S(t) be as defined in Section 3. The following theorem is the analogue of Theorem 1.

Theorem 2. Let d be a positive number. Suppose that $S(t) \not= \emptyset$ for all t, and let

\begin{equation} \epsilon_t = 12 \cdot A \cdot J^2 \cdot \sqrt{t \cdot d \cdot \max\{n_{{\rm tot}}, \, \lambda_{{\rm tot}}\} \cdot \theta \cdot \ln\theta \,}\end{equation}

for all t. For all $\theta \gt \theta^*$ where

(9)\begin{equation}\frac{\theta^*}{\ln \theta^*} = 144 \cdot A^2 \cdot J^4 \cdot d \cdot \max\{n_{{\rm tot}}, \lambda_{{\rm tot}}\} \cdot \Big[\max_t\frac{t}{b_t} \Big],\label{eqn22}\end{equation}

there exists a constant $M \gt 0$ independent of T and M such that

(10)\begin{equation}\frac{V^{{\rm RAC}}_{\theta} - V^S_{\theta}}{V^S_{\theta}} \le M \cdot \bigg[\frac{T^{3/2}}{\theta^{d/2}} + \sqrt{\frac{d \cdot T \cdot \ln\theta}{\theta}}\bigg].\label{eqn23}\end{equation}

Proof. The proof is similar to that of Theorem 1. First, condition (9) is equivalent to $\theta \cdot b_t \gt \epsilon_t$ for all t. It is important to guarantee that the right-hand side of the budget constraints in linear program (6) is positive and, therefore, the linear program has a feasible solution.

Let E denote the event where (7) is satisfied. By Lemma 4 and our choice of $\boldsymbol{\epsilon}$, $\mathbb{P}(E)$ is at least $1 - \Theta({T}/{\theta^d})$. Now, we consider the event $E^c$. Unlike in the proof of Theorem 1 where we can simply bound $\smash{C^{{\rm RAC}}_{\theta} - V^D_{\theta}(\boldsymbol{\epsilon})}$ with a number that is of order $\Theta(T \cdot \theta \cdot n_{{\rm tot}})$, we now need to take into account new bandits that arrive at each time period. At time t, we have at most $\theta\cdot n_{{\rm tot}} +\smash{\sum_{s=1}^t \sum_{j=1}^J \Lambda^{\theta}_{js}}$ bandits present in the system, where $\Lambda^{\theta}_{js}$ is a Poisson random variable with rate $\theta \cdot \lambda_{js}$. As a result, we can bound:

(11)\begin{align}&\mathbb{E}[(C^{{\rm RAC}}_{\theta} - V^D_{\theta}(\boldsymbol{\epsilon})) \cdot \mathbf{1}\{E^c\}]\nonumber\\& \le \Big[\phi + \max_{a,j} c^a_j\Big] \cdot \mathbb{E}\bigg[\bigg( T\cdot \theta \cdot n_{{\rm tot}} + \sum_{t=1}^T t \cdot \bigg(\sum_{j=1}^J \Lambda^{\theta}_{j,t}\bigg)\bigg) \cdot \mathbf{1}\{E^c\}\bigg].\label{eqn24}\end{align}

Note that $\smash{\sum_{t=1}^T t \cdot (\sum_{j=1}^J\Lambda^{\theta}_{jt})}$ is stochastically dominated by the random variable $X \sim $ Pois$(T^2 \cdot \theta \cdot \lambda_{{\rm tot}})$. By the Cauchy–Schwarz inequality, we can further bound (11) as

(12)\begin{align}\mathbb{E}[(C^{{\rm RAC}}_{\theta} - V^D_{\theta}(\boldsymbol{\epsilon})) \cdot \mathbf{1}\{E^c\}] &\le \Big[\phi + \max_{a,j} c^a_j\Big] \cdot \mathbb{E}[(T \cdot \theta \cdot n_{{\rm tot}} + X)^2]^{1/2} \cdot \mathbb{P}(E^c)^{1/2} \nonumber\\& = O\bigg(\frac{T^{5/2}}{\theta^{d/2 - 1}}\bigg)\label{eqn25}\end{align}

where the last equality follows since

\begin{equation}\mathbb{P}(E^c)^{1/2} = O\bigg(\frac{T^{1/2}}{\theta^{d/2}}\bigg)\quad\mbox{and}\quad \mathbb{E}[X^2]^{1/2} = T^2 \cdot \theta \cdot \lambda_{{\rm tot}},\end{equation}

which implies that $\mathbb{E}[(T \cdot \theta \cdot n_{{\rm tot}} +X)^2]^{1/2} = O(T^2 \cdot \theta)$.

Putting (12) together with

\begin{equation}\mathbb{E}[(C^{{\rm RAC}}_{\theta} - V^D_{\theta}(\boldsymbol{\epsilon})) \cdot \mathbf{1}\{E\}] = O(T^{3/2} \cdot \sqrt{d \cdot \theta \cdot \log\theta})\end{equation}

(by Lemma 4 and our choice of $\boldsymbol{\epsilon}$) and the fact that $\smash{V^D_{\theta}(\mathbf{0})} = \Omega(T \cdot \theta)$ (because $S(t) \not= \emptyset$ for all t) yields the desired results. This completes the proof of Theorem 2.□

Proposition 2. Let d be a positive number. Suppose that we multiply the cost at time t with $\delta^{t-1}$ for some discount factor $\delta \in (0,1)$. Suppose that $S(t) \not= 0$ for all t. Let

\begin{equation}\epsilon_t = 12 \cdot A \cdot J^2 \cdot \sqrt{t \cdot \log(t + {\rm e} - 1)\cdot d \cdot \max\{n_{{\rm tot}}, \lambda_{{\rm tot}}\} \cdot \theta\cdot \ln \theta }\end{equation}

for all t (note that $\epsilon_t$ now depends on t). For all $\theta> \theta^*$ where

\begin{equation}\frac{\theta^*}{\ln \theta^*} = 144 \cdot A^2 \cdot J^4 \cdot d \cdot \max\{n_{{\rm tot}}, \lambda_{{\rm tot}}\} \cdot \log(T+{\rm e}-1) \cdot \Big[\max_t \frac{t}{b_t} \Big],\end{equation}

there exists a constant $M \gt 0$ independent of T and $\theta$ such that

\begin{equation}\frac{V^{{\rm RAC}}_{\theta} - V^S_{\theta}}{V^S_{\theta}} \le\frac{V^{{\rm RAC}}_{\theta} - V^D_{\theta}(\mathbf{0})}{V^D_{\theta}(\mathbf{0})} \le M \cdot \bigg[\frac{1}{\theta^{d/2}} +\sqrt{\frac{d \cdot \ln \theta}{\theta}}\bigg].\end{equation}

The proof of Proposition 2 is very similar to the proof of Theorem 2 and, therefore, is omitted. Unlike in Remark 2, where the bound in Proposition 1 holds for $T= o(\exp \{{\theta}/{(d \cdot \ln \theta)} \})$, the bound in Proposition 2 holds so long as

\begin{equation}T \cdot \log T = o\bigg(\frac{\theta}{\log \theta}\bigg).\end{equation}

In what follows, we discuss a special case of the dynamic population model in which bandits can arrive and leave randomly over time.

6.1. Bandit departure and abandonment

Suppose now that each bandit waits to receive a proper treatment ($a \ge1$) and immediately leaves the system once it receives a proper treatment (i.e. a service completion). Moreover, suppose that, if a bandit has not received a proper treatment up until the end of the current time period, it will stay in the system for the next period with probability $\alpha\in (0,1)$ and will leave (i.e. abandon) the system with probability $1 - \alpha$. These features can be included in the current modeling framework by introducing an extended state variable that includes the type of treatment a bandit received in the previous time period and a tracking variable that indicates whether the bandit is still in the system at the current time period. For example, we can let $j_t = (j'_t,a_{t-1}, \ell_t)$ be our extended state variable, where $j'_t$ is the current state/condition of the bandit, $a_{t-1}$ is the treatment received in the previous time period, and $\ell_t \in \{0, 1\}$ an indicator of the bandit’s presence in the system with 0 (1) denoting it left (remains in) the system at the beginning of period t. The one-step transition probabilities for this model now satisfy the following conditions:

\begin{align}\label{eqn26}\end{align}

(13)\begin{align}p^{a''}_{(j', a', 0), (j'', a'', 1)} = 0 \quad\text{for all } j',a', j\prime\prime,a\prime\prime, \label{eqn27}\end{align}

(14)\begin{align}p^{a''}_{(j',0,1), (j'', a'',1)} = 0 \quad\text{for all } j', j\prime\prime, a\prime\prime\ge 1, \label{eqn28}\end{align}

(15)\begin{align}\sum_{j''} p^0_{(j', 0, 1), \, (j'', 0, 1)} = \alpha \quad\text{for all }j'.\label{eqn29}\end{align}

Condition (13) simply says that a bandit that leaves the system never returns; condition (14) says that a bandit immediately leaves the system after receiving a proper treatment; and condition (15) says that the probability a bandit stays in the system for the next time period given that it has not received a proper treatment is $\alpha$.

For simplicity, suppose that, at the end of the horizon, we only get penalized for the number of bandits who leave the system without being properly treated (i.e. abandonments). This requires keeping track of the number of abandonments at each time period. To do this, we set $\smash{c^a_{(j',0,1)}} = M$ for all j’ and $a \ge 1$, where M is a sufficiently large number, and ${c^0_{(j',0,1)}} = 0$ for all j’. Since applying a proper action to a nonexisting bandit is very costly, this setup implies that, once a bandit leaves the system without being properly treated, the last two components of its state remains (0,1) throughout the horizon. As a result, the total number of abandonments throughout the horizon simply equals $\smash{\sum_{a=0}^A \sum_{j'}X^{a}_{(j', 0, 1)}}(T+1)$ and its corresponding total penalty costs is

\begin{equation}\phi \cdot \bigg(\sum_{a=0}^A \sum_{j'} X^{a}_{(j', 0, 1)}(T+1) -m\bigg)^+,\end{equation}

where $m \ge 0$ can be interpreted as the maximum number of acceptable abandonments.

The setting with service completions and abandonments discussed above is a special case of the more general setting considered in Theorem 2. As a result, the bound in (10) still holds in this setting. In fact, as previously noted, this bound can be significantly tightened by exploiting the special structure of the problem. That is, we have the following result.

Proposition 3. Let d be a positive number. Consider the setting with bandit departure and abandonment described above. Suppose that $S(t) \not= \emptyset$ for all t. If we let

\begin{equation}\epsilon_t = \frac{12 \cdot A \cdot J^2}{\min_t b_t} \cdot \sqrt{\bigg(\sum_{s=0}^{t-1} \alpha^s \bigg) \cdot d \cdot \max\{n_{{\rm{tot}}}, \lambda_{{\rm tot}}\} \cdot \theta \cdot \ln \theta }\end{equation}

for all t, then

\begin{equation}\frac{V^{{\rm RAC}}_{\theta} - V^S_{\theta}}{V^S_{\theta}} =O\bigg(\frac{1}{1 - \alpha} \cdot \frac{T^{1/2}}{\theta^{d/2}} +\sqrt{\frac{d \cdot \ln \theta}{(1 - \alpha) \cdot \theta}}\bigg).\end{equation}

Proof. The proof of Proposition 3 is very similar to the proof of Theorem 2. Here, we will only provide its outline. Note that we can strengthen the bound in Lemma 5 by replacing the term $\smash{\sum_{s=1}^t \sum_{j=1}^J \lambda_{jt}}$ in the third summation in (8) with $\smash{\sum_{s=1}^t \alpha^{t-s}\cdot [\sum_{j=1}^J \lambda_{jt}]}$ (see Remark 5 in Appendix B.1). Let E denote the event where (7) is satisfied. By our choice of $\boldsymbol{\epsilon}$, we have

\begin{equation}\mathbb{E}[(C^{{\rm RAC}}_{\theta} - V^D_{\theta}(\boldsymbol{\epsilon})) \cdot \mathbf{1}\{E\}] = O\bigg(T \cdot \sqrt{\frac{d \cdot \theta \cdot \ln\theta}{1-\alpha}} \bigg).\end{equation}

As for $\smash{\mathbb{E}[(C^{{\rm RAC}}_{\theta} -V^D_{\theta}(\boldsymbol{\epsilon})) \cdot \mathbf{1}\{E^c\}]}$, since each bandit stays in the system with probability at most $\alpha$, we can upper bound (in the stochastic dominance sense) the number of initial bandits that are still present in the system at time t by Binomial $(\theta \cdot n_{{\rm tot}}, \alpha^{t-1})$. Similarly, we can also upper bound the number of bandits that arrived in the system at time s in state j and are still present in the system at time t by Poisson $(\theta \cdot \lambda_{js} \cdot \alpha^{t-s})$. Since $\mathbb{P}(E^c)= O({T}/{\theta^d})$, by similar arguments as in the proof of Theorem 2 (using the Cauchy–Schwarz inequality), we have

\begin{equation}\mathbb{E}[(C^{{\rm RAC}}_{\theta} - V^D_{\theta}(\boldsymbol{\epsilon})) \cdot \mathbf{1}\{E^c\}] = O\bigg(\frac{T \cdot \theta}{1 - \alpha}\bigg) \cdot O\bigg(\frac{T^{1/2}}{\theta^{d/2}}\bigg).\end{equation}

The proof is complete by noting that $\smash{V^D_{\theta}(\mathbf{0})} =\Omega(T \cdot \theta)$ (because $S(t) \not= \emptyset$ for all t).

Note that, if $T \sim \theta^n$, where n is a natural number, then we can choose $d \gt n$ and RAC is asymptotically optimal. This illustrates that RAC performs well for a large decision horizon T with Poisson arrivals and geometrically distributed waiting times.

7. Closing remarks

In this paper, we consider discrete-time, finite-horizon, RMABPs for fixed and dynamic populations of bandits. To our knowledge, we are the first to simultaneously allow for multiple actions and a nonstationary number of active bandits. We proposed a heuristic control, RAC, that is based on an LP relaxation of the original stochastic control formulation and showed that it is asymptotically optimal. Similar to Verloop (2016), our heuristic does not depend on any indexability properties. In contrast to Verloop (Reference Verloop2016), it does not require assumptions on the underlying structure of the Markov process generated by each policy or assumptions on the dynamics on the associated deterministic approximation.

As mentioned, we can include nonstationary transition probability matrices and nonstationary costs. We can also generalize to the setting when there are several constraints for several subsets of bandits or one constraint on the total number of activated bandits over the entire finite-time horizon. Since we consider an asymptotic regime where the initial budgets and the number of initial bandits are proportionally scaled, it can be shown that these constraints (one or multiple constraints) will not be violated with a ‘high’ probability. Also, having state-dependent terminal costs will not change the asymptotic order. Finally, the proposed approach would also work for the multi-action setting where each action contributes a different weight to the budget since once the linear program is solved, we can show that the budget constraints will not be violated with a very high probability. This holds regardless of whether each action contributes the same or different weights to the budget.

There are several avenues for further research. The approach in Hu and Frazier (Reference Hu and Frazier2017) does not directly translate to our setting, but our approach may be applicable to their model since we make no idexability assumptions. One topic for further exploration is to test our heuristic using their data in the specific examples considered in their study. Developing and analyzing alternative policies is of clear interest. Comparing alternative policies with RAC in a numerical study can be considered. One might consider the same model but allow for states to be observed only when a treatment is applied. This would extend the model considered by Deo et al. (Reference Deo2013) that is motivated by capacity management in healthcare. This model requires keeping track of additional information (i.e. the time between interventions) and, hence, will require a larger state space. In this case, RAC may no longer be asymptotically optimal and, more broadly, the model would present significant technical challenges for the analysis and computation of good control policies. Another consideration is the possibility of allowing more general (i.e. nonlinear) cost structure and constraints. For instance, to capture the fact that, for many settings (e.g. EMS response, humanitarian logistics), the number of available servers is random, the budget $b_t$ may be assumed to be a random variable for each t. Notwithstanding, RMABPs are a very useful modeling framework that can be applied to a broad range of problems; they provide significant technical challenges for optimal control, which are a bright research direction.