Drift Control of High-Dimensional Reflected Brownian Motion: A Computational Method Based on Neural Networks

1. Introduction

Beginning with the seminal work of Iglehart and Whitt (1970a, b), a large literature has been developed over the last 50+ years that justifies the use of reflected Brownian motions (RBMs) as approximate models of queueing systems under “heavy traffic” conditions. In particular, a limit theorem proved by Reiman (1984) justifies the use of d-dimensional RBM as an approximate model of a d-station queueing network. Reiman’s theory is restricted to networks of the generalized Jackson type, also called single-class networks, or networks with homogeneous customer populations, but it has been extended to more complex multiclass networks under certain restrictions, most notably by Peterson (1991) and Williams (1998b). The survey papers by Harrison and Nguyen (1993) and by Williams (1996) provide an overview of heavy traffic limit theory through its first 25 years.

Many authors have commented on the compactness and simplicity of RBM as a mathematical model, at least in comparison with the conventional discrete-flow models that it replaces. For example, in the preface to the book by Kushner (2001) on heavy traffic analysis, one finds the following passage.

These approximating [Brownian] models have the basic structure of the original problem, but are significantly simpler. Much inessential detail is eliminated … They greatly simplify analysis, design, and optimization, [yielding] good approximations to problems that would otherwise be intractable. (Kushner 2001)

Of course, having adopted RBM as a system model, one still confronts the question of how to do performance analysis, and in that regard, there has been an important recent advance; Blanchet et al. (2021) have developed a simulation-based method to estimate steady-state performance measures for RBM in dimensions up to 200, and those estimates come with performance guarantees.

1.1. Descriptive Performance Analysis vs. Optimal Control

Early work on heavy traffic approximations, including the papers cited above, focused on descriptive performance analysis under fixed operating policies. Harrison (1988, 2000) expanded the framework to include consideration of dynamic control using informal arguments to justify Brownian approximations for queueing network models where a system manager can make sequencing, routing, and/or input control decisions. Early papers in that vein by Harrison and Wein (1989, 1990) and by Wein (1991) dealt with Brownian models simple enough that their associated control problems could be solved analytically. But, for larger systems and/or more complex decisions, the Brownian control problem that approximates an original queueing control problem may only be solvable numerically. Such stochastic control problems may be of several different types depending on context.

At one end of the spectrum are drift control problems, in which the controlling agent can effect changes in system state only at bounded finite rates. At the other end of the spectrum are impulse control problems, in which the controlling agent can effect instantaneous jumps in system state, usually with an associated fixed cost. In between are singular control problems, in which the agent can effect instantaneous state changes of any desired size, usually at a cost proportional to the size of the displacement; see, for example, Karatzas (1983). In this paper, we develop a computational method for the first of those three problem classes, and then, we illustrate its use on selected test problems. Our method is a variant of the one developed by Han et al. (2018) for the solution of semilinear partial differential equations (PDEs), and in its implementation, we have reused substantial amounts of the code provided by Han et al. (2018) and Zhou et al. (2021a).

1.2. Literature Review

Two of the most relevant streams of literature are (i) drift rate control problems and (ii) solving PDEs using deep learning. Ata et al. (2005) considers a one-dimensional drift rate control problem on a bounded interval under a general cost of control but no state costs. The authors characterize the optimal policy in closed form, and they discuss the application of their model to a power control problem in wireless communication. Ormeci Matoglu and Vande Vate (2011) consider a drift rate control problem where a system controller incurs a fixed cost to change the drift rate. The authors prove that a deterministic, nonoverlapping control band policy is optimal; also see Vande Vate (2021). Ghosh and Weerasinghe (2007, 2010) extend Ata et al. (2005) by incorporating state costs and abandonments and by optimally choosing the interval where the process lives.

Drift control problems arise in a broad range of applications in practice. Rubino and Ata (2009) studies a dynamic scheduling problem for a make-to-order manufacturing system. The authors model order cancellations as abandonments from their queueing system. This model feature gives rise to a drift rate control problem in the heavy traffic limit. Ata et al. (2019) uses a drift control model to study a dynamic staffing problem in order to determine the number of volunteer gleaners (who sign up to help but may not show up) for harvesting leftover crops donated by farmers for the purpose of feeding food-insecure individuals. Bar-Ilan et al. (2007) use a drift control model to study international reserves. All of the papers mentioned above study one-dimensional drift rate control problems.

The recent working paper by Ata and Kasikaralar (2023) studies dynamic scheduling of a multiclass queue motivated by the call center industry. Focusing on the Halfin–Whitt asymptotic regime, the authors derive a (limiting) drift rate control problem whose state space is ${\mathbb{R}}^d$, where d is the number of buffers in their queueing model. Like us, those authors build on the earlier work by Han et al. (2018) to solve their (high-dimensional) drift rate control problem. However, our work differs from their work significantly because their control problem has no state-space constraints.

As mentioned earlier, our work builds on the seminal paper by Han et al. (2018). In the last five years, there have been many other papers written on solving PDEs using deep neural networks; see the recent surveys by E et al. (2022) and Beck et al. (2023).

1.3. The Remainder of This Paper

Section 2 recapitulates essential background knowledge from the RBM theory, after which Section 3 states in precise mathematical terms the discounted control and ergodic control problems that are the object of our study. In each case, the problem statement is expressed in probabilistic terms initially, and then, it is re-expressed analytically in the form of an equivalent Hamilton–Jacobi–Bellman (HJB) equation. Section 4 derives key identities that significantly contribute to the subsequent development of our computational method. Section 5 describes our computational method in detail.

Section 6 specifies three families of drift control test problems, each of which has members of dimensions $d=1,2,\dots$. The first two families arise as heavy traffic limits of certain queueing network control problems, and we explain that motivation in some detail. Drift control problems in the third family have a separable structure that allows them to be solved exactly by analytical means, which is of obvious value for assessing the accuracy of our computational method. Section 7 presents numerical results obtained with our method for all three families of test problems. In that admittedly limited context, our computed solutions are accurate to within a fraction of 1%, and our method remains computationally feasible up to at least dimension $d=30$ and in some cases, up to dimension 100 or more. In Section 8, we describe variations and generalizations of the problems formulated in Section 3 that are of interest for various purposes and that we expect to be addressed in future work. Finally, there are a number of appendices that contain proofs or other technical elaboration for arguments or procedures that have only been sketched in the body of the paper.

2. RBM Preliminaries

We consider here a reflected Brownian motion $Z=\{Z(t),t\ge 0\}$ with state space ${\mathbb{R}}_{+}^d,$ where $d\ge 1.$ The data of Z are a (negative) drift vector $\mu \in {\mathbb{R}}^d,$ a $d\times d$ positive-definite covariance matrix $A=(a_{ij}),$ and a $d\times d$ reflection matrix R of the form

The restriction to reflection matrices of the form (1) is not essential for our purposes, but it simplifies the technical development and is consistent with usage in the related earlier paper by Blanchet et al. (2021). Denoting by $W=\{W(t),t\ge 0\}$ a d-dimensional Brownian motion with zero drift, covariance matrix A, and $W(0)=0$, we then have the representation

Harrison and Reiman (1981) showed that the relationships (1)–(4) determine Y and Z as path-wise functionals of W and that the mapping $W\to (Y,Z)$ is continuous in the topology of uniform convergence. We interpret the ith column of R as the direction of the reflection on the boundary surface $S_i=\{z\in {\mathbb{R}}_{+}^d:z_i=0\},$ and we call $Y_i=\{Y_i(t),t\ge 0\}$ the “pushing process” on that boundary surface.

In preparation for future developments, let f be an arbitrary $C^2$ (that is, twice continuously differentiable) function ${\mathbb{R}}^d\to \mathbb{R}$, and let $\nabla f$ denote its gradient vector as usual. Also, we define a second-order differential operator $\mathcal{L}$ via

and a first-order differential operator $\mathcal{D}={({\mathcal{D}}_1,\dots ,{\mathcal{D}}_d)}^{\top }$ viawhere $\top$ in (6) denotes transpose. Thus, ${\mathcal{D}}_{i}f(·)$ is the directional derivative of f in the direction of the reflection on the boundary surface $S_i=\{z\in {\mathbb{R}}_{+}^d:z_i=0\}.$ With these definitions, an application of Ito’s formula now gives the following identify (cf. Harrison and Reiman 1981, section 3):

In the obvious way, the first inner product on the right side of (7) is shorthand for a sum of d Ito differentials, whereas the last one is shorthand for a sum of d Riemann–Stieltjes differentials.

3. Problem Statements and HJB Equations

Let us now consider a stochastic control problem whose state space is ${\mathbb{R}}_{+}^d$$(d\ge 1).$ The controlled process Z has the form

where (i) $W=\{W(t),t\ge 0\}$ is a d-dimensional Brownian motion with zero drift, covariance matrix A, and $W(0)=0$ as in Section 2, (ii) $\theta =\{\theta (t),t\ge 0\}$ is a nonanticipating control or nonanticipating drift process chosen by a system manager and taking values in a bounded set $\mathrm{\Theta }\subset {\mathbb{R}}^d,$ and (iii) $Y=\{Y(t),t\ge 0\}$ is a d-dimensional pushing process with components $Y_i$ that satisfy (3) and (4). Note that our sign convention on the drift in the basic system Equation (8) is not standard. That is, we denote by $\theta (t)$ the negative drift vector at time t.

The control $\theta$ is chosen to optimize an economic objective (see below), and attention will be restricted to stationary Markov controls or stationary control policies, by which we mean that

Hereafter, the set $\mathrm{\Theta }$ of drift vectors available to the system manager will be referred to as the action space for our control problem, and a function $u:{\mathbb{R}}_{+}^d\to \mathrm{\Theta }$ will simply be called a policy. We denote by $Z^u$ the controlled RBM defined via (8) and (9), and we denote by $Y^u$ the associated d-dimensional boundary pushing process.

Before specifying the system manager’s economic objective, we establish the following terminology; for $m,n\ge 1$, a function $g:D\subset {\mathbb{R}}^m\to {\mathbb{R}}^n$ is said to have polynomial growth if there exist constants ${\alpha }_1,\hspace{0.17em}{\beta }_1>0$ such that

Because the action space $\mathrm{\Theta }$ is bounded, for a function $g:{\mathbb{R}}_{+}^d\times \mathrm{\Theta }\to \mathbb{R}$, the polynomial growth assumption reduces to the following:

where ${\alpha }_2,\hspace{0.17em}{\beta }_2$ are positive constants.

With regard to the system manager’s objective, we take as given a continuous cost function $c:{\mathbb{R}}_{+}^d\times \mathrm{\Theta }\to \mathbb{R}$ with polynomial growth and a vector $\kappa \in {\mathbb{R}}_{+}^d$ of penalty rates associated with pushing at the boundary. (As it happens, the boundary penalty rates are all zeros for the numerical examples considered in this paper, but positive penalty rates will be needed in future applications.) The cumulative cost incurred over the time interval [0, t] under policy u is

Because our action space $\mathrm{\Theta }$ is bounded by assumption, the controlled RBM $Z^u$ has bounded drift under any policy u, from which one can prove the following mild but useful property; see Appendix A for its proof.

Proposition 1.

Under any policy u and for any integer $n=1,2,\dots$, the function

$$g_n(z,t)={\mathbb{E}}_z\left\{|Z^u(t){|}^n\right\},\hspace{0.17em}t\ge 0,$$

has polynomial growth in t for each fixed $z\in {\mathbb{R}}_{+}^d$.

3.1. Discounted Control

In our first problem formulation, an interest rate $r>0$ is taken as given, and we adopt the following discounted cost objective; choose a policy u to minimize

where ${\mathbb{E}}_z(·)$ denotes a conditional expectation given that $Z(0)=z.$ Given the polynomial growth Condition (10), it follows from Proposition 1 that the moments of $Z(t)$ have polynomial growth as functions of t for each fixed initial state z. Also, because $\mathrm{\Theta }$ is bounded by assumption, one can easily derive an affine bound for ${\mathbb{E}}_z\{\kappa ·Y^u(t)\}$ viewed as a function of t. Given the assumed positivity of the interest rate r, the expectation in (12) is, therefore, well defined and finite for each $z\in {\mathbb{R}}_{+}^d$.

Hereafter, we refer to $V^u(·)$ as the value function under policy u, and we define the optimal value function

where $\mathcal{U}$ is the set of stationary Markov control policies.

Let u be an arbitrary policy, and let $f:{\mathbb{R}}_{+}^d\to \mathbb{R}$ be a $C^2$ function with polynomial growth. Also, we define the second-order differential operation associated with policy u:

where $\mathcal{L}$ is defined via (5). The identity (7) continues to hold for the controlled process $Z^u$ if one replaces the constant (reference) drift vector $\mu$ with the state-dependent drift function $\theta (z)=u(z)$, and then, standard arguments as in Harrison (2013, section 6.3) lead to the following identity:

Comparing (14) with (12), one is led to the following PDE for $V^u$:

with boundary conditions

That is, if there is a $C^2$ solution of the PDE (15) and (16) that has polynomial growth, then it is equal to the value function under policy u.

The corresponding HJB equation, to be solved for the optimal value function $V(·),$ is

with boundary conditions

Moreover, the policy

is optimal, meaning that $V^{u^{\ast }}(z)=V(z)$ for $z\in {\mathbb{R}}_{+}^d$.

There will be no attempt here to prove the existence of $C^2$ solutions, but our computational method proceeds as if that was the case, striving to compute a $C^2$ function V that satisfies (17) and (18) as closely as possible in a certain sense. As an aside, whenever we refer to a $C^2$ function on a closed set, we mean that it is $C^2$ in an open neighborhood of the set.

In Appendix B.1, we use (7) to verify that a sufficiently regular solution of the PDE (15) and (16) does, in fact, satisfy (12) as intended and similarly, that a sufficiently regular solution of (17) and (18) does, in fact, satisfy (13).

3.2. Ergodic Control

For our second problem formulation, it is assumed that

Readers will see that our analysis can be extended to cost functions that take on negative values in at least some states, but to do so, one must deal with certain irritating technicalities. To be specific, the issue is whether the expected values involved in our formulation are well defined.

In preparation for future developments, let us recall that a square matrix R of the form (1), called a Minkowski matrix in linear algebra (or just M-matrix for brevity), is nonsingular, and its inverse is given by the Neumann expansion

Hereafter, we assume that

It is known that an RBM with a nonsingular covariance matrix, reflection matrix R, and negative drift vector $\theta$ has a stationary distribution if and only if the inequality in (21) holds (cf. Harrison and Williams 1987, section 6). Of course, our statement of this “stability condition” reflects the nonstandard sign convention used in this paper. That is, $\theta$ denotes the negative drift vector of the RBM under discussion.

For our ergodic control problem, a policy function $u:{\mathbb{R}}_{+}^d\to \mathrm{\Theta }$ is said to be admissible if, first, the corresponding controlled RBM $Z^u$ has a unique stationary distribution ${\pi }^u$ and if, moreover,

for any function $f:{\mathbb{R}}_{+}^d\to \mathbb{R}$ with polynomial growth. For an admissible policy u, as done in Harrison and Williams (1987) and Dai and Harrison (1991), one can define the boundary measures ${\nu }_1^u,\dots ,{\nu }_d^u$ as follows:where A is any Borel subset of the boundary surface $S_i=\{z\in {\mathbb{R}}_{+}^d:z_i=0\}$ and ${\mathbb{E}}_{{\pi }^u}$ denotes a conditional expectation given that $Z(0)\sim {\pi }^u$. Our assumption (21) ensures the existence of at least one admissible policy u, as follows. Let $\theta \in \mathrm{\Theta }$ be a negative drift vector satisfying (21), and consider the constant policy $u(·)\equiv \theta$. The corresponding controlled process $Z^u$ is then an RBM having a unique stationary distribution ${\pi }^u$, as noted above. It has been shown in Budhiraja and Lee (2007) that the moment-generating function of ${\pi }^u$ is finite in a neighborhood of the origin, from which it follows that ${\pi }^u$ has finite moments of all orders. Thus, ${\pi }^u$ satisfies (22) for any function f with polynomial growth, so u is admissible.

Because our cost function $c(z,\theta )$ has polynomial growth and our action space $\mathrm{\Theta }$ is bounded, the steady-state average cost

is well defined and finite under any admissible policy u. The objective in our ergodic control problem is to find an admissible policy u for which ${\xi }^u$ is minimal.

Let u be an arbitrary admissible policy, and consider the following PDE:

with boundary conditions

If $(\xi ,v^u)$ solve this PDE and $v^u$ is $C^2$ with polynomial growth, then one can show that $\xi ={\xi }^u$, the steady-state average cost under policy u, and $v^u$ is the relative value function corresponding to policy u.

The HJB equation for ergodic control is again of a standard form, involving a constant $\xi$ (interpreted as the minimum achievable steady-state average cost) and a relative value function $v:{\mathbb{R}}_{+}^d\to \mathbb{R}$. To be specific, the HJB equation is

with boundary conditions

Paralleling the previous development for discounted control, we show the following in Appendix B.2; if a $C^2$ function v and a constant $\xi$ jointly satisfy (26) and (27), then

where $\mathcal{U}$ denotes the set of admissible controls for the ergodic cost formulation. Moreover, the policyis optimal, meaning that ${\xi }^{u^{\ast }}=\xi .$ Again, paralleling the previous development for discounted control, we make no attempt to prove that such a solution for (26) and (27) exists. In Appendix B.2, we use (7) to verify that a sufficiently regular solution of the PDE (24) and (25) does, in fact, satisfy (23) as intended and similarly, that a sufficiently regular solution of (26) and (27) does, in fact, satisfy (28).

4. Equivalent SDEs

In this section, we prove two key identities, Equations (32) and (44) below, that are closely patterned after results used by Han et al. (2018) to justify their “deep [backward stochastic differential equation] method” for the solution of certain nonlinear PDEs. That earlier work provided both inspiration and detailed guidance for our study, but we include these derivations to make the current account as nearly self-contained as possible. Sections 4.1 and 4.2 treat the discounted and ergodic cases, respectively.

Our method begins by specifying what we call a reference policy. This is a nominal or default policy, specified at the outset but possibly revised in light of computational experience, that we use to generate sample paths of the controlled RBM Z. Roughly speaking, one wants to choose the reference policy so that paths of the state process tend to occupy parts of the state space thought to be most frequently visited by an optimal policy.

4.1. Discounted Control

Our reference policy for the discounted case chooses a constant action $u(z)=\tilde{\theta }>0$ in every state $z\in {\mathbb{R}}_{+}^d$. (Again, we stress that given the nonstandard sign convention embodied in (8) and (20), this means that $\tilde{Z}$ has a constant drift vector $-\tilde{\theta }$, with all components negative.) Thus, the corresponding reference process $\tilde{Z}$ is a d-dimensional RBM, which in combination with its d-dimensional pushing process $\tilde{Y}$ and the d-dimensional Brownian motion W defined in Section 2, satisfies

plus the obvious analogs of Equations (3) and (4). For the key identity (32) below, let

Remark 1.

In our numerical examples, the cost functions $c(z,\theta )$ allow for a closed-form expression of $F(z,x)$, simplifying our method. However, when the cost function is complex, computing $F(z,x)$ and its partial derivatives with respect to x (i.e., $\partial F(z,x)/\partial x_i$ ($i=1,\dots ,d$)) can be challenging. In such cases, several approaches are possible. First, one can generate a data set consisting of inputs ($z^l,x^l$) and outputs $F(z^l,x^l)$ for $l=1,\dots ,L$. With a sufficiently large L, a neural network can be trained offline to approximate F. The partial derivatives of F can then be approximated via autodifferentiation, as demonstrated in Ata and Zhou (2024). Alternatively, an optimization solver can be used to compute $F(z,x)$, and the envelope theorem can be invoked to determine its partial derivatives $\partial F(z,x)/\partial x_i$ ($i=1,\dots ,d$). However, this approach is computationally expensive. Lastly, the actor-critic method developed in Zhou et al. (2021a) can be employed.

Proposition 2.

If $V(·)$ is a $C^2$ function with polynomial growth that satisfies the HJB Equations (17) and (18), then it also satisfies the following identity almost surely for any $T>0$:

$$\begin{array}{l}{e}^{-rT}V(\tilde{Z}(T))-V(\tilde{Z}(0))={\displaystyle {\int }_0^T{e}^{-rt}}\nabla V(\tilde{Z}(t))·\text{d}W(t)\\ -{\displaystyle {\int }_0^T{e}^{-rt}}\kappa ·\text{d}\tilde{Y}(t)-{\displaystyle {\int }_0^T{e}^{-rt}}F(\tilde{Z}(t),\nabla V(\tilde{Z}(t)))\text{d}t.\end{array}$$(32)

Proof.

Applying Ito’s formula to $e^{-rt}V(\tilde{Z}(t))$ and using Equation (7) yield

$$\begin{array}{l}{e}^{-rT}V(\tilde{Z}(T))-V(\tilde{Z}(0))={\displaystyle {\int }_0^T{e}^{-rt}}\nabla V(\tilde{Z}(t))·\text{d}W(t)+{\displaystyle {\int }_0^T{e}^{-rt}}\mathcal{D}V(\tilde{Z}(t))·\text{d}\tilde{Y}(t)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\displaystyle {\int }_0^T{e}^{-rt}}(\mathcal{L}V(\tilde{Z}(t))-\tilde{\theta }·\nabla V(\tilde{Z}(t))-rV(\tilde{Z}(t)))\text{d}t.\end{array}$$(33)

Using the boundary Condition (18) plus the complementarity Condition (4) for $\tilde{Y}$ and $\tilde{Z}$, one has

$${\displaystyle {\int }_0^T{e}^{-rt}}\mathcal{D}V(\tilde{Z}(t))·\text{d}\tilde{Y}(t)=-{\displaystyle {\int }_0^T{e}^{-rt}}\kappa ·\text{d}\tilde{Y}(t).$$(34)

Furthermore, substituting $z=\tilde{Z}(t)$ in the HJB Equation (17), multiplying both sides by $e^{-rt}$, rearranging the terms, and integrating over [0, T] yields

$${\displaystyle {\int }_0^T{e}^{-rt}}(\mathcal{L}V(\tilde{Z}(t))-rV(\tilde{Z}(t)))\text{d}t={\displaystyle {\int }_0^T{e}^{-rt}}\underset{\theta \in \mathrm{\Theta }}{\max }(\theta ·\nabla V(\tilde{Z}(t))-c(\tilde{Z}(t),\theta ))\text{d}t.$$(35)

Substituting Equations (34) and (35) into Equation (33) gives Equation (32). □

Proposition 3 provides the motivation for the loss function that we strive to minimize in our computational method (see Section 5). Before developing that approach, we prove the following, which can be viewed as a converse of Proposition 3.

Proposition 3.

Suppose that $V:{\mathbb{R}}_{+}^d\to \mathbb{R}$ is a $C^2$ function; that $G:{\mathbb{R}}_{+}^d\to {\mathbb{R}}^d$ is continuous; and that $V,\hspace{0.17em}\nabla V$, and G all have polynomial growth. Also, assume that the following identity holds almost surely for some fixed $T>0$ and every $\tilde{Z}(0)=z\in {\mathbb{R}}_{+}^d$:

$$e^{-rT}V(\tilde{Z}(T))-V(\tilde{Z}(0))={\displaystyle {\int }_0^T{e}^{-rt}}G(\tilde{Z}(t))·\text{d}W(t)-{\displaystyle {\int }_0^T{e}^{-rt}}\kappa ·\text{d}\tilde{Y}(t)-{\displaystyle {\int }_0^T{e}^{-rt}}F(\tilde{Z}(t),G(\tilde{Z}(t)))\hspace{0.17em}\text{d}t.$$(36)

Then, $G(·)=\nabla V(·)$, and V satisfies the HJB Equations (17) and (18).

Remark 2.

The surprising conclusion that (36) implies $\nabla V(·)=G(·)$, without any a priori relationship between G and $\nabla V$ being assumed, motivates the “double-parametrization” method in Section 5.

Proof Sketch for Proposition 3.

Because $\tilde{Z}$ is a time-homogeneous Markov process, we can express (36) equivalently as follows for any $k=0,1,\dots$:

$$\begin{array}{l}{e}^{-rT}V(\tilde{Z}((k+1)T))-V(\tilde{Z}(kT))={\displaystyle {\int }_{kT}^{(k+1)T}{e}^{-r(t-kT)}}G(\tilde{Z}(t))·dW(t)-{\displaystyle {\int }_{kT}^{(k+1)T}{e}^{-r(t-kT)}}\kappa ·\text{d}\tilde{Y}(t)\\ -{\displaystyle {\int }_{kT}^{(k+1)T}{e}^{-r(t-kT)}}F(\tilde{Z}(t),G(\tilde{Z}(t)))\hspace{0.17em}dt.\end{array}$$(37)

Now, multiply both sides of (37) by $e^{-\mathit{\text{rkT}}}$, and then, add the resulting relationships for $k=0,1,\dots ,n-1$ to arrive at the following:

$$e^{-\mathit{\text{rnT}}}V(\tilde{Z}(nT))=V(\tilde{Z}(0))+{\displaystyle {\int }_0^{nT}{e}^{-rt}}G(\tilde{Z}(t))·dW(t)-{\displaystyle {\int }_0^{nT}{e}^{-rt}}\kappa ·\text{d}\tilde{Y}(t)-{\displaystyle {\int }_0^{nT}{e}^{-rt}}F(\tilde{Z}(t),G(\tilde{Z}(t)))dt.$$(38)

Because G has polynomial growth, one can show that

$${\mathbb{E}}_z\left[{\displaystyle {\int }_0^{nT}{e}^{-2rt}}G{(\tilde{Z}(t))}^{2}dz\right]<\infty$$

for all $n\ge 1$. Thus, when we take ${\mathbb{E}}_z$ of both sides of (38), the stochastic integral (that is, the second term) on the right side vanishes, and then, rearranging terms gives the following:

$$V(z)=e^{-\mathit{\text{rnT}}}{\mathbb{E}}_z\left[V(\tilde{Z}(nT))\right]+{\mathbb{E}}_z\left[{\displaystyle {\int }_0^{nT}{e}^{-rt}}F(\tilde{Z}(t),G(\tilde{Z}(t)))\text{d}t+{\displaystyle {\int }_0^{nT}{e}^{-rt}}\kappa ·\text{d}\tilde{Y}(t)\right],$$

for arbitrary positive integer n.

By Proposition 1 and the polynomial growth condition of V, we have $e^{-\mathit{\text{rnT}}}{\mathbb{E}}_z[V(\tilde{Z}(nT))]\to 0$ as $n\to \infty$. Therefore,

$$V(z)=\underset{n\to \infty }{\lim }{\mathbb{E}}_z\left[{\displaystyle {\int }_0^{nT}{e}^{-rt}}F(\tilde{Z}(t),G(\tilde{Z}(t)))\text{d}t+{\displaystyle {\int }_0^{nT}{e}^{-rt}}\kappa ·\text{d}\tilde{Y}(t)\right] \text{for} z\in {\mathbb{R}}_{+}^d.$$

Similarly, because F and G have polynomial growth, we conclude that

$$\begin{array}{l}{\mathbb{E}}_z\left[{\displaystyle {\int }_0^{\infty }{e}^{-rt}}|F(\tilde{Z}(t),G(\tilde{Z}(t)))|\text{d}t\right]<+\infty \text{for} z\in {\mathbb{R}}_{+}^d \text{and} \hfill \\ {\displaystyle {\int }_0^{nT}{e}^{-rt}}F(\tilde{Z}(t),G(\tilde{Z}(t)))\text{d}t\le {\displaystyle {\int }_0^{\infty }{e}^{-rt}}|F(\tilde{Z}(t),G(\tilde{Z}(t)))|\text{d}t<+\infty \text{for} z\in {\mathbb{R}}_{+}^d.\hfill \end{array}$$

Thus, by dominated convergence and monotone convergence, we have

$$V(z)={\mathbb{E}}_z\left[{\displaystyle {\int }_0^{\infty }{e}^{-rt}}F(\tilde{Z}(t),G(\tilde{Z}(t)))\text{d}t+{\displaystyle {\int }_0^{\infty }{e}^{-rt}}\kappa ·\text{d}\tilde{Y}(t)\right] \text{for} z\in {\mathbb{R}}_{+}^d.$$(39)

In other words, $V(z)$ can be viewed as the expected discounted cost associated with the RBM under the reference policy starting in state $\tilde{Z}(0)=z,$ where $F(·,G(·))$ is the state-cost function. Now, consider the following PDE:

$$\mathcal{L}f(z)-\tilde{\theta }·\nabla f(z)+F(z,G(z))=rf(z),\hspace{0.17em}z\in {\mathbb{R}}_{+}^d,$$

with boundary conditions

$${\mathcal{D}}_{i}f(z)=-{\kappa }_i \text{if} z_i=0\hspace{0.17em}(i=1,2,\dots ,d).$$

As assumed for the other PDEs considered in this paper, we assume that this PDE has a $C^2$ solution with polynomial growth. Using Ito’s lemma, one can then show that $f(z)=V(z)$. Thus, V satisfies the following PDE:

$$\mathcal{L}V(z)-\tilde{\theta }·\nabla V(z)+F(z,G(z))=rV(z),\hspace{0.17em}z\in {\mathbb{R}}_{+}^d,$$(40)

with boundary conditions

$${\mathcal{D}}_{i}V(z)=-{\kappa }_i \text{if} z_i=0\hspace{0.17em}(i=1,2,\dots ,d).$$(41)

Suppose that $G(·)=\nabla V(·)$ (which we will prove later). Substituting this into Equation (40) and using the definition of F, it follows that

$$\mathcal{L}V(z)-\underset{\theta \in \mathrm{\Theta }}{\max }\left\{\theta ·\nabla V(z)-c(z,\theta )\right\}=rV(z),\hspace{0.17em}z\in {\mathbb{R}}_{+}^d,$$(42)

which along with the boundary Condition (16), gives the desired result.

To complete the proof, it remains to show that $G(·)=\nabla V(·).$ By applying Ito’s formula to $e^{-rt}V(\tilde{Z}(t))$ and using Equations (3), (4), and (18), we conclude that

$$\begin{array}{l}{e}^{-rT}V(\tilde{Z}(T))-V(\tilde{Z}(0))={\displaystyle {\int }_0^T{e}^{-rt}}(\mathcal{L}V(\tilde{Z}(t))-\tilde{\theta }·\nabla V(\tilde{Z}(t))-rV(\tilde{Z}(t)))\text{d}t\\ \hspace{0.17em}+{\displaystyle {\int }_0^T{e}^{-rt}}\nabla V(\tilde{Z}(t))·\text{d}W(t)+{\displaystyle {\int }_0^T{e}^{-rt}}\mathcal{D}V(\tilde{Z}(t))·\text{d}Y(t).\end{array}$$

Then, using Equations (40) and (41), we rewrite the preceding equation as follows:

$$e^{-rT}V(\tilde{Z}(T))-V(\tilde{Z}(0))={\displaystyle {\int }_0^T{e}^{-rt}}\nabla V(\tilde{Z}(t))·\text{d}W(t)-{\displaystyle {\int }_0^T{e}^{-rt}}\kappa ·\text{d}\tilde{Y}(t)-{\displaystyle {\int }_0^T{e}^{-rt}}F(\tilde{Z}(t),G(\tilde{Z}(t)))\text{d}t.$$

Comparing this with Equation (36) yields

$${\displaystyle {\int }_0^T{e}^{-rt}}(G(\tilde{Z}(t))-\nabla V(\tilde{Z}(t)))·\text{d}W(t)=0,$$

which yields the following:

$${\mathbb{E}}_z\left[{\left({\displaystyle {\int }_0^T{e}^{-rt}}(G(\tilde{Z}(t))-\nabla V(\tilde{Z}(t)))·\text{d}W(t)\right)}^2\right]=0.$$(43)

Thus, provided that $e^{-rt}(G(\tilde{Z}(t))-\nabla V(\tilde{Z}(t)))$ is square integrable, Ito’s isometry (Zhang et al. 2020, lemma D.1) yields the following:

$${\mathbb{E}}_z\left[{\left({\displaystyle {\int }_0^T{e}^{-rt}}(G(\tilde{Z}(t))-\nabla V(\tilde{Z}(t)))·\text{d}W(t)\right)}^2\right]={\mathbb{E}}_z\left[{\displaystyle {\int }_0^T{\Vert e^{-rt}(G(\tilde{Z}(t))-\nabla V(\tilde{Z}(t)))\Vert }_A^2}\text{d}t\right]=0,$$

where $\Vert x{\Vert }_A≔x^{\top }Ax$. The square integrability of $e^{-rt}(G(\tilde{Z}(t))-\nabla V(\tilde{Z}(t)))$ follows because G and $\nabla V$ have polynomial growth and the action space $\mathrm{\Theta }$ is bounded. Because A is a positive definite matrix, we then have $\nabla V(\tilde{Z}(t))=G(\tilde{Z}(t))$ almost surely. By the continuity of $\nabla V(·)$ and $G(·),$ we conclude that $\nabla V(·)=G(·)$. □

4.2. Ergodic Control

Again, we use a reference policy with the constant (negative) drift vector $\tilde{\theta }$, and now, we assume that $R^{-1}\tilde{\theta }>0$, which ensures that the reference policy is admissible for our ergodic control formulation.

Proposition 4.

If $v(·)$ and $\xi$ solve the HJB Equations (26) and (27) and $v(·)$ is a $C^2$ function with polynomial growth, then they also satisfy the following identity almost surely for any $T>0$:

$$v(\tilde{Z}(T))-v(\tilde{Z}(0))={\displaystyle {\int }_0^T\nabla }v(\tilde{Z}(t))·\text{d}W(t)+T\xi -{\displaystyle {\int }_0^T\kappa }·\text{d}\tilde{Y}(t)-{\displaystyle {\int }_0^{T}F}(\tilde{Z}(t),\nabla v(\tilde{Z}(t)))\hspace{0.17em}\text{d}t.$$(44)

Proof.

Applying Ito’s formula to $v(z)$ yields

$$v(\tilde{Z}(T))-v(\tilde{Z}(0))\begin{array}{l}={\displaystyle {\int }_0^T\nabla }v(\tilde{Z}(t))·\text{d}W(t)+{\displaystyle {\int }_0^T\mathcal{D}}v(\tilde{Z}(t))·\text{d}\tilde{Y}(t)+{\displaystyle {\int }_0^T(\mathcal{L}v(\tilde{Z}(t))-\tilde{\theta }·\nabla v(\tilde{Z}(t)))}\text{d}t.\hfill \end{array}$$(45)

Recall that the boundary condition of the HJB equation is ${\mathcal{D}}_{j}v(z)=-{\kappa }_j$ if $z_j=0$. Thus, Equations (3) and (4) jointly imply

$${\displaystyle {\int }_0^T\mathcal{D}}v(\tilde{Z}(t))·\text{d}\tilde{Y}(t)=-{\displaystyle {\int }_0^T\kappa }·\text{d}\tilde{Y}(t).$$

Then, substituting the HJB Equation (26) into Equation (45) gives (44). □

Proposition 5.

Suppose that $v:{\mathbb{R}}_{+}^d\to \mathbb{R}$ is a $C^2$ function; that $g:{\mathbb{R}}_{+}^d\to {\mathbb{R}}^d$ is continuous; and that $v,\hspace{0.17em}\nabla v,\hspace{0.17em}g$ all have polynomial growth. Also, assume that the following identity holds almost surely for some fixed $T>0$, a scalar $\xi$, and every $Z(0)=z\in {\mathbb{R}}_{+}^d$:

$$v(\tilde{Z}(T))-v(\tilde{Z}(0))={\displaystyle {\int }_0^{T}g}(\tilde{Z}(t))·\text{d}W(t)+T\xi -{\displaystyle {\int }_0^T\kappa }·\text{d}\tilde{Y}(t)-{\displaystyle {\int }_0^{T}F}(\tilde{Z}(t),g(\tilde{Z}(t)))\hspace{0.17em}\text{d}t.$$(46)

Then, $g(·)=\nabla v(·)$, and $(v,\xi )$ satisfies the HJB Equations (17) and (18).

Proof Sketch.

Let $\tilde{\pi }$ be the stationary distribution of the RBM $\tilde{Z}$ under the reference policy and $\tilde{Z}(\infty )$ be a random variable with the distribution $\tilde{\pi }$. Then, assuming that the initial distribution of the RBM under the reference policy is $\tilde{\pi }$ (i.e., $\tilde{Z}(0)\sim \tilde{\pi }$), its marginal distribution at time t is also $\tilde{\pi }$ (i.e., $\tilde{Z}(t)\sim \tilde{\pi }$ for every $t\ge 0$).

Because g has polynomial growth, one can show that the expectation of the stochastic integral (that is, the first term) on the right side of (46) vanishes. Then, by taking the expectation over $\tilde{Z}(0)\sim \tilde{\pi },$Equation (46) implies

$${\mathbb{E}}_{\tilde{\pi }}\left[v(\tilde{Z}(0))\right]={\mathbb{E}}_{\tilde{\pi }}\left[v(\tilde{Z}(T))\right]+{\mathbb{E}}_{\tilde{\pi }}\left[{\displaystyle {\int }_0^{T}F}(\tilde{Z}(t),g(\tilde{Z}(t)))\text{d}t\right]+{\mathbb{E}}_{\tilde{\pi }}\left[{\displaystyle {\int }_0^T\kappa }·\text{d}\tilde{Y}(t)\right]-T\xi .$$(47)

By observing that ${\mathbb{E}}_{\tilde{\pi }}[v(\tilde{Z}(0))]={\mathbb{E}}_{\tilde{\pi }}[v(\tilde{Z}(T))]$ and

$$\begin{array}{l}{\mathbb{E}}_{\tilde{\pi }}\left[F(\tilde{Z}(t),g(\tilde{Z}(t)))\right]=\mathbb{E}\left[F(\tilde{Z}(\infty ),g(\tilde{Z}(\infty )))\right] \text{for} t\ge 0,\hfill \\ {\mathbb{E}}_{\tilde{\pi }}\left[{\displaystyle {\int }_0^T\kappa }·\text{d}\tilde{Y}(t)\right]=T{\mathbb{E}}_{\tilde{\pi }}\left[{\displaystyle {\int }_0^1\kappa }·\text{d}\tilde{Y}(t)\right],\hfill \end{array}$$

we conclude that

$$\begin{array}{l}\xi =\mathbb{E}[F(\tilde{Z}(\infty ),g(\tilde{Z}(\infty )))]+{\mathbb{E}}_{\tilde{\pi }}\left[{\displaystyle {\int }_0^1\kappa }·\text{d}\tilde{Y}(t)\right]\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\mathbb{E}[F(\tilde{Z}(\infty ),g(\tilde{Z}(\infty )))]+{\displaystyle \sum _{i=1}^d{\kappa }_i}{\tilde{\nu }}_i(S_i),\end{array}$$

where ${\tilde{\nu }}_i(·)$ is the boundary measure on the boundary surface $S_i=\{z\in {\mathbb{R}}_{+}^d:z_i=0\}$ associated with the RBM under the reference policy. In other words, $\xi$ can be viewed as the expected steady-state cost associated with the RBM under the reference policy, where $F(·,g(·))$ is the state-cost function. Now, consider the following PDE:

$$\mathcal{L}\tilde{v}(z)-\tilde{\theta }·\nabla \tilde{v}(z)+F(z,g(z))=\xi ,\hspace{0.17em}\hspace{0.17em}z\in {\mathbb{R}}_{+}^d,$$(48)

with boundary conditions ${\mathcal{D}}_i\tilde{v}(z)=-{\kappa }_i \text{if} z_i=0\hspace{0.17em}\hspace{0.17em}(i=1,\dots ,d)$.

We assume that (48) has a $C^2$ solution with polynomial growth. Using Ito’s lemma, one can then show that $\tilde{v}$ is the relative value function corresponding to the reference process $\tilde{Z}$ (i.e., for policy $u(z)=\tilde{\theta }$ for $z\in {\mathbb{R}}_{+}^d$) under the state cost function $F(z,g(z))$. Furthermore, applying Ito’s formula to $\tilde{v}(\tilde{Z}(t))$ on the interval [0, nT] for $n\ge 1$ yields

$$\tilde{v}(\tilde{Z}(nT))-\tilde{v}(\tilde{Z}(0))\begin{array}{l}={\displaystyle {\int }_0^{nT}\nabla }\tilde{v}(\tilde{Z}(t))·\text{d}W(t)+{\displaystyle {\int }_0^{nT}\mathcal{D}}\tilde{v}(\tilde{Z}(t))·\text{d}\tilde{Y}(t)+{\displaystyle {\int }_0^{nT}(\mathcal{L}\tilde{v}(\tilde{Z}(t))-\tilde{\theta }·\nabla \tilde{v}(\tilde{Z}(t)))}\text{d}t.\hfill \end{array}$$(49)

Because $\tilde{v}(z)$ also satisfies the boundary conditions ${\mathcal{D}}_i\tilde{v}(z)=-{\kappa }_i \text{if} z_i=0\hspace{0.17em}\hspace{0.17em}(i=1,\dots ,d)$, it follows from Equations (3) and (4) that

$${\displaystyle {\int }_0^{nT}\mathcal{D}}\tilde{v}(\tilde{Z}(t))·\text{d}\tilde{Y}(t)=-{\displaystyle {\int }_0^{nT}\kappa }·\text{d}\tilde{Y}(t).$$

Then, substituting Equation (48) into Equation (49) gives

$$\tilde{v}(\tilde{Z}(nT))-\tilde{v}(\tilde{Z}(0))={\displaystyle {\int }_0^{nT}\nabla }\tilde{v}(\tilde{Z}(t))·\text{d}W(t)+nT\xi -{\displaystyle {\int }_0^{nT}\kappa }·\text{d}\tilde{Y}(t)-{\displaystyle {\int }_0^{nT}F}(\tilde{Z}(t),g(\tilde{Z}(t)))\hspace{0.17em}\text{d}t,\hspace{0.17em}n\ge 1.$$(50)

In the proof of Proposition 3, we first showed that because $\tilde{Z}$ is a time-homogeneous Markov process, the assumed stochastic relationship (36) can be extended to the more general form (38), with n being an arbitrary positive integer. In the current context, one can argue in exactly the same way to establish the following. First, the assumed stochastic relationship (46) actually holds in the more general form where T is replaced by nT, with n being an arbitrary positive integer. Then, after taking expectations on both sides of the generalized version of (46) and (50), we arrive at the following:

$$v(z)={\mathbb{E}}_z\left[v(\tilde{Z}(nT))\right]-nT\xi +{\mathbb{E}}_z\left[{\displaystyle {\int }_0^{nT}\kappa }·\text{d}\tilde{Y}(t)\right]+{\mathbb{E}}_z\left[{\displaystyle {\int }_0^{nT}F}(\tilde{Z}(t),g(\tilde{Z}(t)))\text{d}t\right] \text{and}$$(51)

$$\begin{array}{l}\tilde{v}(z)={\mathbb{E}}_z\left[\tilde{v}(\tilde{Z}(nT))\right]-nT\xi +{\mathbb{E}}_z\left[{\displaystyle {\int }_0^{nT}\kappa }·\text{d}\tilde{Y}(t)\right]+{\mathbb{E}}_z\left[{\displaystyle {\int }_0^{nT}F}(\tilde{Z}(t),g(\tilde{Z}(t)))\text{d}t\right],\hfill \end{array}$$(52)

for $z\in {\mathbb{R}}_{+}^d$ and an arbitrary positive integer n. Note that the expectation of the stochastic integral vanishes because $\nabla \tilde{v}$ has polynomial growth. Subtracting (52) from (51) further yields

$$v(z)-\tilde{v}(z)={\mathbb{E}}_z\left[v(Z(nT))\right]-{\mathbb{E}}_z\left[\tilde{v}(Z(nT))\right].$$

Because the solution of (48) is determined only up to an additive constant, we assume $\mathbb{E}[\tilde{v}(\tilde{Z}(\infty ))]=\mathbb{E}[v(\tilde{Z}(\infty ))]$. Because $v(·)$ and $\tilde{v}(·)$ have polynomial growth, we conclude the following from Budhiraja and Lee (2007, theorem 4.12):

$$\begin{array}{ll}\underset{n\to +\infty }{\lim }{\mathbb{E}}_z\left[v(Z(nT))\right]\hfill & =\mathbb{E}\left[v(\tilde{Z}(\infty ))\right] \text{and}\hfill \\ \underset{n\to +\infty }{\lim }{\mathbb{E}}_z\left[\tilde{v}(Z(nT))\right]\hfill & =\mathbb{E}\left[\tilde{v}(\tilde{Z}(\infty ))\right].\hfill \end{array}$$

Therefore, we have

$$v(z)-\tilde{v}(z)=\underset{n\to +\infty }{\lim }({\mathbb{E}}_z\left[v(Z(nT))\right]-{\mathbb{E}}_z\left[\tilde{v}(Z(nT))\right])=0 \text{for} z\in {\mathbb{R}}_{+}^d,$$

which means that $v(·)$ also satisfies the PDE (48) and the associated boundary conditions. That is,

$$\mathcal{L}v(z)-\tilde{\theta }·\nabla v(z)+F(z,g(z))=\xi \text{for} z\in {\mathbb{R}}_{+}^d.$$(53)

$$\begin{array}{l}{\mathcal{D}}_{i}v(z)=-{\kappa }_i \text{if} z_i=0\hspace{0.17em}\hspace{0.17em}(i=1,\dots ,d).\hfill \end{array}$$(54)

Suppose that $g(·)=\nabla v(·)$ (which we will prove later). Substituting this into Equation (53) and using the definition of F, it follows that

$$\mathcal{L}v(z)-\underset{\theta \in \mathrm{\Theta }}{\max }\left\{\theta ·\nabla v(z)-c(z,\theta )\right\}=\xi ,\hspace{0.17em}\hspace{0.17em}z\in {\mathbb{R}}_{+}^d,$$

which along with the boundary Condition (54), gives the desired result.

To complete the proof, it remains to show that $g(·)=\nabla v(·).$ By applying Ito’s formula to $v(\tilde{Z}(t))$ and using Equations (3), (4), and (54), we conclude that

$$v(\tilde{Z}(T))-v(\tilde{Z}(0))={\displaystyle {\int }_0^T(\mathcal{L}v(\tilde{Z}(t))-\tilde{\theta }·\nabla v(\tilde{Z}(t)))}\text{d}t-{\displaystyle {\int }_0^T\kappa }·\text{d}\tilde{Y}(t)+{\displaystyle {\int }_0^T\nabla }v(\tilde{Z}(t))·\text{d}W(t).$$

Then, using Equation (53), we rewrite the preceding equation as follows:

$$v(\tilde{Z}(T))-v(\tilde{Z}(0))=T\xi -{\displaystyle {\int }_0^{T}F}(\tilde{Z}(t),g(\tilde{Z}(t)))\text{d}t-{\displaystyle {\int }_0^T\kappa }·\text{d}\tilde{Y}(t)+{\displaystyle {\int }_0^T\nabla }v(\tilde{Z}(t))·\text{d}W(t).$$

Comparing this with Equation (46) yields

$${\displaystyle {\int }_0^T(g(\tilde{Z}(t))-\nabla v(\tilde{Z}(t)))}·\text{d}W(t)=0,$$

which yields the following:

$${\mathbb{E}}_z\left[{\left({\displaystyle {\int }_0^T(g(\tilde{Z}(t))-\nabla v(\tilde{Z}(t)))}·\text{d}W(t)\right)}^2\right]=0.$$(55)

Thus, provided that $g(\tilde{Z}(t))-\nabla v(\tilde{Z}(t))$ is square integrable, Ito’s isometry (Zhang et al. 2020, lemma D.1) gives the following:

$${\mathbb{E}}_z\left[{\left({\displaystyle {\int }_0^T(g(\tilde{Z}(t))-\nabla v(\tilde{Z}(t)))}·\text{d}W(t)\right)}^2\right]={\mathbb{E}}_z\left[{\displaystyle {\int }_0^T{\Vert g(\tilde{Z}(t))-\nabla v(\tilde{Z}(t))\Vert }_A^2}\text{d}t\right]=0.$$

The square integrability of $g(\tilde{Z}(t))-\nabla v(\tilde{Z}(t))$ follows because g and $\nabla v$ have polynomial growth, and ${\mathbb{E}}_z(|\tilde{Z}(nT){|}^k)$ is finite for all k because our action space $\mathrm{\Theta }$ is bounded. Then, because A is positive definite, $\nabla v(\tilde{Z}(t))=g(\tilde{Z}(t))$ almost surely. By the continuity of $\nabla v(·)$ and $g(·),$ we conclude that $\nabla v(·)=g(·)$. □

5. Computational Method

We follow in the footsteps of Han et al. (2018), who developed a computational method to solve semilinear parabolic partial differential equations. Those authors focused on a backward stochastic differential equation associated with their PDE, and in similar fashion, we focus on the stochastic differential Equations (36) and (46) that are associated with our two stochastic control formulations (see Section 4). Our method differs from that of Han et al. (2018) because they consider PDEs on a finite time interval with an unbounded state space and a specified terminal condition, whereas our stochastic control problem has an infinite time horizon and state-space constraints. As such, it leads to a PDE on a polyhedral domain with oblique derivative boundary conditions. We modify the approach of Han et al. (2018) to incorporate those additional features, treating the discounted and ergodic formulations in Sections 5.1 and 5.2, respectively.

5.1. Discounted Control

We approximate the value function $V(·)$ and its gradient $\nabla V(·)$ by the deep neural networks $V_{w_1}(·)$ and $G_{w_2}(·)$, respectively, with associated parameter vectors $w_1$ and $w_2$. Seeking an approximate solution of the stochastic Equation (36), we define the loss function

The initial state $\tilde{Z}(0)$ can be randomly chosen, and the expectation in the previous equation is calculated with respect to the sample path distribution of the reference process $\tilde{Z}(·)$ (see Algorithm 1 for details). As mentioned in step 6 of Algorithm 1, the process $\tilde{Z}(·)$ is run continuously, meaning that the terminal state of one iteration becomes the initial state of the next iteration. This approach can be seen as an approximation to starting the reference process with its steady-state distribution.

Remark 3.

The loss function $\mathcal{\ell }(w_1,w_2)$ defined in Equation (56) corresponds to the $L_1(V)$ loss defined in Zhou et al. (2021a, equation 2.31). The authors consider this a “variance reduced” loss function because of the additional term ${\int }_0^T{e}^{-rt}{G}_{w_2}(\tilde{Z}(t))·\text{d}W(t)$. This term can be interpreted as an approximating martingale process or a control variate, a common approach in the simulation literature (see, for example, Andradóttir et al. 1993, Henderson and Glynn 2002, Dai and Gluzman 2022).

Our definition (56) of the loss function does not explicitly enforce the consistency requirement $\nabla V_{w_1}(·)=G_{w_2}(·)$, but Proposition 3 provides the justification for this separate parametrization. This type of double parametrization has also been implemented by Zhou et al. (2021a).

Subroutine 1

(Euler Discretization Scheme)

• 1: Input: The drift vector $-\tilde{\theta }$, the covariance matrix A, the reflection matrix R, the time horizon T, a step size h (for simplicity, we assume that $N\triangleq T/h$ is an integer), and a starting point $\tilde{Z}(0)=z$.

• 2: Output: A discretized reflected Brownian motion with the boundary pushing process increments and the Brownian increments at times $h,2h,\dots ,Nh$.

• 3: function Discretize(T, h, z)

• 4:  For time interval [0, T] and $N=T/h$, construct the partition $0=t_0<t_1<\cdots <t_N=T$, where $\mathrm{\Delta }{t}_n=t_{n+1}-t_n=h$ for $n=0,1,\dots ,N-1$.

• 5:  Generate N i.i.d. d-dimensional Gaussian random variables with mean of zero and covariance matrix hA, denoted by ${\delta }_0,\dots ,{\delta }_{N-1}$.

• 6:  for $k\leftarrow 0$ to $N-1$, do

• 7:   $x\leftarrow \tilde{Z}(kh)+{\delta }_k-\tilde{\theta }h$

• 8:   $\tilde{Z}((k+1)h),\mathrm{\Delta }\tilde{Y}(kh)$$\leftarrow$ Skorokhod(x)

• 9:  end for

• 10:  return $\tilde{Z}(h),\tilde{Z}(2h),\dots ,\tilde{Z}(Nh)$; $\mathrm{\Delta }\tilde{Y}(0)$, $\mathrm{\Delta }\tilde{Y}(h),\dots ,\mathrm{\Delta }\tilde{Y}((N-1)h)$; and ${\delta }_0$,…,${\delta }_{N-1}$.

• 11: end function

Subroutine 2

(Solve the Skorokhod Problem (Linear Complementarity Problem))

• 1: Input: A vector $x\in {\mathbb{R}}^d$ and the reflection matrix R.

• 2: Output: A solution to the Skorokhod problem $y\in {\mathbb{R}}_{+}^d$

• 3: Set $ϵ={10}^{-8}$;

• 4: function Skorokhod(x)

• 5:  $y=x$;

• 6:  $u=0$;

• 7:  while exists $y_i<-ϵ$, do

• 8:   Compute the set $B=\{i:y_i<ϵ$};

• 9:   Compute $L_B=-R_{B,B}^{-1}{x}_B$;

• 10:   Compute $y=x+R_{:,B}\times L_B$;

• 11:  end while

• 12:  $u_B=L_B$;

• 13:  return y, u.

• 14: end function

Our computational method seeks a neural network parameter combination $(w_1,w_2)$ that minimizes an approximation of the loss defined via (56). Specifically, we first simulate multiple discretized paths of the reference RBM $\tilde{Z}$ with the boundary pushing process $\tilde{Y}$, restricted to a fixed and finite time domain [0, T]. To do that, we sample discretized paths of the underlying Brownian motion W, and then, we solve a discretized Skorohod problem for each path of W (this is the purpose of Subroutine 2) to obtain the corresponding path of $\{\tilde{Z},\tilde{Y}\}$. Thereafter, our method computes a discretized version of the loss (56), summing over sampled paths to approximate the expectation and over discrete time steps to approximate the integral over [0, T], and our method minimizes it using stochastic gradient descent; see Algorithm 1. The discretization inevitably introduces bias into our method. Although we do not provide a formal convergence proof as the partition becomes finer, we direct readers to Han and Long (2020) for a rigorous analysis in a similar context. Furthermore, our numerical examples indicate that the impact of discretization on the learned control is small.

In Subroutine 2, given the index set B, $R_{B,B}$ is the submatrix derived by deleting the rows and columns of R with indices in $\{1,\dots ,d\}$\B. Similarly, $R_{:,B}$ is the matrix that one arrives at by deleting the columns of R whose indices are in the set $\{1,\dots ,d\}$\B. One has considerable latitude in choosing $\tilde{\theta }$ (i.e., the reference policy), provided that it can explore the state space sufficiently. For our numerical examples, we set $\tilde{\theta }=1$.

Algorithm 1

(Method for the Discounted Control Case)

• 1: Input: The number of iteration steps M, a batch size B, a learning rate $\alpha$, a time horizon T, a discretization step size h (for simplicity, we assume that $N\triangleq T/h$ is an integer), a starting point z, and an optimization solver (SGD, ADAM, RMSProp, etc.).

• 2: Output: A neural network approximation of the value function $V_{w_1}$ and the gradient function $G_{w_2}$.

• 3: Initialize the neural networks $V_{w_1}$ and $G_{w_2}$; set $z_0^{(i)}=z$ for $i=1,2,\dots ,B$.

• 4: for $k\leftarrow 0$ to $M-1$, do

• 5:  Simulate B discretized RBM paths and the Brownian increments $\{{\tilde{Z}}^{(i)},\mathrm{\Delta }{\tilde{Y}}^{(i)},{\delta }^{(i)}\}$ with a time horizon T and a discretization step size h starting from ${\tilde{Z}}^{(i)}(0)=z_k^{(i)}$ by invoking Discretize($T,h,z_k^{(i)}$) for $i=1,2,\dots ,B$.

• 6:  Compute the empirical loss

  $$\begin{array}{l}\widehat{\mathcal{\ell }}(w_1,w_2)=\frac{1}{B}{\displaystyle \sum _{i=1}^B(e^{-rT}{V}_{w_1}({\tilde{Z}}^{(i)}(T))-V_{w_1}({\tilde{Z}}^{(i)}(0))+{\displaystyle \sum _{j=0}^{N-1}{e}^{-\mathit{\text{rhj}}}}\kappa ·\mathrm{\Delta }{\tilde{Y}}^{(i)}(hj)}\\ {-{\displaystyle \sum _{j=0}^{N-1}{e}^{-\mathit{\text{rhj}}}}{G}_{w_2}({\tilde{Z}}^{(i)}(hj))·{\delta }_j^{(i)}+{\displaystyle \sum _{j=0}^{N-1}{e}^{-\mathit{\text{rhj}}}}F({\tilde{Z}}^{(i)}(hj),G_{w_2}({\tilde{Z}}^{(i)}(hj)))h)}^2.\end{array}$$(57)

• 7:  Compute the gradient $\partial \widehat{\mathcal{\ell }}(w_1,w_2)/\partial w_1,\partial \widehat{\mathcal{\ell }}(w_1,w_2)/\partial w_2$, and update $w_1,w_2$ using the chosen optimization solver.

• 8:  Update $z_{k+1}^{(i)}$ as the end point of the path ${\tilde{Z}}^{(i)}$: $z_{k+1}^{(i)}\leftarrow {\tilde{Z}}^{(i)}(T)$.

• 9: end for

• 10: return functions $V_{w_1}(·)$ and $G_{w_2}(·)$.

After the parameter values $w_1$ and $w_2$ have been determined, our proposed policy is as follows:

Remark 4.

One can also consider the policy using $\nabla V_{w_1}(·)$ instead of $G_{w_2}(·).$ That is,

$$\arg \underset{\theta \in \mathrm{\Theta }}{\max }\left\{\theta ·\nabla V_{w_1}(z)-c(z,\theta )\right\},\hspace{0.17em}z\in {\mathbb{R}}_{+}^d.$$(59)

However, our numerical experiments suggest that this policy is inferior to (58).

5.2. Ergodic Control

We parametrize $v(·)$ and $\nabla v(·)$ using deep neural networks $v_{w_1}(·)$ and $g_{w_2}(·)$ with parameters $w_1$ and $w_2$, respectively, and then, we use Equation (46) to define an auxiliary loss function:

Then, we define the loss function $\mathcal{\ell }(w_1,w_2)={\min }_{\xi }\tilde{\mathcal{\ell }}(w_1,w_2,\xi )$. Letting X denote a random variable with a finite second moment, we note that

Thus, we arrive at the following expression for the loss function

We present our method for the ergodic control case formally in Algorithm 2.

Algorithm 2

(Method for the Ergodic Control Case)

• 1: Input: The number of iteration steps M, a batch size B, a learning rate $\alpha$, a time horizon T, a discretization step size h (for simplicity, we assume that $N\triangleq T/h$ is an integer), a starting point z, and an optimization solver (SGD, ADAM, RMSProp, etc.).

• 2: Output: A neural network approximation of the value function $v_{w_1}$ and the gradient function $g_{w_2}$.

• 3: Initialize the neural networks $v_{w_1}$ and $g_{w_2}$; set $z_0^{(i)}=z$ for $i=1,2,\dots ,B$.

• 4: for $k\leftarrow 0$ to $M-1$, do

• 5:  Simulate B discretized RBM paths and the Brownian increments $\{{\tilde{Z}}^{(i)},\mathrm{\Delta }{\tilde{Y}}^{(i)},{\delta }^{(i)}\}$ with a time horizon T and a discretization step size h starting from ${\tilde{Z}}^{(i)}(0)=z_k^{(i)}$ by invoking Discretize($T,h,z_k^{(i)}$), for $i=1,2,\dots ,B$.

• 6:  Compute the empirical loss

  $$\begin{array}{l}\widehat{\mathcal{\ell }}(w_1,w_2)=\widehat{\text{Var}}(v_{w_1}({\tilde{Z}}^{(i)}(T))-v_{w_1}({\tilde{Z}}^{(i)}(0))-{\displaystyle \sum _{j=0}^{N-1}{g}_{w_2}}({\tilde{Z}}^{(i)}(hj))·{\delta }_j^{(i)}\\ +{\displaystyle \sum _{j=0}^{N-1}\kappa }·\mathrm{\Delta }{\tilde{Y}}^{(i)}(hj)+{\displaystyle \sum _{j=0}^{N-1}f}({\tilde{Z}}^{(i)}(hj),g_{w_2}({\tilde{Z}}^{(i)}(hj)))h).\end{array}$$(62)

• 7:  Compute the gradient $\partial \widehat{\mathcal{\ell }}(w_1,w_2)/\partial w_1,\partial \widehat{\mathcal{\ell }}(w_1,w_2)/\partial w_2$, and update $w_1,w_2$ using the chosen optimization solver.

• 8:  Update $z_{k+1}^{(i)}$ as the end point of the path ${\tilde{Z}}^{(i)}$: $z_{k+1}^{(i)}\leftarrow {\tilde{Z}}^{(i)}(T)$.

• 9: end for

• 10: return functions $v_{w_1}(·)$ and $g_{w_2}(·)$.

After the parameters values $w_1$ and $w_2$ have been determined, our proposed policy is the following:

6. Three Families of Test Problems

Here, we specify three families of test problems for which numerical results will be presented later (see Section 7). Each family consists of RBM drift control problems indexed by $d=1,2,\dots ,$ where d is the dimension of the orthant that serves as the problem’s state space. The first of the three problem families, specified in Section 6.1, is characterized by a feed-forward network structure and linear cost of control. Recapitulating the earlier work by Ata (2006), Section 6.2 explains the interpretation of such problems as “heavy traffic” limits of input control problems for certain feed-forward queueing networks.

Our second family of test problems is identical to the first one except that now the cost of control is quadratic rather than linear. The exact meaning of that phrase will be spelled out in Section 6.3, where we also explain the interpretation of such problems as heavy traffic limits of dynamic pricing problems for queueing networks. In Section 6.4, we describe two parametric families of policies with special structure that will be used later for comparison purposes in our numerical study. Finally, Section 6.5 specifies our third family of test problems, which have a separable structure that allows them to be solved exactly by analytical means. Such problems are of obvious value for evaluating the accuracy of our computational method. For all of our test problems, the penalty rates associated with pushing at the boundary are set to zero. That is, $\kappa =0$.

6.1. Main Example with Linear Cost of Control

We consider a family of test problems with parameters $K=0,1,\dots$, attaching to each such a problem as the index d (mnemonic for dimension) $=K+1$. Problem d has state space ${\mathbb{R}}_{+}^d$ and the $d\times d$ reflection matrix

where $p_1,\dots ,p_K>0$ and $p_1+\cdots +p_K=1$. Also, the set of drift vectors available in each state iswhere the lower limit ${\underset{¯}{\theta }}_k$ and the upper limit ${\overline{\theta }}_k$ are as specified in Section 6.2 below. Similarly, the $d\times d$ covariance matrix A for problem d is as specified in Section 6.2. Finally, the cost function for problem d has the linear form

That is, the cost rate $c(Z(t),u(Z(t)))$ that the system manager incurs under policy u at time t is linear in both the state vector $Z(t)$ and the chosen drift rate $u(Z(t))$.

In either the discounted control setting or the ergodic control setting, inspection of the HJB equation displayed earlier in Section 3 shows that given this linear cost structure, there exists an optimal policy $u^{*}(·)$ such that

for each state $z\in {\mathbb{R}}_{+}^{K+1}$ and each component $k=0,1,\dots ,K$.

In the next section, we explain how drift control problems of the form specified here arise as heavy traffic limits in queueing theory. Strictly speaking, however, that interpretation of the test problems is inessential to the main subject of this paper; the computational results presented in Section 7 can be read without reference to the queueing-theoretic interpretations of our test problems.

6.2. Interpretation as Heavy Traffic Limits of Queueing Network Control Problems

Let us consider the feed-forward queueing network model of the make-to-order production system portrayed in Figure 1. There are $d=K+1$ buffers, represented by the open-ended rectangles in Figure 1, indexed by $k=0,1,\dots ,K.$ Each buffer has a dedicated server, represented by the circles in Figure 1. Arriving jobs wait in their designated buffer if the server is busy. There are two types of jobs arriving to the system: regular versus thin streams. Thin-stream jobs have the same service time distributions as the regular jobs, but they differ from the regular jobs in two important ways. First, thin-stream jobs can be turned away upon arrival. That is, a system manager can exercise admission control in this manner, but in contrast, she must admit all regular jobs arriving to the system. Second, the volume of thin-stream jobs is smaller than that of the regular jobs; see Assumption 1.

Regular jobs enter the systems only through buffer 0, as shown by the solid arrow pointing to buffer 0 in Figure 1. A renewal process $E=\{E(t):t\ge 0\}$ models the cumulative number of regular jobs arriving to the system over time. We let $\lambda$ denote the arrival rate and $a^2$ denote the squared coefficient of variation of the interarrival times for the regular jobs. The thin-stream jobs arrive to buffer k (as shown by the dashed arrows in Figure 1) according to the renewal process $A_k=\{A_k(t):t\ge 0\}$ for $k=0,1,\dots ,K$. We let ${\eta }_k$ denote the arrival rate and $b_k^2$ denote the squared coefficient of variation of the interarrival times for renewal process $A_k$.

Jobs in buffer k have i.i.d. general service time distributions with mean $m_k$ and squared coefficient of variation $s_k^2\ge 0,\hspace{0.17em}k=0,1,\dots ,K$; ${\mu }_k=1/m_k$ is the corresponding service rate. We let $S_k=\{S_k(t):t\ge 0\}$ denote the renewal process associated with the service completions by server k for $k=1,\dots ,K$. To be specific, $S_k(t)$ denotes the number of jobs server k processes by time t if it incurs no idleness during [0, t]. The jobs in each buffer are served on a first come, first served basis, and servers work continuously unless their buffer is empty. After receiving service, jobs in buffer 0 join buffer k with probability $p_k,\hspace{0.17em}k=1,2,\dots ,K$, independently of other events. This probabilistic routing structure is captured by a vector-valued process $\mathrm{\Phi }(·)$, where ${\mathrm{\Phi }}_k(\mathcal{\ell })$ denotes the total number of jobs routed to buffer k among the first $\mathcal{\ell }$ jobs served by server 0 for $k=1,\dots ,K$ and $\mathcal{\ell }\ge 1$. We let $p=(p_k)$ denote the K-dimensional vector of routing probabilities. Jobs in buffers $1,\dots ,K$ leave the system upon receiving service.

As stated earlier, the system manager makes admission control decisions for thin-stream jobs. Turning away a thin-stream job arriving to buffer k (externally) results in a penalty of $c_k$. For mathematical convenience, we model admission control decisions as if the system manager can simply “turn off” each of the thin-stream arrival processes as desired. In particular, we let ${\mathrm{\Delta }}_k(t)$ denote the cumulative amount of time that the (external) thin-stream input to buffer k is turned off during the interval [0, t]. Thus, the vector-valued process $\mathrm{\Delta }=({\mathrm{\Delta }}_k)$ represents the admission control policy. Similarly, we let $T_k(t)$ denote the cumulative amount of time that server k is busy during the time interval [0, t], and $I_k(t)=t-T_k(t)$ denotes the cumulative amount of idleness that server K incurs during [0, t].

Letting $Q_k(t)$ denote the number of jobs in buffer k at time t, the vector-valued process $Q=(Q_k)$ will be called the queue-length process. Given a control $\mathrm{\Delta }=({\mathrm{\Delta }}_k)$, assuming $Q(0)=0$, it follows that

Moreover, the following must hold:

The system manager also incurs a holding cost at rate $h_k$ per job in buffer k per unit of time. We use the processes $\xi =\{\xi (t),t\ge 0\}$ as a proxy for the cumulative cost under a given admission control policy $\mathrm{\Delta }(·),$ where

This is an approximation of the realized cost because the first term on the right-hand side replaces the admission control penalties actually incurred with their means.

In order to derive the approximating Brownian control problem, we consider a sequence of systems indexed by a system parameter $n=1,2,\dots ;$ we attach a superscript of n to various quantities of interest. Following the approach used by Ata (2006), we assume that the sequence of systems satisfies the following heavy traffic assumption.

Assumption 1.

For $n\ge 1,$ we have that

$${\lambda }^n=n\lambda ,{\eta }_k^n={\eta }_k\sqrt{n} \mathit{\text{and}} {\mu }_k^n=n{\mu }_k+\sqrt{n}{\beta }_k,\hspace{0.17em}\hspace{0.17em}k=0,1,\dots ,K,$$

where $\lambda ,$${\mu }_k,{\eta }_k$, and ${\beta }_k$ are nonnegative constants. Moreover, we assume that

$$\lambda ={\mu }_0=\frac{{\mu }_k}{p_k}\hspace{0.17em} \mathit{\text{for}} \hspace{0.17em}k=1,\dots ,K.$$

One starts the approximation procedure by defining suitably centered and scaled processes. For $n\ge 1,$ we define

In what follows, we assume

where $I_k(·)$ is the limiting idleness process for server k; see Harrison (1988) for an intuitive justification of (73).

Then, defining

and using Equations (67), (68), and (73), it is straightforward to derive the following for $t\ge 0$ and $k=1,\dots ,K:$

Moreover, it follows from Equation (71) that ${\mathrm{\Delta }}_k(t)$ is absolutely continuous. We denote its density by ${\delta }_k(·);\text{that} \text{is}.$

where ${\delta }_k(t)\in [0,1].$ Using this, we write

Then, passing to the limit formally as $n\to \infty$ and denoting the weak limit of $({\widehat{Q}}^n,{\widehat{X}}^n,{\widehat{\xi }}^n)$ by $(Z,\chi ,\xi ),$ where $\chi$ is a $(K+1)$-dimensional driftless Brownian motion with covariance matrix (see Appendix C for its derivation)

we deduce from (74), (75), and (76) that

In order to streamline the notation, we make the following change of variables:

and we let

Lastly, we define the set of negative drift vectors available to the system manager as in Equation (64). As a result, we arrive at the following Brownian system model:

which can be written as in Equation (2) with $d=K+1,$ where the reflection matrix R is given by Equation (63). Moreover, the processes Y, Z inherit properties in Equations (67)–(77) from their prelimit counterparts in the queueing model (cf. Equations (69) and (70)).

To minimize technical complexity, we restrict attention to the stationary Markov control policies as done in Section 3. That is, $\theta (t)=u(Z(t))$ for $t\ge 0$ for some policy function $u:{\mathbb{R}}_{+}^d\to \mathrm{\Theta }.$ Then, defining $c={(c_0,c_1,\dots ,c_K)}^{\top },\hspace{0.17em}h={(h_0,h_1,\dots ,h_K)}^{\top }$, and

as in Equation (65), the cumulative cost incurred over the time interval [0, t] under policy u can be written as in Equation (11). Note that $C^u(t)$ and $\xi (t)$ differ only by a term that is independent of the control. Given $C^u(t),$ one can formulate the discounted control problem as done in Section 3.1. Similarly, the ergodic control problem can be formulated as done in Section 3.2.

6.2.1. Interpreting the Solution of the Drift Control Problem in the Context of the Queueing Network Formulation.

Because the instantaneous cost rate $c(z,\theta )$ is linear in the control, inspection of the HJB equation reveals that the optimal control is of a bang-bang nature. That is, ${\theta }_k(t)\in \{{\underset{¯}{\theta }}_k,{\overline{\theta }}_k\}$ for all k, t as stated in Equation (66). This can be interpreted in the context of the queueing network displayed in Figure 1 as follows. For $k=0,1,\dots ,K$, whenever ${\theta }_k(t)={\overline{\theta }}_k$, the system manager turns away the thin-stream jobs arriving to buffer k externally (i.e., she shuts off the renewal process $A_k(·)$ at time t). Otherwise, she admits them to the system. Of course, the optimal policy is determined by the gradient $\nabla V(z)$ of the value function through the HJB equation, which we solve for using the method described in Section 5.

6.3. Related Example with Quadratic Cost of Control

Çelik and Maglaras (2008) and Ata and Barjesteh (2023) advance formulations where a system manager controls the arrival rate of customers to a queueing system by exercising dynamic pricing. One can follow a similar approach for the feed-forward queueing networks displayed in Figure 1 with suitable modifications (e.g., the dashed arrows also correspond to the arrivals of regular jobs). This ultimately results in a problem of drift control for RBM with the cost of control

where $\underset{¯}{\theta }$ is the drift rate vector corresponding to a nominal price vector.

6.4. Two Parametric Families of Benchmark Policies

Recall that the optimal policy can be characterized as

6.4.1. The Benchmark Policy for the Main Test Problem.

In our main test problem (see Section 6.1), we have $c(z,\theta )=h^{\top }z+c^{\top }\theta$. Therefore, it follows from (80) that for $k=0,1,\dots ,K$,

Namely, the optimal policy is of a bang-bang type. Therefore, we consider the following linear boundary policies as our benchmark polices. For $k=0,1,\dots ,K$,

where ${\beta }_0,{\beta }_1\dots ,{\beta }_K\in {\mathbb{R}}^{K+1}$ are vectors of policy parameters to be tuned.

In our numerical study, we primarily focus attention on the symmetric case, where

The symmetry allows us to limit the number of parameters needed for the benchmark policy. To be more specific, because of this symmetry, the downstream buffers look identical. As such, we restrict attention to parameter vectors of the following form:

The parameter vector ${\beta }_0$, which is used to determine the benchmark policy for buffer 0, has two distinct parameters: ${\varphi }_1$ and ${\varphi }_2$. In considering the policy for buffer 0, ${\varphi }_1$ captures the effect of its own queue length, whereas ${\varphi }_2$ captures the effects of the downstream buffers $1,\dots ,K$. We use a common parameter for the downstream buffers because they look identical from the perspective of buffer 0. Similarly, the parameter vector ${\beta }_i\hspace{0.17em}(i=1,\dots ,K)$ has three distinct parameters: ${\varphi }_3,\hspace{0.17em}{\varphi }_4$, and ${\varphi }_5$, where ${\varphi }_3$ is used as the multiplier for buffer 0 (the upstream buffer), ${\varphi }_5$ is used to capture the effect of buffer i itself, and ${\varphi }_4$ is used for all other downstream buffers. Note that all ${\beta }_i$ use the same three parameters ${\varphi }_3,\hspace{0.17em}{\varphi }_4$, and ${\varphi }_5$ for $i=1,\dots ,K$. They only differ with respect to the position of ${\varphi }_5$ (i.e., it is in the $i+1^{\text{st}}$ position for ${\beta }_i$).

In summary, the benchmark policy uses five distinct parameters in the symmetric case. This allows us to do a brute-force search via simulation on a five-dimensional grid, regardless of the number of buffers.

6.4.2. The Benchmark Policy for the Test Problem with the Quadratic Cost of Control.

In this case, substituting Equation (79) into Equation (80) gives the following characterization of the optimal policy:

Namely, the optimal policy is affine in the gradient. Therefore, we consider the following affine rate policies as our benchmark polices. For $k=0,1,\dots ,K$,

where ${\beta }_0,{\beta }_1\dots ,{\beta }_K\in {\mathbb{R}}^{K+1}$ are vectors of policy parameters to be tuned. We truncate this at the upper bound ${\overline{\theta }}_k$ if needed.

We focus attention on the symmetric case for this problem formulation too. To be specific, we assume

Because of this symmetry, the downstream buffers look identical. As such, we restrict attention to parameter vectors of the following form:

As done for the first benchmark policy above, this particular form of the parameter vectors can be justified using the symmetry as well.

6.5. Parallel-Server Test Problems

In this section, we consider a problem whose solution can be derived analytically by considering a one-dimensional problem. To be specific, we consider the parallel-server network that consists of K identical single-server queues as displayed in Figure 2. Clearly, this network can be decomposed into K separate single-server queues, leading to K separate one-dimensional problem formulations, which can be solved analytically; see Appendix D for details. For this example, we have that $R=I_{d\times d}$ and $A=I_{d\times d}$. In addition, we assume that the action space $\mathrm{\Theta }$ and the cost function $c(z,\theta )$ are the same as above.

7. Computational Results

For the test problems introduced in Section 6, we now compare the performance of policies derived using our method (see Section 5) with the best benchmark we could find. The results show that our method performs well, and it remains computationally feasible up to at least dimension $d=30$. We implement our method using three-layer or four-layer neural networks with the “elu” activation function (Rasamoelina et al. 2020) in Tensorflow 2 (Abadi et al. 2016) and using code adapted from that of Han et al. (2018) and Zhou et al. (2021b); see Appendix E for further details of our implementation.¹ In our numerical experiments, we observed that the performance of our method is robust to most hyperparameters. For tuning our algorithm, the activation function proved to be the most critical element. The “elu” activation function resulted in superior performance (Rasamoelina et al. 2020). Additionally, we found that decaying the learning rate to 0.0001 helped achieve good performance.

For our main test problem with the linear cost of control (introduced previously in Section 6.1) and also, for its variant with the quadratic cost of control (Section 6.3), the following parameter values are assumed: $h_0=2$, $h_k=1.9$ for $k=1,\dots ,K$; $c_k=1$ for $k=0,\dots ,K$; and $p_k=1/K$ for $k=1,\dots ,K$. Also, the reflection matrix R and the covariance matrix A for those families of problems are as follows:

We also consider a variation of our main test problem that has asymmetric routing probabilities.

As stated previously in Section 6.5, the reflection matrix and the covariance matrix for our parallel-server test problems are $R=I_{d\times d}$ and $A=I_{d\times d}$. Problems in that third class have $K=d$ buffers indexed by $k=1,\dots ,K$, and we set $h_1=2$ and $h_k=1.9$ for $k=2,\dots ,K$.

7.1. Main Test Problem with Linear Cost of Control

For our main test problem with linear cost of control (Section 6.1), we take

Also, interest rates $r=0.1$ and $r=0.01$ will be considered in the discounted case.

To begin, let us consider the simple case where $K=0$ (that is, there are no downstream buffers in the queueing network interpretation of the problem), and hence, $d=1$. In this case, one can solve the HJB equation analytically; see Appendix D for details. For the discounted formulation with $r=0.1$, Figure 3 compares the derivative of the value function computed using the analytical solution, which is shown in blue, with the neural network approximation for it that we computed using our method, which shown in red. (The comparisons for $r=0.01$ and the ergodic control case are similar.)

Combining Figure 3 with Equation (80), one sees that the policy derived using our method is close to the optimal policy. Table 1 reports the simulated performance with standard errors of these two policies based on 4 million sample paths and using the same discretization of time as in our computational method. Specifically, we report the long-run average cost under each policy in the ergodic control case, and we report the simulated value $V(0)$ in the discounted case. To repeat, the benchmark policy in this case is the optimal policy determined analytically but not accounting for the discretization of the timescale. Of course, all of the performance figures reported in Table 1 are subject to simulation errors. Finally, it is worth noting that our method took less than one hour to compute its policy recommendations using a 10-CPU core computer.

Table 1. Performance Comparison of Our Proposed Policy with the Benchmark Policy in the One-Dimensional Case ($K=0$)

Let us consider now the two-dimensional case ($K=1$), where the optimal policy is unknown. Therefore, we compare our method with the best benchmark we could find: the linear boundary policy described in Section 6.4. In the two-dimensional case, the linear boundary policy reduces to the following:

Through simulation, we perform a brute-force search to identify the best values of ${\beta }_0$ and ${\beta }_1$. The policies for $b=2$ and $b=10$ are shown in Figures 4 and 5, respectively, for the discounted control case with $r=0.1$. Our proposed policy sets the drift to b in the red regions and to zero in the blue regions in Figures 4 and 5, whereas the best linear boundary policy is represented by the white dashed lines in Figures 4 and 5. That is, the benchmark policy sets the drift to b in the region above and to the right of the dashed lines in Figures 4 and 5 and sets it to zero below and to the left of the dashed lines in Figures 4 and 5. Table 2 presents the costs with standard errors of the benchmark policy and our proposed policy obtained in a simulation study. The two policies have similar performance. Our method takes about one hour to compute policy recommendations using a 10-CPU core computer.

Table 2. Performance Comparison of Our Proposed Policy with the Benchmark Policy in the Two-Dimensional Case ($K=1$)

We then consider the six-dimensional case ($K=5$), where the linear boundary policy reduces to

Although there appears to be 36 parameters to be tuned, recall that we reduced the number of parameters to 5 in Section 6.4 by exploiting symmetry. This makes the brute-force search computationally feasible. Table 3 compares the performance with standard errors of our proposed policies with the benchmark policies. They have similar performance. In this case, the running time for our method is several hours using a 10-CPU computer.

Table 3. Performance Comparison of Our Proposed Policy with the Benchmark Policy in the Six-Dimensional Case $d=6$ ($K=5$)

To illustrate the scalability of our approach, we next consider the 21-dimensional case ($K=20$), where $p_1=\cdots =p_{20}=0.05$. Because of symmetry, we can perform a brute-force search for the best linear boundary policy as done earlier. The results, presented in Table 4, demonstrate that our method’s performance is comparable with that of the best benchmark. The run time for our method is about one day in this case using a 20-CPU core computer.

Table 4. Performance Comparison of Our Proposed Policy with the Benchmark Policy in the 21-Dimensional Case $d=21$ ($K=20$)

To further demonstrate the effectiveness of our approach, we consider a six-dimensional test problem with asymmetric routing probabilities. Specifically, for the example shown in Figure 1, we set

All other problem parameters remain the same as in the earlier six-dimensional symmetric test problem; see Appendix F for its reflection matrix R and the covariance matrix A. However, for the asymmetric problem, tuning its 36 parameters for the linear boundary policy becomes computationally prohibitive. Therefore, we propose an alternative approach in Appendix F to identify an effective boundary policy. Table 5 presents the performance, along with standard errors, of our proposed policies compared with benchmark policies. The two policies exhibit similar performance, and the run time of our method is comparable with that of the earlier symmetric six-dimensional example.

Table 5. Performance Comparison of Our Proposed Policy with the Benchmark Policy in the Six-Dimensional Case $d=6$ ($K=5$) with Asymmetric Routing Probabilities

7.2. Test Problems with Quadratic Cost of Control

In this section, we consider the test problem introduced in Section 6.3, for which we set ${\alpha }_k=1$ and ${\underset{¯}{\theta }}_k=1$ for all k. As in the previous treatment of our main test example, we report results for the cases of $d=1,2,6$ in Tables 6–8, respectively, where the benchmark policies are the affine rate policies discussed in Section 6.4, with policy parameters optimized via simulation through a brute-force search. We observe that our proposed policies outperform the best affine rate policies by very small margins in all cases.

Table 6. Performance Comparison of Our Proposed Policy with the Benchmark Policy in the Case of Quadratic Cost of Control and $d=1$

Table 7. Performance Comparison of Our Proposed Policy with the Benchmark Policy in the Case of Quadratic Cost of Control and $d=2$$(K=1)$

Table 8. Performance Comparison of Our Proposed Policy with the Benchmark Policy in the Case of Quadratic Cost of Control and $d=6$$(K=5)$

In the one-dimensional ergodic control case ($K=0$), we obtain analytical solutions to the RBM control problem in closed form by solving the HJB equation directly, which reduces to a first-order ordinary differential equation in this case; see Appendix D for details. Figure 6 compares the derivative of the optimal value function (derived in closed form) with its approximation via neural networks in the ergodic case. Combining Figure 6 with Equation (81) reveals that our proposed policy is close to the optimal policy.

In the two-dimensional case, our proposed policy is shown in Figure 7 for the ergodic case, with contour lines showing the state vectors ($z_0,z_1$) for which the policy chooses successively higher drift rates. The white dashed lines in Figure 7 similarly show the states ($z_0,z_1$) for which our benchmark policy (that is, the best affine rate policy) chooses the drift rate ${\theta }_k=1.5$ (for $k=0$ in panel (a) of Figure 7 and $k=1$ in panel (b) of Figure 7).

7.3. Parallel-Server Test Problems

This section focuses on parallel-server test problems (see Section 6.5) to demonstrate our method’s scalability. As illustrated in Figure 2, the parallel-server networks are essentially K independent copies of the one-dimensional case. We present the results in Table 9 for $d=30$ and linear cost of control. When $b=2$, our policies perform almost equally as well as the optimal policy, whereas for $b=10$, our policies perform within 1% of the optimal policy. The run time for our method is about one day in this case using a 20-CPU core computer.

Table 9. Performance Comparison Between Our Proposed Policy and the Benchmark Policy for 30-Dimensional Parallel-Server Test Problems with Linear Cost of Control

For quadratic cost of control, we are able to solve the test problems up to at least 100 dimensions. The results for $d=100$ are given in Table 10, where the benchmark policies are the best affine rate policies (see Section 6.4). The performance of our policy is within 1% of the benchmark performance. The run time for our method is several days in this case using a 30-CPU core computer.

Table 10. Performance Comparison Between Our Proposed Policy and the Benchmark Policy for 100-Dimensional Parallel-Server Test Problems with Quadratic Cost of Control

8. Concluding Remarks

Consider the general drift control problem formulated in Section 3, assuming specifically that the instantaneous cost rate $c(z,\theta )$ is linear in $\theta$ and further assuming that the set of available drift vectors is a rectangle $\mathrm{\Theta }=[0,b_1]\times \cdots \times [0,b_d]$. If one relaxes such a problem by letting $b_i\uparrow \infty$ for one or more i, then one obtains what is called a singular control problem (cf. Kushner and Martins 1991). Optimal policies for such problems typically involve the imposition of endogenous reflecting barriers (that is, reflecting barriers imposed by the system controller in order to minimize cost) in addition to exogenous reflecting barriers that may be imposed to represent physical constraints in the motivating application.

There are many examples of queueing network control problems whose natural heavy traffic approximations involve singular control; see, for example, Martins and Kushner (1990), Krichagina and Taksar (1992), and Martins et al. (1996). In our follow-up paper (Ata et al. 2024), we extend the method developed in this paper for drift control in a natural way to treat singular control, and we illustrate that extension by means of queueing network applications.

Separately, the following are three desirable generalizations of the problem formulations propounded in Section 3 of this paper. Each of them is straightforward in principle, and we expect to see these extensions implemented in future work, perhaps in combination with mild additional restrictions on problem data. (a) Instead of requiring that the reflection matrix R have the Minkowski form (1), require only that R be a completely $\mathcal{S}$ matrix, which Taylor and Williams (1993) showed is a necessary and sufficient condition for an RBM to be well defined. (b) Allow a more general state space for the controlled process Z, such as the convex polyhedrons characterized by Dai and Williams (1996). (c) Remove the requirement that the action space $\mathrm{\Theta }$ be bounded.

Lastly, we have considered in this paper the PDEs that arise in performance analysis and optimal control of RBMs, assuming that those PDEs admit $C^2$ solutions. Borkar and Budhiraja (2005) have established the existence and uniqueness of viscosity solutions for such PDEs, extending the earlier work by Dupuis and Ishii (1991). To the best of our knowledge, there is no theory currently available concerning the existence and uniqueness of classical $C^2$ solutions, either exact or approximate, and we leave that exploration as a topic for future research.