\documentclass[11pt]{article}
% Packages
\usepackage{fullpage}
\usepackage{latexsym}
\usepackage{graphicx}
\usepackage{psfrag}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{subfigure}
\usepackage{enumerate}
\usepackage{setspace}

% Assaf's definitions
\input{mydefinitions.tex}

% My definitions
\newcommand{\Field}{\mathcal{F}}
\newcommand{\M}{\mathcal{M}}
\newcommand{\Real}{\mathbb{R}}
\newcommand{\RDPlus}{\Real^d_{+}}
\newcommand{\RD}{\Real^d}
\newcommand{\Prob}{\mathbb{P}}
\newcommand{\set}[1]{\left\{#1\right\}}
\newcommand{\norm}[1]{\|#1\|_1}
\newcommand{\crs}[1]{\left\langle#1\right\rangle}
\newcommand{\eqq}{\equiv}

% Page delimiters
\setlength{\oddsidemargin}{-.05in} \setlength{\evensidemargin}{.25in} \setlength{\textwidth}{6.5in}
\setlength{\topmargin}{-0.1in} \setlength{\textheight}{9in}

%\renewcommand{\baselinestretch}{1.3}  % increases spacing x1.3
\parindent=18pt
\parskip=6pt

\title{A Linear Programming Algorithm for Computing the Stationary Distribution of Semimartingale Reflected Brownian Motion}
\author{{\sf Denis  Saure}\thanks{Columbia Graduate School of Business, e-mail: {\tt
dsaure05@gsb.columbia.edu}}
\\Columbia University
\and {\sf Peter W. Glynn}\thanks{Management Science and Engineering, e-mail: {\tt glynn@stanford.edu}}
\\Stanford University
\and {\sf Assaf  Zeevi}\thanks{Columbia Graduate School of Business, e-mail: {\tt assaf@gsb.columbia.edu}. Research partially supported by NSF grant
DMI-0447562 and BSF grant \#2006075 }
\\Columbia University
}

\date{March  2008}
\begin{document}

\maketitle

\begin{abstract}
This paper proposes a linear programming algorithm for computing the stationary distribution of semimartingale reflected Brownian motion (SRBM),
which arises as an approximation of certain queueing networks operating in heavy-traffic. Our algorithm is based on a Basic Adjoint Relationship
(BAR) which characterizes the stationary distribution. Approximating the state space with a finite grid of points and using a finite set of ``test''
functions, the BAR reduces to a set of linear equations that can be solved using standard linear programming techniques. As the set of test functions
increases in complexity and the grid becomes finer, the sequence of stationary distribution estimates which arise as solutions to the algorithm
converges, in a suitable sense, to the stationary distribution of the SRBM. The algorithm is seen to produce good estimates of the stationary moments
as well as the entire distribution. Extensions to settings with parameter uncertainty, and to Ornstein-Uhlenbeck-type diffusions arising in the
many-server heavy-traffic regime, are discussed as well.
\end{abstract}

\begin{center}
\rule{5in}{2pt}\\
{\sc Preliminary Draft Copy -- Not for Distribution}\\
\rule{5in}{2pt}
\end{center}

\vspace*{1em}

\noi {\bf Short Title:} Computing the stationary distribution of SRBM\\
%\vspace*{1em}

\noi {\bf Keywords:} queueing, diffusion approximations, reflected Brownian motion,
 Halfin-Whitt regime, numerical methods, linear programming, stationary distribution.\\
%\vspace*{1em}

\noi {\bf 2008 AMS Subject Classification (Primary): } XXXX - XXXX - XXXX

\renewcommand{\baselinestretch}{1.3}  % increases spacing x1.3
\small \normalsize
%---------------------------------
%\newpage
%\tableofcontents
%\newpage

\section{Introduction}\label{sec:intro}

\textbf{Background and motivation.} Multiclass queueing networks are widely used as mathematical models of real-world complex systems, e.g.,
communication networks, manufacturing operations, and service systems. At the same time, relevant queueing network models are mostly intractable for
purposes of exact performance analysis. A possible remedy for this is to approximate the queueing network model with a more simplified structure that
supports the derivation of the pertinent performance measures. Perhaps the most prevalent approach in this context has focused on a class of
multidimensional diffusion processes, known as semimartingale reflected Brownian motions (SRBM) which serve as a proxy for the dynamics of the
underlying queueing network. Loosely speaking, a $d$-dimensional SRBM living on the nonnegative orthant $S := \Real^d_+$ behaves like an ordinary
Brownian motion while in the interior of $S$, and is confined to the orthant by ``instantaneous'' reflection at the boundary faces; a formal
definition will be given in $\S$\ref{sec:definitions}.

The connection between the queuing network and its ``Brownian'' counterpart was first established in Reiman's seminal paper \cite{Reiman:84}. There,
it is shown that for a certain family of open single-class queueing networks, the normalized vector of queue length processes converges in
distribution to a SRBM on the positive orthant, as the traffic intensity approaches one. The scope of queueing network models for which such weak
limit theory exists has been greatly expanded since; see, for example, the influential papers of Bramson \cite{Bramson:98} and Williams
\cite{Williams:98}.

Regardless of whether supporting limit theory exists, one can always informally replace the original network with a suitable SRBM and pursue
performance analysis of the latter as an approximation for the former. This is in essence the proposal advocated by Harrison and Nguyen
\cite{Harrison-Nguyen:90}, an approach that was dubbed QNET. In that paper, Harrison and Nguyen explain how to ``map'' the parameters of the
underlying queueing network to those of the approximating Brownian model (the SRBM), and then use the latter to derive estimates of steady-state
performance measures such as mean queue length, workload, etc. (It should be noted that even if the SRBM can be arrived at as a rigorous weak limit,
it is still unclear whether it is justified to approximate the stationary distribution of the underlying queueing network with that of the SRBM; for
cases when this interchange-of-limits has been established, see e.g., Gamarnik and Zeevi \cite{Gamarnik-Zeevi:06} and Gurvich and Zeevi
\cite{Gurvich-Zeevi}). Assuming that the approximating SRBM indeed admits a steady-state, the main task, in order to carry out the program in
Harrison and Nguyen \cite{Harrison-Nguyen:90}, is to compute the stationary distribution of the SRBM. This is the focal point of the present paper.

Consider an SRBM $Z= (Z(t):t\geq 0)$ with state-space $S:=\Real^d_+$, whose drift vector $\mu=(\mu_1,\hdots,\mu_d)$ and covariance matrix $\sigma\in
\Real^{d\times d}$ (assumed to be positive definite), dictate the dynamics of $Z$ in the interior of the orthant, and with matrix $R \in
\Real^{d\times d}$ which dictates the directions of reflection at the boundary. Suppose $Z$ admits a unique stationary distribution $\pi$. In most
cases, it is impossible to compute $\pi$ in closed form, and a plausible alternative is to compute it numerically. To this end, define

\begin{equation}\label{eq:operatorG}
 Gf := \frac{1}{2}\sum_{i=1}^d \sum_{j=1}^d \sigma_{ij} \frac{\partial^2f}{\partial x_i
\partial x_j} + \sum_{i=1}^d \mu_i \frac{\partial f}{\partial x_i},
\end{equation}
\begin{equation}\label{eq:operationD}
D_if(x) := R_i \cdot \nabla f(x) \mbox{  for }x \in F_i:= \{x\in S:x_i = 0\},
\end{equation}
where $R_i$ is the $i$-th column of the reflection matrix, and for two vectors $a,b\in\Real^d$, $a\cdot b$ denotes their scalar product, and
$\nabla(\cdot)$ denotes the gradient operator. Using Ito's lemma it is straightforward to see that the stationary distribution $\pi$ together with a
\emph{boundary measure} $\nu$ must satisfy
\begin{equation}\label{eq:bar-intro}
\int_S \left(Gf\right)(x)\, \pi(dx)  + \sum_{i=1}^d \int_{F_i} \left(D_if\right)(x)\, \nu_i(dx) = 0,
\end{equation}
for all $f\in C^2_b$ (the space of twice differentiable functions with bounded second derivative). The equation above is known as the {\it basic adjoint relationship} (BAR), and the aforementioned necessity was established in
Harrison and Williams \cite{Harrison-Williams:87a} (a more formal definition of the BAR will be advanced in $\S$\ref{sec:definitions}). The first
term on the left hand side of (\ref{eq:bar-intro}) is driven by the behavior of $Z$ in the interior of $S$, while the second term is related to the
behavior on the boundary faces $\set{F_i}$, and it is associated with the $d$-dimensional pushing process keeping $Z$ in the orthant. Dai an Kurtz
\cite{Dai-Kurtz:Unpublished}  established that this relationship is also sufficient to characterize the stationary distribution $(\pi)$ and boundary
measure $(\nu)$.

With this characterization at hand, Dai and Harrison \cite{Dai:90} developed an algorithm for computing the stationary distribution of SRBM, by
viewing the BAR as an orthogonality condition between an infinite dimensional space of test functions and the stationary distribution. By considering an
increasing sequence of finite dimensional approximations of $C^2_b$, they obtain a sequence of approximating densities of $\pi$ by means of
orthogonal projections. They prove that this sequence of densities converges in $L_2$ to the stationary density
of the SRBM, but their proof relies on a certain conjecture regarding densities that arise as solutions to the BAR. Later, Shen et al.
\cite{Shen-Chen:02} considered a variant of the Dai-Harrison algorithm for the case of an SRBM in a hypercube, using a finite element method (piecewise
polynomials) to form the finite dimensional approximations of the space $C^2$ (Shen and Chen \cite{Shen-Chen:03} extend
this algorithm to the positive orthant).

\textbf{Our approach.} We view (\ref{eq:bar-intro}) as characterizing $\pi$ via a solution to an infinite dimensional system of linear equations.
Obviously it is impractical to try to solve this directly, but one can try approximating $S$ and $C^2_b$ by sequences $\set{S_n:
n=1,\hdots}$ and $\set{\mathcal{F}_m : m=1,\hdots}$ of finite cardinality subsets of $S$ and finite-dimensional subspaces of $C^2_b$, respectively,
such that $S_n\uparrow S$ and $\mathcal{F}_m\uparrow C^2_b$ in a suitable sense. For fixed values of $n$ and $m$, (\ref{eq:bar-intro}) reduces to a
finite-dimensional linear program (LP), that can be solved efficiently. The hope is then that, roughly speaking, $\pi$ will emerge as a limit of a
sequence of optimal solutions to the aforementioned LP's.

To establish the convergence of the sequence of approximations to the stationary distribution $\pi$, we need to ensure that this sequence is
\emph{tight}, and this consideration leads us to impose additional constraints in the LP formulation. This, in turn, necessitates estimates of
interior and boundary moments. One of the contributions of this paper is to provide means for enforcing such constraints, and in doing so to
establish a certain relationship between the interior and boundary measure (Harrison and Williams \cite{Harrison-Williams:87a} conjectured that the
latter is a restriction of the former to the faces of the orthant). It is worth emphasizing that our proof of convergence does not require the
conjecture advanced in Dai and Harrison \cite{Dai-Harrison:92} and follow up papers. The LP formulation is also quite flexible and it is
straightforward to incorporate additional constraints beyond those that enforce tightness. We indicate how this can be used to get more stable
numerical results, and obtain accurate approximations to the \emph{entire distribution} as opposed to \emph{first moments} as reported in Dai an
Harrison \cite{Dai-Harrison:92} and most follow up papers. Finally, we also illustrate the flexibility of the LP approach by: (i) indicating how it
can be adapted to account for parameter uncertainty in the primitives; and (ii) explaining how it can be used to compute approximations to the
stationary distribution of diffusion processes arising in the many-server heavy-traffic regime.

\textbf{Related literature.} The linear programming approach pursued in the present paper originates with the work of Manne \cite{Manne:60} in the
context of discrete-time and finite state-space Markov chains. Hernandez and Lasserre \cite{Hernandez-Lasserre:98} extend this to analyze the
convergence of linear-programming approximations for discrete-time general state-space controlled Markov chains. In \cite{Hernandez-Lasserre:98}, the
authors approximate a discrete time analog of the BAR using a discrete probability distribution and a finite subspace of test functions. The main
focus there is not on the stationary distribution, but rather on minimizing a steady-state cost function (see also \cite{Hernandez-Lasserre} for
further discussion and references). Mendiondo and Stockbridge \cite{Mendiondo-stockbrige:98} extend this work to the continuous time setting in the
context of long-run average and discounted control problems, where the state and control spaces are assumed to be compact. The key features that
distinguish our work relative to these papers are: (i) feasibility in our sequence of approximating problems is guaranteed, while the work cited
above assumes feasibility in the respective measure-control space; (ii) we deal with particulars of the SRBM problem that involve a non-compact state
space, and reflection at the boundaries; and (iii) our focus is on directly computing the stationary distribution.

\textbf{The remainder of the paper.}  Section \ref{sec:definitions} provides a formal definition of the SRBM, a summary of known existence,
uniqueness and stability results, as well as the BAR condition. Section \ref{sec:algorithm} formulates our algorithm, and establishes its
convergence. Section \ref{sec:app-SRBM} presents applications to several illustrative SRBM instances, with extensive computational results. Section
\ref{sec:robust} discusses an extension of the algorithm to a setting where there is uncertainty with regard to certain parameters in the queueing
network, that constitutes the (heavy-traffic) pre-limit of the SRBM. Section \ref{sec:app-twoclassqueue} provides an application to an
Ornstein-Uhlenbeck-type diffusion that arises in the so-called many server heavy-traffic regime. Finally, \S\ref{sec:prob-meas-choice} discusses
generalizations of the algorithm and choice of approximating spaces $\set{S_n}$ and $\set{\mathcal{F}_m}$.

\section{Semi-martingale Reflected Brownian Motion in the Orthant} \label{sec:definitions}

This section reviews some basic proporties of SRBM in the orthant, and then derives a result which pertains to the boundary behavior of SRBM, that is
essential for the  purposes of the algorithm we develop in $\S$\ref{sec:algorithm}.

\subsection{Review and Basic Properties}
We start with a formal definition of SRBM, adapted from Williams \cite{Williams:98}.
\begin{definition}\label{def:SRBM} For an integer $d\geq 1$, let $S:=\RDPlus$ (the positive $d$-dimensional orthant). Let $\mu$ be a constant vector
in $\RD$, $\sigma$ a $d\times d$ non-degenerate covariance matrix (symmetric and strictly positive definite), and $R$ a $d \times d$ matrix. For each
$x\in S$, an SRBM  associated with the data $(S,\mu,\sigma,R,x)$ is a $\Field_t$-adapted, $d$-dimensional process $Z = (Z(t):t\geq 0)$ defined on
some filtered probability space $(\Omega,\Field,\Field_t,\Prob)$ such that:
\begin{enumerate}[(i)]
\item  $Z = X + RY$, $\Prob_x$-a.s.,
\item $\Prob_x$-a.s., $Z$ has continuous paths and $Z(t)\in S$ for all $t\geq 0$,
\item $X$ is a $d$-dimensional Brownian motion with drift vector $\mu$, covariance matrix $\sigma$
and $X(0)=x$. In addition $X(t)-\mu t$ is a $\Field_t$-martingale,
\item $Y$ is an $\Field_t$-adapted $d$-dimensional process such that under $\Prob$ it satisfies for
each $j=1,\hdots,d:$
\begin{enumerate}[a.)]
\item $Y(0)=0$,
\item $(Y_i(t):t\geq 0)$ is continuous and non-decreasing,
\item $Y_i(t)$ can increase only when $Z$ hits the face $F_i = \{x\in S:x_i=0\}$.
\end{enumerate}
\end{enumerate}
\end{definition}

Loosely speaking, an SRBM $Z$ behaves like the Brownian motion $X$ in the interior of $S$, and is confined to the orthant by instantaneous
``reflection'' at the boundary faces: when $Z$ hits $F_i$, the process $Y_i$ increases, pushing $Z$ in the direction $R_i$ and keeping it in $S$. The
most general condition currently known to ensure existence and uniqueness (in law) of SRBM in the orthant is that the matrix $R$ be \emph{completely
$\mathcal{S}$}.

\begin{definition}
A $d\times d$ matrix R is said to be $\mathcal{S}$ if there exits a $d$-dimensional vector $u\geq 0$ such that $R u>0$, and to be a completely
$\mathcal{S}$ if each of its principal submatrices is an $\mathcal{S}$ matrix.
\end{definition}
This condition is in fact both necessary (Reiman and Williams \cite{Reiman-Williams:88}) and sufficient (Taylor and Williams
\cite{Taylor-Williams:93}) for existence of SRBM.

With regard to the stationary distribution of SRBM, which constitutes the main focus of the present paper, Dupuis and Williams
\cite{Dupuis-Williams:94} proved that a sufficient condition for existence is that all solutions of an associated deterministic Skorohod problem are
attracted to the origin in finite time. In addition, any stationary distribution $\pi$ of the SRBM $Z$ is unique and equivalent to Lebesgue measure
on $S$. To complete the picture, for each $i \in \set{1,\hdots,d}$ there exists a finite Borel measure $\nu_i$ on $F_i$ such that $\nu_i$ is
absolutely continuous with respect to the ($d-1$)-dimensional Lebesgue measure on $F_i$, and for each bounded Borel function $f$ on $F_i$ and $t\geq
0$,
\begin{equation}\label{prop:boundary}
\mathbb{E}_{\pi}\left[\int_0^t f(Z(s))dY_i(s) \right] = t\int_{F_i}f(x) \nu_i(dx) ,
\end{equation}
where $\mathbb{E}_{\pi}$ denotes expectation with respect to the stationary distribution $\pi$ (i.e., when $Z(0)$ is distributed according to $\pi$).
As alluded to earlier, the algorithm proposed in this paper is based on the BAR condition (restated here):
\begin{equation}\label{eq:bar}
\int_S \left(Gf\right)(x)\, \pi(dx)  + \sum_{i=1}^d \int_{F_i} \left(D_if\right)(x)\, \nu_i(dx) = 0 \quad \mbox{for all } f\in C^2_b\, ,
\end{equation}
where the second order differential operator $G$ and the directional derivative $D_i$ are defined in (\ref{eq:operatorG}) and (\ref{eq:operationD}),
respectively. Necessity of the BAR was first shown by Harrison and Williams \cite{Harrison-Williams:87a} when the matrix $R$ is Minkowski ($I-R\geq
0$ and $I-R$ transient), and then by Dai and Harrison \cite{Dai-Harrison:92} for the completely $\mathcal{S}$ case. Sufficiency was later proven by
Dai and Kurtz \cite{Dai-Kurtz:Unpublished}.

\subsection{Relationship of the Interior and Boundary Measures}
As mentioned in $\S$\ref{sec:intro}, our plan is to approximate the BAR (\ref{eq:bar}) by a finite system of linear equations with a finite number of
unknowns, use this to construct a sequence of probability measures converging to $\pi$. Since $S$ is not compact, we will need to ensure the
tightness of such a sequence, which can be done by bounding expectations for both interior and boundary measures.

Using an explicit Lyapunov function argument, Glynn and Zeevi \cite{Glynn-Zeevi:07} prove the existence of exponential moments of the interior
measure for the case of $R$ being symmetric and positive definite. Budhiraja and Lee \cite{Budhiraja-Lee:06} provide a more general result which
hinges on the existence of a suitable Lyapunov function for the completely $\mathcal{S}$ case. Using this property of the interior measure, we can
show the following.

\begin{theorem}\label{thm:finite-boundary}
Assume R is completely $\mathcal{S}$, and suppose $\pi$ is the stationary distribution for $Z$ associated with boundary measures $\nu_i\,,\,
i=1,\hdots,d$. Then, there exits a $d$-dimensional vector $u$ with all components being strictly positive such that
\begin{equation*}
\int_{F_i} e^{u\cdot x}\, \nu_i(dx) <\infty \quad ,\mbox{ for } i=1,\hdots,d.
\end{equation*}
\end{theorem}
\begin{proof}
Consider a $d$-dimensional vector $v>0$ such that $R^{\top}v>0$ (we know such a vector exists, since the completely $\mathcal{S}$ class is closed
under transposition; see Reiman and Williams \cite{Reiman-Williams:88}). Fix $a>0$ and let $f(x)=\exp(a x\cdot v)$. Applying It\^{o}'s formula, we
have that
$$f(Z(t))-f(Z(0)) -\int_0^t (Gf)(Z(s))ds - \sum_{i=1}^d \int_0^t (D_i f)(Z(s)) dY_i(s)$$
is a local $P_x$-martingale with respect to a sequence of localizing stopping times $T_n := \inf\set{t\geq 0: \norm{Z(t)}\geq n}$ $n=1,\hdots$, where
$\norm{x}$ denotes the sum of absolutes values of the vector $x$. It follows that
$$\mathbb{E}_x[f(Z(t\wedge T_n))] - f(x) = \mathbb{E}_x[\int_0^{t\wedge T_n} (Gf)(Z(s))ds] + \mathbb{E}_x[\sum_{i=1}^d \int_0^{t\wedge T_n} (D_i f)(Z(s)) dY_i(s)],$$
where $\mathbb{E}_x[\cdot] = \mathbb{E}[\cdot|Z(0) =x]$ and $a\wedge b:=\min\set{a,b}$. For this choice of function $f$, $(Gf)(x) = f(x)\left[a (\mu
\cdot v) + \frac{a^2}{2}\sum_{i=1}^d \sum_{j=1}^d \sigma_{ij} v_i v_j \right]$ and $(D_i f)(x) = a f(x) [R_{i}\cdot v]$. Note that
\begin{equation}\label{eq:muv}
\mu\cdot v = \mu (R^{\top})^{-1}R^{\top}v < 0,
\end{equation}
since $R^{-1}\mu<0$ must hold for $\pi$ to exist. Define $\alpha := - a (\mu\cdot v)>0$  and $\beta_i:= a R_i\cdot v>0\,,\, i=1,\hdots,d$. The next
inequality follows from the sign of $\alpha$ and the nonnegativity of $f$,
\begin{eqnarray*}
%\mathbb{E}_x[f(Z(t\wedge T_n))] - f(x) - \mathbb{E}_x[\int_0^{t\wedge T_n} (Gf)(Z(s))ds] =  \mathbb{E}_x[\sum_{i=1}^d \beta_i \int_0^{t\wedge T_n}
%f(Z(s)) dY_i(s)]\\
\mathbb{E}_x[f(Z(t\wedge T_n))] - f(x) + \alpha \mathbb{E}_x[\int_0^{t} f(Z(s))ds] \geq  \sum_{i=1}^d \beta_i \mathbb{E}_x[\int_0^{t\wedge T_n}
f(Z(s)) dY_i(s)].\\
\end{eqnarray*}
 Now note that
\begin{eqnarray*}
\mathbb{E}_x[f(Z(t\wedge T_n))] & = & \mathbb{E}_x[f(Z(t\wedge T_n))\,;\, T_n > t] + \mathbb{E}_x[f(Z(t\wedge T_n))\,;\, T_n \leq t]\\
 & \leq & \mathbb{E}_x[f(Z(t))\,;\, T_n > t] + \phi(n) \Prob[T_n \leq t]\\
 & \leq  & \mathbb{E}_x[f(Z(t))] + e^{cn} \Prob\left\{\max_{s\in[0\,,\,t]}\norm{Z(s)} \geq n\right\},\\
\end{eqnarray*}
where $\phi(n) = \max\set{f(x) : \norm{x} = n}\leq e^{c n}$ for some $c>0$. Using the El Kharroubi oscillation inequality (Lemma 1 in
\cite{Bernard-Kharroubi:91}) we have that
\begin{eqnarray*}
\max_{s\in[0\,,\,t]}\norm{Z(s)}  & \leq & K \max_{s\in[0\,,\,t]}\norm{x + \mu s + A B(s)} \\
                            & \leq & K(\norm{x} + c_1 t + c_2 \max_{s\in[0\,,\,t]}\norm{ B(s)}) \\
                            & \leq & K(\norm{x} + c_1 t + (c_2 d) \max_{s\in[0\,,\,t]}\max_{i=1,\hdots,d}|B_i(s)|), \\
\end{eqnarray*}
where $0< K<\infty$, $B(s)$ is a $d$-dimensional standard Brownian motion, $A$ is the Cholesky decomposition of $\sigma$ with $\sigma:=AA^{\top}$,
$c_1$ is a constant such that $\norm{\mu}\leq c_1$, and $c_2$ is such that $\sum_{ij}|A_{ij}| \leq c_2$. Using this, we have that
\begin{eqnarray*}
\Prob\set{\max_{s\in[0\,,\,t]}\norm{Z(s)} \geq n} & \leq & \Prob\set{\max_{i=1,\hdots,d} \max_{s\in[0\,,\,t]}|B_i(s)| \geq (n-K(\norm{x}+c_1t))/c_2}\\
                                         & \leq & d\Prob\set{\max_{s\in[0\,,\,t]}|B_1(s)| \geq (n-K(\norm{x}+c_1t))/K c_2}\\
                                         & \leq & 2d\Prob\set{\max_{s\in[0\,,\,t]}B_1(s) \geq (n-K(\norm{x}+c_1t))/K c_2}.\\
\end{eqnarray*}
A standard Gaussian tail bound gives $\Prob\set{\max_{s\in[0\,,\,t]}B_1(s) \geq x} \leq \sqrt{\frac{t}{2 \pi }}\frac{4}{x}e^{-x^2/2t}$, therefore we
can conclude that
\begin{equation*}
\Prob\set{\max_{s\in[0\,,\,t]}\norm{Z(s)} \geq n} \leq k_1 \sqrt{t}e^{-k_2 n^2/t}
\end{equation*}
for some $0<k_1,k_2<\infty$, and $n>K(\norm{x} + c_1t)$. Combining these results we have that for $n$ sufficiently large,
\begin{equation*}
\mathbb{E}_x\left[\sum_{i=1}^d \beta_i \int_0^{t\wedge T_n} f(Z(s)) dY_i(s)\right]\leq \mathbb{E}_x[f(Z(t))] + k_1\sqrt{t}e^{c n-k_2n^2/t} - f(x) +
\alpha \mathbb{E}_x\int_0^{t} f(Z(s))ds.
\end{equation*}
Letting $n\rightarrow \infty$, and invoking monotone convergence on the right-hand-side (since $f\geq0$) we have;
\begin{equation*}
\mathbb{E}_x\left[\sum_{i=1}^d \beta_i \int_0^{t} f(Z(s)) dY_i(s)\right]\leq  \mathbb{E}_{x}[f(Z(t)] -f(x)  +  \alpha \mathbb{E}_x\int_0^{t}
f(Z(s))ds.
\end{equation*}
Taking expectations with respect to $\pi$, and taking $a$ sufficiently small so that $\mathbb{E}_{\pi}[f(Z(0))]$ is finite, we have that
$$\sum_{i=1}^d \beta_i \mathbb{E}_{\pi}\left[ \int_0^{t} f(Z(s)) dY_i(s)\right] \leq t\alpha
\mathbb{E}_{\pi}[f(Z(0))].$$ Now, for $M>0$ consider $f_M(x) = f(x) \mathbf{1} \set{f(x)\leq M}$. Using (\ref{prop:boundary}) we have
\begin{equation*}
\mathbb{E}_{\pi}\left[\int_0^t f_M(Z(s))dY_i(s) \right] = t\int_{F_i}f_M(x)\nu_i(dx).
\end{equation*}
Leting $M\uparrow\infty$ and using monotone convergence, we conclude that
\begin{equation*}
\mathbb{E}_{\pi}\left[\int_0^t f(Z(s))dY_i(s) \right] = t\int_{F_i}f(x) \nu_i(dx), \quad i=1,\hdots,d.
\end{equation*}
It follows that
\begin{equation}\label{eq:connection}
\sum_{i=1}^d \beta_i \int_{F_i} f(x) \nu_i(dx) \leq \alpha \int_S f(x) \pi(dx).
\end{equation}
This concludes the proof. \qed
\end{proof}

%\begin{remark} {\rm A direct consequence of Theorem \ref{thm:finite-boundary} is the finiteness of the moment generating function of the stationary boundary
%measures in a neighborhood of zero. The nature of the result displayed in (\ref{eq:connection}) is somewhat reminiscent of a conjecture made by
%Harrison and Williams \cite{Harrison-Williams:87a} regarding the connection between the interior and boundary measures. In particular, based on their
%analysis of the skew-symmetric case, they conjecture that $\pi$ should have a density function $p$ on $S$ such that the boundary measure $\nu_i$ has
%density $\sigma_{ii}/(2R_{ii}) p$ on $F_i$.}
%\end{remark}

\section{The Proposed Algorithm}\label{sec:algorithm}

\subsection{LP Formulation}\label{sec:LPformulation}
The BAR condition (\ref{eq:bar}) which characterizes the stationary distribution $\pi$, can be visualized as an infinite dimensional system of linear
equations, and as such it is impractical to solve it directly. The approach we adopt here it to develop a sequence of approximations which involve:
(a) a finite system of linear equations; and (b) a suitable discretization of $\pi$. Let $\widetilde{S}$ be a dense countable subset of $S$, and
consider an increasing sequence of finite cardinality sets $S_n \subseteq \widetilde{S}$ such that $S_n\uparrow \widetilde{S}$, and $0\in S_n$ for
all $n\geq 1$. Also, let $\set{\mathcal{F}_m}$ be an increasing sequence of finite dimensional subspaces of $C^2_b$, such that for all $f\in C^2_b$
and for any compact set $A\subseteq S$, there exists a sequence $\set{f_m}$ with $f_m\in\mathcal{F}_m$ for each $m=1,\hdots$, such that
\begin{equation*}
\lim_{m\rightarrow\infty}\sup_A | (G(f_m-f))(x) |= 0,
\end{equation*}
\begin{equation*}
\lim_{m\rightarrow\infty}\sup_A | (D_j(f_m-f))(x) |= 0 \,,\,\mbox{ for all }j=1,\hdots,d,
\end{equation*}
where $G$ is the second order operator defined in (\ref{eq:operatorG}), and $D_i$ are the differential operators defined in (\ref{eq:operationD}).
Consider the following linear program.

$\hspace{-0.6cm}(\mathcal{P}_{n m})$

$\hspace{0.5cm }\min \quad\quad u $
\begin{eqnarray}\label{eq:Discrete-Bar}
 \mbox{s.t.} & & \left| \sum_{i: \,x_i\in S_n} (Gf)(x_i) \lambda_i + \sum_{j=1}^d
 \sum_{i: \,x_i\in F_j^n} (D_jf)(x_i) \gamma_{ji} \right|\leq u\,\, ,\,\, \mbox{for all } f\in \mathcal{F}_m \hspace{0.2cm} [\approx \mbox{BAR}]\\\label{eq:prob-dist}
 & & \sum_{i: \, x_i\in S_n} \lambda_i = 1 \hspace{7.9cm} [\mbox{Normalization}]\\\label{eq:tightness-dist}
 & & \sum_{i: \,x_i\in S_n} g(x_i)\lambda_i \leq M \hspace{4.9cm} [\mbox{Interior measure tightness}]\\\label{eq:tightness-boundary}
 & & \sum_{i: \,x_i\in F^n_j} g_i(x_i)\gamma_{ji} \leq M_j \,, \mbox{ for all } j\in\set{1\hdots,d} \hspace{0.45cm} [\mbox{Boundary measure tightness}]\\\label{eq:nonnegative}
 & & \lambda_i,\, \gamma_{ij} \geq 0 \,\,,\, \mbox{ for all }\set{i:\,x_i\in S_n} \mbox{ and } j\in\set{1,\hdots,d}.\hspace{1.25cm} [\mbox{Nonnegativity}]
 \end{eqnarray}
In this formulation, $\set{\lambda_i\,:\,i \mbox{ s.t. }x_i\in S_n}$ is a discrete distribution, with $S_n$ as its support, that approximates $\pi$,
and $\set{\gamma_{ji}\,:\,i \mbox{ s.t. }x_i\in F^n_j}$ is a discrete approximation to $\nu_j\, , \, j=1,\hdots,d$. The constants $M_j$,
$j=1,\hdots,d$, and $M$, and the functions $g_j(\cdot)$, $j=1\hdots,d$ and $g(\cdot)$ in (\ref{eq:tightness-dist})) and (\ref{eq:tightness-boundary})
are specified as follows. First, the function $g:S\rightarrow\Real_+$ is chosen such that $g(x)\leq g(y)$ for $x\leq y$ componentwise, and such that
$g(x)\rightarrow\infty$ as $\norm{x}\rightarrow\infty$, and with $M>0$ such that $\mathbb{E}_{\pi}(g(Z))\leq M$. The functions
$g_i:F_i\rightarrow\Real_+$ are chosen such that $g_j(x)\leq g_j(y)$ for $x\leq y$ componentwise, and such that $g_i(x)\rightarrow\infty$ as
$\norm{x}\rightarrow\infty$, with $M_i>0$ such that $\int_{F_i}g_i(x)\nu_i(dx)\leq M_i$. The constraint (\ref{eq:Discrete-Bar}) is a discrete
approximation for the BAR, while constraints (\ref{eq:prob-dist}) and (\ref{eq:nonnegative}) impose the required regularity on $\set{\lambda_i}$ and
$\set{\gamma_{ji}}$, $j=1,\hdots,d$. With a forward view towards the convergence result in $\S$\ref{sec:convergence}, constraints
(\ref{eq:tightness-dist}) and (\ref{eq:tightness-boundary}) are put in place to impose tightness on the sequence of solutions to these LP's. The
actual specification of the functions and constants in (\ref{eq:tightness-dist}) and (\ref{eq:tightness-boundary}) will be discussed in
$\S$\ref{sec:issues}.

\begin{remark}\label{re:convex-distr-interpretation}
{\rm Our LP formulation constructs an approximation of $\pi$ using discrete distributions with support in $S_n$. That is, this approximation is a
convex combination of point mass distributions, where $\lambda_i$ is the weight given to the distribution function that assigns probability one to
the point $x_i\,\in\, S_n$. It is not difficult to generalize this construction so that for any given $n$, $\pi$ $(\nu)$ will be approximated using a
convex combination of primitive distribution functions (finite measures). For this purpose, $\{S_n\}$ can be viewed as a sequence of function spaces
such that for any probability measure $\mu$ and finite measures $\nu_i$ $i=1,\hdots,d$, there exist sequences $\{\mu_n\}$ and $\{\nu_{ni}\}$,
$i=1,\hdots,d$, with $(\mu_n,\nu_n)\in S_n$, for which $\mu_n\Rightarrow \mu$ and $\nu_n \Rightarrow \nu$. For further details and discussion see
\S\ref{sec:prob-meas-choice}.}
\end{remark}

\subsection{Convergence Result}\label{sec:convergence}

Let $(\lambda,\underline{\gamma})$ be a feasible solution to $(\mathcal{P}_{n m})$, where $\underline{\gamma}:= (\gamma_{1},\hdots,\gamma_{d})$.
Define
\begin{equation*}
\pi^{n}_{\lambda}(x) := \sum_{i:\,x_i\in S_n} \lambda_i \delta_{x_i}(x),
\end{equation*}
\begin{equation*}
\nu^{n}_{\underline{\gamma} j}(x) := \sum_{i:\,x_i\in F_j^n} \gamma_{ji} \delta_{x_i}(x)\,\,\quad j=1,\hdots,d,
\end{equation*}
where $\delta_y(x)=1$ if $x=y$ and $0$ otherwise. These discrete measures provide an approximation to the true underlying probability measure $\pi$
and boundary measures $\nu:=(\nu_1,\hdots,\nu_d)$. The following is the main result of this paper, and characterizes the properties of the sequence
of approximating distributions generated by the LP ($\mathcal{P}_{nm}$). We use the notation ``$\Rightarrow$'' to mean weak convergence.

\begin{theorem}\label{theo:convergence}
Consider an SRBM $Z$ as given in Definition \ref{def:SRBM}. Let $\pi$ be the stationary distribution for $Z$, with boundary measures $\nu_j\,,
\,j=1,\hdots,d$. Then, there exist sequences $\set{n_k}$ and $\set{m_k}$, $k=1,2,\hdots$, of positive integers that diverge to infinity as
$k\rightarrow\infty$, such that $\pi_{\lambda^{k}}^{n_km_k}\Rightarrow \pi$ and $\nu^{n_km_k}_{\underline{\gamma}^{k}j}\Rightarrow \nu_j$ as
$k\rightarrow \infty$, for $j=1,\hdots,d$, where $(\lambda^k, \underline{\gamma}^k)$ are optimal solutions to the LP ($\mathcal{P}_{n_km_k}$).
\end{theorem}

\begin{proof} First, consider sequences $\set{\phi^n}$, and $\{\varphi^n_1\},\hdots, \{\varphi^n_d\}$, of probability and finite
measures on $S$ and $F_1,\hdots,F_d$, respectively, such that
\begin{itemize}
\item[(a)] for every $n\geq 1$, $\phi^n$ ($\varphi^n_j$) has support on $S_n$ ($F_j^n$),
that is, $\phi^n$ ($\varphi^n_j$) is of the form $\phi^n(x) = \pi^{nm}_{\theta^n}(x)$ ($\varphi^n_j(x) = \nu_{\vartheta^n j}^{nm}(x)$) for some
vector $\theta^n\geq 0$ ($\vartheta^n\geq 0$) such that  $\sum_i \theta^n_i = 1$ ($\sum_i \vartheta^n_i < \infty$);
%$\phi^n(x)=\sum_{i:\,x_i\in S_n} \theta^n_i\delta_{x_i}(x)$, ($\phi^n_j=\sum_{i:\,x_i\in F_j^n}
%\vartheta^n_{ji}\delta_{x_i}(x)$), with $\theta^n_i\geq 0$ $(\vartheta^n_{ji}\geq 0)$, and $\sum_{i:\,x_i\in S_n} \theta^n_i =1$,
\item[(b)] the sequences $\set{\phi^n}$ and $\{\varphi^n_1\},\hdots, \{\varphi^n_d\}$ converge weakly to $\pi$ and $\nu_1,\hdots,\nu_d$,
respectively;
\item[(c)] $\int_S g(x) \phi^n(dx) \leq \int_S g(x) \pi(dx)$, and $\int_{F_j} g(x) \varphi_j^n(dx) \leq \int_{F_j} g_j(x) \nu_j(dx)$, for
$j=1,\hdots,d$, for all $n\geq 1$.
\end{itemize}
It is straightforward to ensure the first two properties above hold. We next explain how the third can be made to hold as well. Fix $n\geq 1$.
Without lost of generality, assume the elements of $S_n= \set{x_1,x_2,\hdots,x_{|S_n|}}$ are such that $x_1<_lx_2<_l\hdots<_lx_{|S_n|}$, where $<_l$
represents lexicographic order, and $|\cdot|$ denotes the cardinality of a set. For $x_i\,\in S_n$ define $A_{x_i} := \set{x\in S \,:\, x_i <_l
x}\cup \set{x_i}$. Put $\widehat{A}_i:=\set{x\in A_{x_i}\,:\,i=\max\set{k: x\in A_{x_k}}}$. Notice that $\set{\widehat{A}_i}$ defines a partition of
$S$. Set $\theta_i := \pi(\widehat{A}_i)$, i.e., $\int_{x\in\widehat{A}_i}\pi(dx)$. Put $\theta^n:= (\theta_1,\hdots,\theta_{|S_n|})$, and $\phi^n :=
\pi^{nm}_{\theta^n}$.
 Consider a bounded continuous function $f:S\rightarrow\Real$, and the simple function approximation $f_n(x) := f(x_i)$ for $x\in\widehat{A}_i$.
Since $S_n\uparrow S$ we have that $f_n\rightarrow f$ by continuity of $f$. Hence,
$$ \int_{S} f(x)\phi_n(x) =  \int_{S} f_n(x)\pi(x) \rightarrow \int_{S} f(x)\pi(x) \quad\mbox{ as } n\rightarrow \infty,$$
by bounded convergence. Since this holds for all bounded continuous functions $f$, we conclude that $\phi_n\Rightarrow \pi$ as $n\rightarrow\infty$.
On the other hand, by construction $g(x_i)\leq g(x)$ for all $x\in \widehat{A}_i$, and $i=1,\hdots,|S_n|$. Therefore
$$ \int_S g(x)\phi_n(dx) \leq \int_S g(x) \pi(dx).$$ The same argument applies to each of the boundary measure approximations.

Notice that (c) implies that $(\theta^n, \vartheta^n_j)$ is ($\mathcal{P}_{nm}$)-feasible for all $m$ and $n$. Let $\mathcal{M}$ denote the space of
probability measures on $S$, and $\mathcal{M}_j$ the space of finite measures on $F_j \,,\, j=1\hdots,d$. For $\phi\in\mathcal{M}$,
$\varphi_j\in\mathcal{M}_j \,,\, j=1\hdots,d,$ and $f\in C^2_b$, let $B:\mathcal{M}\times\mathcal{M}_1\times\hdots\times\mathcal{M}_d\times C^2_b
\rightarrow \Real$ denote the following operator:
$$ B(\phi,\varphi,f) :=\int_S \left(Gf\right)(x)\, \phi(dx)  + \sum_{j=1}^d \int_{F_j} \left(D_jf\right)(x)\, \varphi_j(dx),$$
where $\varphi := (\varphi_1,\dots,\varphi_d)$. Fix $n$, and consider the sequence $\set{(\lambda^n_{m},\, \underline{\gamma}^n_m)}$ of optimal
solutions to ($\mathcal{P}_{nm}$), $m=1,2\hdots$. This sequence is tight, so it converges weakly, along a subsequence $\set{m^n_k}$, to a probability
distribution $\hat{\pi}^n$ and finite measures $\hat{\nu}^n := (\hat{\nu}^n_1,\dots,\hat{\nu}^n_d)$. Now, for any function $f \in C^2_b$ Consider the
sequence $\set{f_m}$ with $f_m\in \mathcal{F}_m$ such that $\lim_{m\rightarrow\infty} \max_{x\in S_n}|(G(f-f_m))(x)| = 0$, and
$\lim_{m\rightarrow\infty} \max_{x\in F_j^n}|(D_j(f-f_m))(x)| = 0$. Put $\pi^n_m:= \pi^n_{\lambda^n_{m}}$, and $\nu^n_m :=
\nu^n_{\underline{\gamma}^n_{m}}$. We have that
\begin{eqnarray*}
|B(\hat{\pi}^n,\, \hat{\nu}^n,\,f)| & = & |\lim_{k\rightarrow\infty}B(\pi^n_{m^n_k},\, \nu^n_{m^n_k},f)|\\
& = & |\lim_{k\rightarrow\infty}B(\pi^n_{m^n_k},\, \nu^n_{m^n_k},f-f_{m^n_k})+ \lim_{k\rightarrow\infty}B(\pi^n_{m^n_k},\, \nu^n_m,f_{m^n_k})|\\
& \leq & |\lim_{k\rightarrow\infty} B(\phi^n,\, \varphi^n,f_{m^n_k})| \\
& = & |B(\phi^n,\, \varphi^n,f)|.
\end{eqnarray*}
The sequence $\set{(\hat{\pi}^n_,\, \hat{\nu}^n)}$, $n = 1,2,\hdots$ is also tight, so it converges weakly, along a subsequence $\set{n_k}$, to a
probability distribution $\hat{\pi}$ and finite measures $\hat{\nu} := (\hat{\nu}_1,\dots,\hat{\nu}_d).$ From the above we have that
$$|B(\hat{\pi}, \hat{\nu},f)| = |\lim_{k\rightarrow\infty}B(\hat{\pi}^{n_k},\, \hat{\nu}^{n_k},f)| \leq |\lim_{k\rightarrow\infty} B(\phi^{n_k},\, \varphi^{n_k},f)|
= |B(\pi, \nu,f)| \,=\, 0.$$ Applying a diagonal argument, the
%Now, for each $n$ there exists a subsequence $\set{m_h}\subseteq\set{m_k}$ for which $\pi^{m_hn}_{\lambda^n_{m_h}}\Rightarrow \hat{\pi}^n$ and
%$\nu^{m_hn}_{\underline{\gamma}^n{m_h}}\Rightarrow \hat{\nu}^n_j$ for $j=1,\hdots,d$, with $\hat{\pi}^n$($\nu^n_j$, $j=1,\hdots,d$) proper
%probability(finite) measure(s), therefore exists an increasing sequence $\set{n_k}$ such that
%$\pi^{m_kn_k}_{\lambda^{n_k}_{m_k}}$($\nu^{m_kn_k}_{\gamma^{n_k}_{m_k j}}$) converges to $\hat{\pi}$($\hat{\nu}_j$, $j=1,\hdots,d)$.The \
result follows from the uniqueness of the stationary distribution, and hence $\hat{\pi}$ and $\hat{\nu}$ are solutions to the BAR.\qed
\end{proof}

\subsection{Practical Considerations}\label{sec:issues}
In order to implement the algorithm described $\S$\ref{sec:LPformulation}, we need to specify the sequence of subsets of $S$, $\set{S_n}$, the
sequence of finite dimensional spaces $\set{\mathcal{F}_m}$ approximating $C^2_b$, and the tightness conditions (for the latter, we need to specify
feasible values of $M$ and $M_j$ for $j=1\hdots,d$, as well as functions $g$ and $g_j\,, \,j=1,\hdots,d$).

\noindent\textbf{Tightness of the interior distribution.} For illustration purposes we will assume that $R$ is symmetric and that $R^{-1}$ is
positive definite; this facilitates the construction of a simple Lyapunov function (a more general, but less implicit argument can be followed using
a similar Lyapunov function to the one identified in Dupuis and Williams \cite{Dupuis-Williams:94}). Put $f(x) = x\cdot R^{-1}x$. Following the steps
in the proof of Theorem \ref{thm:finite-boundary}, applying Ito's formula and taking expectations with respect to $\Prob_x$, we have
$$\mathbb{E}_x[f(Z(t\wedge T_n))]-f(x)+  \mathbb{E}_x[\int_0^{t\wedge T_n} Z(s)\cdot \varrho \,ds] = C t + \mathbb{E}_x\left[\sum_{i=1}^d \int_0^t Z_i(s) dY_i(s)\right],$$
where $C = \sum_{i j} \sigma_{ij}(R^{-1})_{ij}>0$, $ \sigma_{ij}$ is the $(i,j)$-element of $\sigma$, $\varrho := -R^{-1}\mu>0$, and $T_n :=
\inf\set{t\geq 0\, :\, \norm{Z(t)}\geq n}$. Note that $\mathbb{E}_x[\sum_{i=1}^d \int_0^t Z_i(s) dY_i(s)] = 0$ by Definition \ref{def:SRBM}.
Considering the sign of $\varrho$, and that $f\geq0$, we have that
$$-f(x)+  \mathbb{E}_x[\int_0^{t\wedge T_n} Z(s)\cdot \varrho \,ds] \leq C t.$$
Applying the Monotone Convergence Theorem to the left-hand-side, and taking expectations with respect to $\pi$, we conclude that
$$\int_S (e\cdot x)\, \pi(dx) \leq \frac{C}{\min_j\set{\varrho_j}},$$
with $e = (1,\hdots,1)^{\top}\in\Real^d$. Here we can take $g(x) = e\cdot x$ and $M = C/\min_j\set{\varrho_j}$ (for the case where $R$ is Minkowski,
a bound can be derived directly from lemma 8.4 in Harrison and Williams \cite{Harrison-Williams:87a}).

\noindent\textbf{Tightness of the boundary measures.} Take a $d$-dimensional vector $v>0$ such that $R^{\top}v>0$ (which follows from the completely
$\mathcal{S}$ condition) and construct a $d \times d$ symmetric matrix $V$ by setting $V_{ij}=v_j v_i>0$. It follows that $(RV)_{ij}>0$ for
$i,j=1,\hdots,d$. Let $f(x) = x\cdot Vx$. Note that $(Gf)(x) = 2 \mu\cdot Vx + \sum_{ji} v_iv_j\sigma_{ij}$. Following the steps in the proof of
Theorem \ref{thm:finite-boundary}, for $n\geq1$, we have
\begin{eqnarray*}
\mathbb{E}_x[f(Z(t))] + k_1 \sqrt{t} n^2 e^{-k_2n^2/t} + \mathbb{E}_x[\int_0^{t} Z(s)\cdot \zeta ds]-f(x) &\geq & \mathbb{E}_x\left[\sum_{i=1}^d
\int_0^{t\wedge T_n} Z(s)\cdot (R_iV)\, dY_i(s)\right],
\end{eqnarray*}
where $\zeta := V\mu >0$, following (\ref{eq:muv}). Letting $n\rightarrow\infty$, taking expectations with respect to $\pi$, and using
(\ref{prop:boundary}) as in proof of Theorem \ref{thm:finite-boundary}, we have that
$$\int_{F_i}  (e\cdot x)\, \nu_i(dx) \leq \frac{\alpha}{\beta_i}\int_S  e\cdot x\, \pi(dx),$$
$i=1,\hdots,d,$ where $\alpha := \max_j\set{\zeta_j}>0$ and $\beta_i := \min_j\set{(V R_i)_j}>0,\, i=1,\hdots,d$. Combining this result with the
previous one, we can choose $g_i(x) =e\cdot x$, and $M_i = (C \alpha)/(\beta_i \min_j\set{\varrho_j})$.

\noindent\textbf{Distribution and function space approximation.} Harrison and Williams \cite{Harrison-Williams:87b} proved that the stationary
distribution for a standard SRBM\footnote{A SBRM is said to be standard if $\sigma_{ii} =1 \, i=1,\hdots,d$, $R$ is invertible and
$\|\varrho\|_{\infty}=1.$} has a separable density function (in the usual Cartesian coordinates) if and only if the covariance matrix $\sigma$
satisfies a ``skew symmetry'' condition:
\begin{equation}
2 \sigma_{ij} = R_{ij}+R_{ji} \quad \mbox{for} \quad i\neq j\,\in\set{1,\hdots,d}.\label{eq:skew-symmetric}
\end{equation}
In this case, the marginal distribution for each coordinate $i$ is exponential with rate $2\varrho_i:=(-2R^{-1}\mu)_i$. Moreover, the boundary
measures are the restriction of the joint distribution to the corresponding faces of $S$ (up to a scaling factor). Recently, Budhiraja and Lee
\cite{Budhiraja-Lee:06} established that the moment generating function of the stationary distribution is finite in a neighborhood of zero under the
sufficient stability condition of Dupuis and Williams \cite{Dupuis-Williams:94}.

The above observations motivate taking $\set{S_n}$ to have an ``exponential spacing'' on the marginals, that is for each $n\geq 1 $, the grid can be
chosen to be
\begin{equation}\label{eq:SRBM-grid}
S_n = \set{ x\in S \,:\, x_i \in\set{[\log(n)-\log(k)]/(2\varrho_i) \,, \mbox{ for } k=1,\hdots,n} \,, \mbox{ for } i=1,\hdots,d}.
\end{equation}
For the boundary, we simply project this grid on the corresponding faces, $ F^n_j = \set{ x\in S_n \,:\, x_j =0} \, , j=1,\hdots,d$ (this structure
of $S_n$ is related to the choice of a reference density in the Dai-Harrison algorithm \cite{Dai-Harrison:92}). Regarding the function spaces
approximating $C^2_b$, for each $m\geq 1$ let $\mathcal{H}_m$ be the linear space spanned by $\set{x_1^{p_1}\cdots x_d^{p_d} \,:\, p_1+\hdots+p_d
\leq m \,,\, p_1,\hdots,p_d \in \mathbb{N}}$. Take $\eps>0$ and define
$$\mathcal{H}^l_m:= \{f\in C^2_b \,:\, \exists \, g\in\mathcal{H}_m\mbox{ and } c>0 \mbox{, such that } f(x) = g(x)\mbox{, for }\norm{x}\leq l  \mbox{, and } f(x) = c \mbox{, for } \norm{x}>l\}.$$
For each $m\geq 1$ we take $\mathcal{F}_m = \bigcup_{1\leq l\leq m} \mathcal{H}^l_m$.

\begin{remark}{\rm A practical choice of the exponential spacing can be had using a rate vector $C\varrho$, with $0.5\leq C\leq2$. In the numerical experiments we take $\mathcal{F}_m := \mathcal{H}_m$.}
\end{remark}

\section{Numerical Results}\label{sec:app-SRBM}
In this section we compare the performance of our algorithm to some known analytical results derived for particular instances of SRBM, as well as
other numerical methods proposed in antecedent literature. All results reported in this section were computed using a regular Desktop PC with
processor Intel Pentium D (3.2Ghz), using a MATLAB implementation of the algorithm. The code and documentation can be downloaded from
http://www.columbia.edu/$\sim$drs2114/SRBM.

\textbf{Example 1: A two-dimensional special case.} Consider a two-dimensional SRBM $Z = (Z(t)\,:\,t\geq 0)$ associated to the following data:
$$ R = \left(\begin{array}{cc}1 & 0\\-1&1\\\end{array}\right) \quad, \quad \sigma = \left(\begin{array}{cc}1 & 0\\0&1\\\end{array}\right) \quad, \quad \mu = \left(\begin{array}{c}\mu_1 \\0\\\end{array}\right).$$
For this SRBM the stationary condition reduces to $\mu_1<0$. Harrison \cite{Harrison:78} computed a closed form solution for the density of the
stationary distribution in polar coordinates:
\begin{equation*}
p(x_1,x_2) = (2|\mu_1|)^{3/2}/(\pi^{1/2}) r^{-1/2}\exp(\mu_1r(1+\cos(\theta)))\cos(\theta/2),
\end{equation*}
where $(x_1,x_2) =(r\cos(\theta),r\sin(\theta))$. This SRBM is essentially the only instance, excluding the skew symmetry condition, for which the
stationary distribution can be computed in closed form.

Without loss of generality, set $\mu_1 =-1$. It can be shown (Greenberg \cite{Greenberg:86}) that
\begin{equation*}
\mathbb{E}_{\pi}Z_1 = 0.5 \quad ,\quad  \mathbb{E}_{\pi}Z_2 = 0.75.
\end{equation*}
Taking $S_n$ as in (\ref{eq:SRBM-grid}), with $n=100$, we use the proposed algorithm to get estimates for these moments. The results are displayed in
Table \ref{table:two-dim} for various degrees of polynomials that compose $\mathcal{F}_m$.
\begin{table}
  \centering
    \begin{tabular}{c||c|c|c|c|c|c|c|c}
        $m$ & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \\\hline \hline
        $\mathbb{E}_{\pi}Z_1$ & 0.5 & 0.5 & 0.5 & 0.5 & 0.5 & 0.5 & 0.5 & 0.5 \\\hline
        $\mathbb{E}_{\pi}Z_2$ & 0.793 & 0.821 & 0.764 & 0.768 & 0.762 &0.758 & 0.760 & 0.750 \\
    \end{tabular}

  \caption{\textbf{Moment estimates for the SRBM in Example 1}. The algorithm uses the grid $S_n$ in (\ref{eq:SRBM-grid}) with $n=100$ points on each coordinate, and $m=$ degree of polynomials. True
  moments are $\mathbb{E}_{\pi}Z_1=0.5$ and $\mathbb{E}_{\pi}Z_2=0.75$.}\label{table:two-dim}
\end{table}
We see that the estimates are quite accurate even for low values of $m$. Each one of the instances shown in Table \ref{table:two-dim} had a runtime
of less than $1$ minute. Figure \ref{figure:two-dim} gives the marginal distribution estimates for $n=50$ and $m=6$ (dotted line), compared with the
true marginal distributions (solid line).
\begin{figure}
\includegraphics[width=18cm]{N-figure1.eps}
  \caption{\textbf{Marginal distribution estimates for the SRBM in Example 1.} The algorithm uses the grid $S_n$ in (\ref{eq:SRBM-grid}) with $n=50$ grid points on each coordinate, and $m=6$ degree polynomials.
  The marginal distribution estimates (dotted line) are superimposed on the true marginals (solid line).}\label{figure:two-dim}
\end{figure}

\textbf{Example 2: A Symmetric case.} A standard SRBM is said to be \emph{symmetric} if its data has the following properties: $\sigma_{ij} = \rho$
for $1\leq i<j\leq d$, $\mu_i=-1$ for $i=1,\hdots,d$ and $R_{ij}=R_{ji}=-r\leq 0$ for $1\leq i<j\leq d$. That $\sigma$ is positive implies that
$-1/(d-1)<\rho<1$, and the completely $\mathcal{S}$ condition reduces to $r(d-1)<1$. This type of SRBM arises as a heavy-traffic limit of a symmetric
generalized Jackson network. Manipulating the BAR condition, Harrison and Dai \cite{Dai-Harrison:92} showed that
\begin{equation*}
\mathbb{E}Z_1=\hdots=\mathbb{E}Z_d = \frac{1-(d-2)r+(d-1)r\rho}{2(r+1)}.
\end{equation*}
We use our algorithm to compute estimates for these moments for the case $d=2$. Imitating the work of Harrison and Dai \cite{Dai-Harrison:92}, we let
$\rho$ take values in $\set{-0.9, -0.5,0.0,0.5,0.9}$ and $r$ in $\set{0.2, 0.4,0.6,0.8,0.9, 0.95}$. Table \ref{table:symmetric} shows the relative
errors between our estimates and the exact values for $m=6$ and $n=100$. The relative error is always less than $1\%$ and running times less than 1
minute.

\begin{table}
  \centering
    \begin{tabular}{c||c|c|c|c|c}
        $r\backslash\rho$ & -0.9 & -0.5 & 0.0 & 0.5 & 0.9  \\\hline \hline
        $0.2 $ & 1.98e-4 & 1.65e-4 & 1.26e-3 & 2.32e-3 & 6.66e-4  \\\hline
        $0.4 $ & 3.90e-4 & 4.88e-5 & 2.52e-5 & 2.92e-3 & 2.00e-4  \\\hline
        $0.6 $ & 1.76e-5 & 5.97e-4 & 2.05e-4 & 1.06e-6 & 5.70e-4  \\\hline
        $0.8 $ & 1.53e-5 & 4.27e-5 & 2.22e-5 & 1.17e-3 & 2.16e-3  \\\hline
        $0.9 $ & 1.65e-5 & 4.26e-5 & 5.43e-4 & 1.00e-3 & 3.24e-3  \\\hline
        $0.95$ & 4.57e-6 & 3.42e-5 & 2.46e-4 & 5.98e-4 & 4.82e-3  \\
    \end{tabular}

  \caption{\textbf{Relative errors for the SRBM moment estimates in Example 2.} Table depicts relative errors between moment estimates and true moments
  for various values of $r$ and $\rho$. Here the number of grid points on each coordinate is $n=100$ and the degree of polynomials is $m=6$.}\label{table:symmetric}
\end{table}


\textbf{Example 3: Skew-symmetric SRBM.} As indicated earlier, Harrison and Williams \cite{Harrison-Williams:87b} proved that a standard SRBM has a
product form stationary distribution if and only if $\varrho<0$ and condition (\ref{eq:skew-symmetric}) holds. In this, case we know that the $i$-th
marginal distribution is exponential with mean $1/(2\varrho_i)$, $i=1,\hdots,d$. Consider a two-dimensional SRBM associated with the following data:
$$ R = \left(\begin{array}{cc}1 & -0.6\\-0.25&1\\\end{array}\right) \quad ,\quad \sigma = \left(\begin{array}{cc}1 & -0.425\\-0.425&1\\\end{array}\right) \quad, \quad \mu = \left(\begin{array}{c}-0.85 \\0\\\end{array}\right).$$
One can check that condition (\ref{eq:skew-symmetric}) holds, and that $\varrho^{\top} =(-1,-0.25)$. This implies that $\mathbb{E}_{\pi}Z_1 =0.5$ and
$\mathbb{E}_{\pi}Z_2 =2$. Taking $n=100$ we use our algorithm to get estimates for these moments for various values of $m$, the mayor degree of the
polynomials, and the results are summarized in Table \ref{table:skew}.
\begin{table}
  \centering
    \begin{tabular}{c||c|c|c|c|c|c|c|c}
        $m$ & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \\\hline \hline
        $\mathbb{E}_{\pi}Z_1$ & 0.4942 & 0.4943 & 0.4968 & 0.5014 & 0.5005 & 0.5011 & 0.4996 & 0.5002 \\\hline
        $\mathbb{E}_{\pi}Z_2$ & 2.0095 & 2.0094 & 2.0053 & 1.9976 & 1.9990 &1.9981 & 2.0005 & 1.9996 \\
    \end{tabular}

  \caption{\textbf{Moment estimates for the SRBM in Example 3.} The algorithm uses $n=100$ grid points on each coordinate, and $m =$ degree of polynomials. True
  moments are $\mathbb{E}_{\pi}Z_1=0.5$ and $\mathbb{E}_{\pi}Z_2=2.0$.}\label{table:skew}
\end{table}
We see that our algorithm provides good estimates even for low degree polynomials. Each one of the instances shown in Table \ref{table:skew} had a
runtime of under $1$ minute. Figure \ref{figure:skew} depicts marginal distribution estimates for $n=150$ and $m=10$ (dotted line), compared with the
actual marginal distributions (solid line).
\begin{figure}
\includegraphics[width=18cm]{N-figure2.eps}
  \caption{\textbf{Marginal distribution estimates for the SRBM in Example 3.} The algorithm uses the grid $S_n$ as in (\ref{eq:SRBM-grid}) with $n=150$ grid points on each coordinate, and $m=10$ degree polynomials.
  The marginal distribution estimates (dotted line) are superimposed on the true marginals (solid line).}\label{figure:skew}
\end{figure}

\textbf{Example 4: Two queues in tandem and comparison against QNET and QNA methods.} This section follows closely the analysis presented in Chapter
4 of Dai \cite{Dai:90}. Consider a network of $d$ queues in tandem. Let $\rho_i$ denote the mean service time at station $i$, $C^2_{s_i}$ denote the
squared coefficient of variation of the service time distribution at station $i$, and $C^2_a$ denote the squared coefficient of variation for the
interarrival time distribution. Using the method proposed by Harrison and Nguyen \cite{Harrison-Nguyen:90}, Dai \cite{Dai:90} shows that the
$d$-dimensional current workload process can be approximated by a $d$-dimensional SBRM as in Definition \ref{def:SRBM}. When $d=2$ the parameters
$(\sigma,R,\mu)$ defining $Z$ are given by
$$ R = \left(\begin{array}{cc}1 & 0\\-\rho_2/\rho_1&1\\\end{array}\right) \quad, \quad \sigma = \left(\begin{array}{cc}\rho^2_1(C^2_a+C^2_{s_1}) & -\rho_1\rho_2C^2_{s_1}\\
-\rho_1\rho_2C^2_{s_1}& \rho^2_2(C^2_{s_1}+C^2_{s_2})\\\end{array}\right) \quad ,\quad \mu = \left(\begin{array}{c}\rho_1-1
\\\rho_2/\rho_1-1\\\end{array}\right).$$
The arrival rate is assumed to be 1, so that $\rho_1$ and  $\rho_2$ represents the traffic intensities at station 1 and 2, respectively.

Suresh and Whitt \cite{Suresh-Whitt:90} studied a system of two queues in tandem and considered various variability parameter triples
$(C^2_a,C^2_{s_1},C^2_{s_2})$, for all combinations of the traffic intensities $\rho_1$ and $\rho_2$ in a representative range. For $C^2_{s_1}\neq
C^2_{s_2}$ they considered five variability combinations: (0.5,0.5,2.0), (1.0,0.5,8.0), (1.0,2.0,4.0), (4.5,0.5,1.0) and , (4.0,1.0,4.0). We will
refer to these as \emph{Case 1} to \emph{Case 5}, respectively. For $C^2_{s_1}= C^2_{s_2}$ they consider two variability combinations: (1.0,0.5,0.5)
and (1.0,4.0,4.0). We will refer to these as \emph{Case 6} and \emph{Case 7}, respectively. When  $C^2_{s_1}\neq C^2_{s_2}$, they consider four
values of $\rho_i:$ 0.3, 0.6, 0.8, 0.9. When $C^2_{s_1}\neq C^2_{s_2}$, for each queue they consider five values of $\rho_i:$ 0.1, 0.2, 0.3, 0.6,
0.9, for each queue.

Suresh and Whitt \cite{Suresh-Whitt:90} conducted extensive simulation experiments in order to obtain estimates for the expected steady-state waiting
times. When $C_a^2 = 0.5$, the Erlang distribution was used with 2 degrees of freedom. When $C_a^2 = 1.0$, the exponential distribution was used.
When $C_a^2 > 1.0$, the hyperexponential distribution with balanced means was used. They also compare their simulation results with the results
obtained by using a software package described by Whitt \cite{Whitt:83}, dubbed the QNA method. Dai \cite{Dai:90} uses this set of experiments to
test the algorithm proposed in Harrison and Dai \cite{Dai-Harrison:92}, which they refer to as the QNET method. We use these results to benchmark the
performance of our algorithm against the QNA and QNET methods. For the purposes of our algorithm, we used $n=100$ and $m=6$ for all experiments. Each
scenario described above took less than 3 minutes of running time to yield the LP based estimates.

Tables \ref{Suresh-Whitt-Case1-Table} to \ref{Suresh-Whitt-Case7-Table} (see appendix \ref{sec:numerical}) give detailed simulation estimates, QNET
estimates, QNA estimates, and our LP-based estimates. Table \ref{Suresh-Whitt-Table} summarizes all balanced heavy traffic cases (equal traffic
intensities in each queue, close to 1), and gives an overall comparison of QNA estimates, QNET estimates and the LP-based estimates. Order I refers
to the setting where the arriving customers enter queue 1 first, then proceed to queue 2, and then exit the system. In the same manner, Order II
refers to the setting where the arriving customers enter queue 2 first, then proceed to queue 1, and then exit the system. Each entry in Table
\ref{Suresh-Whitt-Table} represents the minimum between two values: the absolute relative error of the respective method compared to the simulation
estimates; and the absolute error of the respective method, again compared to the simulation estimates.
\begin{table}
  \centering
    \begin{tabular}{c|c|c||c|c|c||c|c|c}
      \multicolumn{3}{c||}{} & \multicolumn{3}{|c||}{Order I} & \multicolumn{3}{|c}{Order II}\\
      \hline
      Case & $\rho_1$ & $\rho_2$ & QNET & QNA & LP-based & QNET & QNA & LP-based\\ \hline
      1 & 0.9 & 0.9 & 0.01 & 0.01 & 0.00 & 0.02 & 0.13 & 0.01\\
        & 0.8 & 0.8 & 0.08 & 0.08 & 0.08 & 0.06 & 0.08 & 0.07\\\hline
      2 & 0.9 & 0.9 & 0.03 & 0.06 & 0.03 & 0.01 & 0.21 & 0.01\\
        & 0.8 & 0.8 & 0.04 & 0.01 & 0.02 & 0.08 & 0.07 & 0.07\\ \hline
      3 & 0.9 & 0.9 & 0.04 & 0.02 & 0.04 & 0.04 & 0.19 & 0.05\\
        & 0.8 & 0.8 & 0.00 & 0.04 & 0.01 & 0.02 & 0.02 & 0.01\\ \hline
      4 & 0.9 & 0.9 & 0.13 & 0.30 & 0.14 & 0.19 & 0.29 & 0.19\\
        & 0.8 & 0.8 & 0.14 & 0.08 & 0.11 & 0.09 & 0.00 & 0.07\\ \hline
      5 & 0.9 & 0.9 & 0.01 & 0.12 & 0.01 & 0.01 & 0.01 & 0.00\\
        & 0.8 & 0.8 & 0.09 & 0.14 & 0.08 & 0.05 & 0.05 & 0.05\\ \hline
      6 & 0.9 & 0.9 & 0.01 & 0.08 & 0.02 & 0.01 & 0.10 & 0.01\\
        & 0.8 & 0.8 & 0.01 & 0.04 & 0.01 & 0.01 & 0.04 & 0.01\\ \hline
      7 & 0.9 & 0.9 & 0.04 & 0.18 & 0.05 & 0.06 & 0.07 & 0.04\\
        & 0.8 & 0.8 & 0.02 & 0.05 & 0.03 & 0.01 & 0.06 & 0.01\\ \hline\hline
      \multicolumn{3}{c||}{Average} & 0.05 & 0.09 & 0.04 & 0.05 & 0.09 & 0.04\\
    \end{tabular}
  \caption{\textbf{Summary of comparisons between QNA, QNET and the LP-based method for Example 4.} The
  table depicts the minimum between the absolute relative error and absolute error with respect to the simulation estimates, for all three methods,
  different orders of the queues, and different values of $\rho_1=\rho_2$.
  The last row depicts the average of the errors for the different combinations of methods/orders.}\label{Suresh-Whitt-Table}
\end{table}
As expected, the QNET and LP-based methods give more accurate estimates compared to the QNA method under balanced heavy traffic conditions. The
LP-based method is seen to perform as least as well as the QNET method, and gives rise to more accurate estimates in most cases.

\section{Extensions: Parameter Uncertainty and a Robust Formulation}\label{sec:robust}

Consider a network of 2 queues in tandem, i.e., the departure process from queue 1 forms the arrival process for queue 2. Jobs arrive to queue 1
according to an exogenous renewal process, and upon completion of processing at queue 2  they depart the system. We can characterize this simple
queueing network  by the 5-tuple
$$(C^2_a,C^2_{s_1},C^2_{s_2},\rho_1,\rho_2),$$
where $C^2_a$ is the squared coefficient of variation for the interarrival time distribution, $\rho_i$ is the mean service time at station $i=1,2$,
and $C^2_{s_i}$ is the squared coefficient of variation of the service time distribution at queue $i=1,2$. Assuming the arrival rate to queue 1 is
one, the value of $\rho_i$ coincides with the traffic intensity at queue $i=1,2$. The 2-dimensional workload process can be approximated by a
2-dimensional SRBM $Z$ in the orthant as given in Definition \ref{def:SRBM}. Following the prescription in Harrison and Nguyen
\cite{Harrison-Nguyen:90}, it can be shown that
$$\mu = \left(\begin{array}{c} \rho_1-1 \\ \rho_2/\rho_1-1\end{array}\right) \,, \quad  \sigma =
\left(\begin{array}{cc} \rho_1^2(C^2_a+C^2_{s_1}) & -\rho_1\rho_2 C^2_{s_1} \\ -\rho_1\rho_2 C^2_{s_1} & \rho_2^2(C^2_{s_1}+C^2_{s_2})
\end{array}\right)\,,\quad R = \left(\begin{array}{cc} 1 & 0 \\ -\rho_2/\rho_1 & 1\end{array}\right).$$
Using this, one can apply the LP-based algorithm to solve for the stationary distribution of $Z$ (see Example 4 in $\S$\ref{sec:app-SRBM}).

One potential issue that arises in practice is that the value of certain input parameters may not be known exactly. In what follows, we will
illustrate how such parameter uncertainty can be dealt with in the context of our algorithm, via a \emph{robust} formulation. To be concrete, suppose
we are interested in incorporating uncertainty with regard to the traffic intensities at each station. Assume that the true values of
$(\rho_1,\rho_2)$ belong to the following set
$$\mathcal{Q}_{\delta}:=\set{\hat{\rho_1}+\delta\rho_1,\hat{\rho_2}+\delta\rho_2\,\,:\,\, \delta\rho_1^2/a_1 + \delta\rho_2^2/a_2\leq r^2},$$
with $a_1$ $a_2>0$, and $a_1 r$, $a_2 r\ll1,$. That is, the values of $(\rho_1,\rho_2)$ are contained in an ellipsoid with ``small'' radius $r$ and
centered at $(\hat{\rho}_1,\hat{\rho}_2)$. This ellipsoidal uncertainty set is motivated by the classical robust counterpart formulation in Ben-Tal
and Nemirovski \cite{Bel-Nemirovski:99}. ``Translating'' this uncertainty to the SRBM data, we have
$$\mu = \left(\begin{array}{c} \hat{\rho}_1 + \delta\rho_1-1 \\ \hat{\rho}_2/\hat{\rho}_1 -(\hat{\rho}_2\delta\rho_1)/\hat{\rho}_1^2 + \delta\rho_2/\hat{\rho}_1-1\end{array}\right) \,,\quad
R = \left(\begin{array}{cc} 1 & 0 \\ -\hat{\rho}_2/\hat{\rho}_1 +(\hat{\rho}_2\delta\rho_1)/\hat{\rho}_1^2 - \delta\rho_2/\hat{\rho}_1 &
1\end{array}\right),$$
$$\sigma =
\left(\begin{array}{cc} (\hat{\rho}_1^2 + 2 \hat{\rho_1}\delta\rho_1) (C^2_a+C^2_{s_1}) & -(\hat{\rho}_1\hat{\rho}_2 + \delta\rho_1\hat{\rho}_2+
\hat{\rho}_1\delta\rho_2 ) C^2_{s_1} \\ -(\hat{\rho}_1\hat{\rho}_2 + \delta\rho_1\hat{\rho}_2+ \hat{\rho}_1\delta\rho_2 ) C^2_{s_1} & (\hat{\rho}_2^2
+ 2 \hat{\rho_2}\delta\rho_2) (C^2_{s_1}+C^2_{s_2})
\end{array}\right).$$
Note that both the second order differential operator $G$ as well as the boundary differential operator $D_i$ are linear in $\delta\rho_1$ and
$\delta\rho_2$. Consider now the BAR constraints in ($\mathcal{P}_{nm}$):
$$\left| \sum_{i:\,x_i\in S_n} (Gf)(x_i) \lambda_i + \sum_{j=1}^d \sum_{i:\,x_i\in F_j^n} (D_if)(x_i) \gamma_{ji} \right|\leq u,$$
for all $f\in \mathcal{F}_m$. Define $y = (\lambda,\gamma_1,\hdots,\gamma_d)$, where $\lambda \in \Real^{|S_n|}$ and $\gamma_i\in\Real^{|F^n_j|},$
$j=1,\hdots,d$. The BAR constraints can then be represented as
$$\left| A^{nm} y\right|\leq e u,$$ for some matrix $A^{nm}\,\in\,\Real^{|\mathcal{F}_m| \times |S_n|+\sum |F^n_j|}$, where $|a|$ denotes componentwise absolute value for a vector $a$, $e$ is a
$|\mathcal{F}_m|$-dimensional vector with all its entries equal to one, and $|\mathcal{F}_m|$ refers to the cardinality of the set of basis functions
spanning $\mathcal{F}_m$. It is easy to see that $A^{nm}$ is linear in $\delta\rho_1$ and $\delta\rho_2$, therefore it can be represented as
$A^{nm}:=\hat{A}^{nm} + \delta\rho_1 A^{nm}_{\delta\rho_1}+ \delta\rho_2 A^{nm}_{\delta\rho_2}$ for some matrices $\hat{A}^{nm}$,
$A^{nm}_{\delta\rho_1}$, $A^{nm}_{\delta\rho_2}$ $\in\,\Real^{|\mathcal{F}_m| \times |S_n|+\sum |F^n_j|}$ independent of $\delta\rho_1$ and
$\delta\rho_2$. As a consequence, when uncertainty with regard to traffic intensities is modeled as above, we have that $A^{nm}$ belongs to the
uncertainty set
$$ \mathcal{U}^{nm} = \set{\hat{A}^{nm} + \delta\rho_1 A^{nm}_{\delta\rho_1}+ \delta\rho_2 A^{nm}_{\delta\rho_2} \,:\, \delta\rho_1^2/a_1 + \delta\rho_2^2/a_2\leq r^2 }.$$
Following Ben-Tal and Nemirovski \cite{Bel-Nemirovski:99}, we consider the {\it robust counterpart} of ($\mathcal{P}_{nm}$)

$\hspace{-0.6cm}(\mathcal{R}_{n m})$

$\hspace{0.5cm}\min \quad\quad u$
\begin{eqnarray*}
 s.t. & & \left| A^{nm} y \right|\leq e u \quad \mbox{ for all } A^{nm}\in \mathcal{U}^{nm} \hspace{2.6cm} [\approx \mbox{BAR + uncertainty set}]\\
 & & \sum_{i:\,x_i\in S_n} \lambda_i = 1 \hspace{8.05cm} [\mbox{Normalization}]\\
 & & \sum_{i:\,x_i\in S_n} g(x_i)\lambda_i \leq M \hspace{5.0cm} [\mbox{Interior measure tightness}]\\
 & & \sum_{i:\,x_i\in F^n_i} g_i(x_i)\gamma_{ji} \leq M_j \quad \mbox{ for all } j\,\in \set{1\hdots,d} \hspace{0.3cm} [\mbox{Boundary measure tightness}]\\
 & & \lambda,\gamma_1,\hdots,\gamma_d \geq 0. \hspace{7.5cm} [\mbox{Nonnegativity}]
\end{eqnarray*}
With this robust formulation we are looking for a probability distribution and finite measures that minimize the violation of the BAR condition, for
all possible realizations of $(\rho_1,\rho_2)\in \mathcal{Q}_{\delta}$. Let us denote by $A_i$ the $i$-th row of a matrix $A$. The point $y =
(\lambda,\gamma)$ is ($\mathcal{R}_{nm}$)-feasible if it is ($\mathcal{P}_{nm}$)-feasible, and
\begin{equation*}
|(\hat{A}^{nm})_i\, y + \delta\rho_1 (A^{nm}_{\delta\rho_1})_i\, y + \delta\rho_2 (A^{nm}_{\delta\rho_2})_i\, y | \leq u \quad \mbox{for
}\delta\rho_1^2/a_1 + \delta\rho_2^2/a_2\leq r^2\, ,\, i=1,\hdots,|\mathcal{F}_{m}|.
\end{equation*}
Noting that
\begin{eqnarray*}
& & \max_{\delta\rho_1^2/a_1 + \delta\rho_2^2/a_2\leq r^2}|(\hat{A}^{nm})_i\, y + \delta\rho_1 (A^{nm}_{\delta\rho_1})_i\, y + \delta\rho_2
(A^{nm}_{\delta\rho_2})_i\, y| \\ & & \hspace{5cm} = |(\hat{A}^{nm})_i\, y| + \| \sqrt{a_1}r(A^{nm}_{\delta\rho_1})_i\, y,
\sqrt{a_2}r(A^{nm}_{\delta\rho_2})_i\, y \|_2,\\
\end{eqnarray*}
we can restate ($\mathcal{R}_{nm}$) as a \emph{conic quadratic program}, with the BAR constraint above replaced by:
\begin{equation*}
\left| (\hat{A}^{nm})_i\, y \right|\leq e u - \| b_1(A^{nm}_{\delta\rho_1})_i\, y, b_2(A^{nm}_{\delta\rho_2})_i\, y \|_2 \, \mbox{ for all }
i=1,\hdots,|\mathcal{F}_m|,
\end{equation*}
where $b_i :=\sqrt{a_i}r$, $i=1,2.$ Solving such conic-quadratic problems can be done in polynomial time using interior point methods, with
essentially the same computational complexity as LP problems of similar size (see \cite{Ben-Zibulevsky:97}).

\section{Extensions: An Illustrative Application to Diffusions Arising in the Many-server Heavy-traffic Regime}\label{sec:app-twoclassqueue}

In this section we illustrate how the proposed algorithm can be used to approximate the stationary distribution of diffusions arising in the so
called many-server (or Halfin-Whitt) heavy-traffic regime. For this purpose we will consider a two class queueing system operating under a first-in
first-out (FIFO) discipline, with class 1 receiving static priority over class 2. Specifically, class 1 arrivals, say customers, occur according to
Poisson process with rate $\lambda_1$. Requests for class 1 service engage one \emph{server} each for i.i.d. exponentially distributed amounts of
time with rate $\mu_1$, provided that the total number of class 1 customers in the system is less than the total number of servers $n$; otherwise,
they wait in queue for an available server. Class 2 customers are assumed to arrive according to a Poisson process with rate $\lambda_2$, independent
of the arrival process for class 1 customer. Requests for class 2 service are i.i.d. exponentially distributed with rate $\mu_2$. Class 2 customers
are served in FIFO manner, but are only processed  when servers are available \emph{and} there are no class 1 customers in queue. Despite its simple
structure, exact analysis of this system is not straightforward.

\subsection{Diffusion Approximation and BAR Condition}\label{sec:2-class-diffapprox}

Maglaras and Zeevi \cite{Maglaras-Zeevi:04} studied the system described above and derived a diffusion approximation under the many-server regime
developed by Halfin and Whitt \cite{Halfin-Whitt:81}. This regime is defined by letting the number of servers and arrival rates grow proportionally,
and concurrently letting the system utilization approach 1 at an appropriate rate. Specifically, consider a sequence of systems with $n$ servers and
arrival rates $\lambda_i^n := n\kappa_i\mu_i-\gamma_i \sqrt{n}\mu_i$ for some constants $\kappa_i>0$ and $\gamma_i\in \Real$ for $i=1,2$ such that
$\kappa_1+\kappa_2 =1$.

Let $Q^n_i(t)$ denote the number of customers in the system, both in service and in queue, in each customer class $i=1,2$. Put $ X_i^n(t) :=
(Q_i^n(t)-\kappa_i n)/\sqrt{n},\, i=1,2$. Then, if $(X^n_1(0),X^n_2(0))\rightarrow \xi\in\Real^2$, then $X^n=((X_1^n(t), X_2^n(t)) \,:\, t\geq 0)$
converges (as $n\rightarrow \infty$) to a diffusion process $X= (X(t)\,:\,t\geq 0)$ being the unique strong solution to
\begin{equation}\label{eq:SDE-two-class}
dX(t) = b(X(t))dt+\Sigma dB(t), \quad \quad X(0) =\xi.
\end{equation}
Here $B =(B(t):t\geq 0)$ is a standard Brownian motion in $\Real^2$, and the infinitesimal drift function $b_i(\cdot)$ for the $i$-th component is
given by
\begin{eqnarray*}
b_1(x) & := & -\mu_1 \gamma_1 - \mu_1 x_1 \, ,\\
b_2(x) & := & \left\{\begin{array}{lr} -\mu_2 \gamma_2 -\mu_2 x_2 & x_1 + x_2\leq 0\\
-\mu_2\gamma_2 + \mu_2 x_1 & x_1 + x_2 >0,\\ \end{array} \right. \\
\end{eqnarray*}
where $\Sigma := $ diag$(\sigma_1, \sigma_2$), with $\sigma_i^2 := 2\mu_i\kappa_i$, $i=1,2$. The first component of $X$ is an Ornstein-Uhlenbeck
process, and the second component has more complicated structure. The process $X$ admits a unique stationary distribution $\pi$ if and only if
$\gamma := \gamma_1+\gamma_2 >0$ (see Maglaras and Zeevi \cite{Maglaras-Zeevi:04} for further details). For $f\in C^2_b$ define the second order
differential operator
\begin{equation*}
Lf := \sum_{i=1}^2 \left(b_i \frac{\partial f}{\partial x_i} + \frac{\sigma^2_i}{2} \frac{\partial^2 f}{\partial x_i^2}\right).
\end{equation*}
Under the above stability condition, a direct application of Ito's lemma shows that for all $f\in C^2_b$, the stationary distribution $\pi$ must
satisfy the following BAR:
\begin{equation}\label{eq:bar-twoserver}
\int_{\Real^2} \left(\sum_{i=1}^2 b_i(x) \frac{\partial f(x)}{\partial x_i} + \frac{\sigma^2_i}{2} \frac{\partial^2 f(x)}{\partial x_i^2}
\right)\pi(dx) = 0.
\end{equation}
The sufficiency of (\ref{eq:bar-twoserver}) is established in the next proposition.
\begin{proposition}\label{prop:bar}
Suppose $\pi$ is a probability measure on $\Real^2$. If it satisfies (\ref{eq:bar-twoserver}), then $\pi$ is the stationary distribution for $X$.
\end{proposition}
\begin{proof}
Consider $f\, \in C^2_b$, and suppose that $x_0\in \Real^2$ is such that $\sup_{x \in \Real^2}f(x) = f(x_0)\geq 0$. Application of Ito' formula gives
\begin{equation*}
f(X(t)) - \int_0^t (Lf)(X(s))ds =  f(x_0) + \int_0^t \nabla f(X(s))\cdot \Sigma dW(s).
\end{equation*}
Taking expectations w.r.t. $\mathbb{P}_{x_0}$, we have
\begin{equation*}
\mathbb{E}_{x_0}[f(X(t))] - f(x_0) = \mathbb{E}_{x_0}\left[\int_0^t (Lf)(X(s))ds\right].
\end{equation*}
Divide by $t$ and take $t\rightarrow 0$. By the continuity of $Lf$ and the path continuity of $X$, we get $(Lf)(x_0)\leq 0$, since the left hand side
above is non-positive. Consider a probability measure $\pi$ on $\Real^2$ satisfying (\ref{eq:bar-twoserver}). Since $f(X(t)) - \int_0^t(Lf)(X(s))ds$
is a $\mathbb{P}_x$-martingale, the measure $\mathbb{P}_{\pi} \equiv \int_{\Real^2} \mathbb{P}_{x} \pi(dx)$ is a solution for the martingale problem
for $(L,\pi)$. Also, $C^2_b$ is an algebra, and is dense in $C(\Real^2$). As inidicated above, the operator ($L,C^2_b$) satisfies the positive
maximum principle, and therefore Echeverria's theorem (see \cite{Ehtier-Kurtz}, Theorem 9.14 of Chapter 4) applies, to assert that $\pi$ is a
stationary solution of the martingale problem for ($L,\pi$). We know that $\mathbb{P}_{\pi}$ is a solution of the martingale problem for ($L,\pi$),
and this solution is unique, since (\ref{eq:SDE-two-class}) admits a unique strong solution; see \cite{Maglaras-Zeevi:04}. This implies that $\pi$ is
a stationary distribution for $X$. The result follows from the uniqueness of the stationary distribution.\qed
\end{proof}

\subsection{LP Formulation and Convergence}
Adopting the approach used in $\S$\ref{sec:algorithm}, consider an increasing sequence of finite cardinality sets $S_n\subseteq \Real^2$ such that
$S_n\uparrow S$, and consider an increasing sequence $\set{\mathcal{F}_m}$ of finite subspaces of $C^2_b$ such that, for all $f\in C^2_b$ and for any
compact set $A\subseteq \Real^2$ there exists a $f_m\in\mathcal{F}_m$, such that $\lim_{m\rightarrow\infty}\sup_{x\in A} | (L(f_m-f))(x) |= 0$. For
fixed values of $n$ and $m$, the following LP mirrors the one in $\S$\ref{sec:LPformulation}:

$\hspace{0.25cm}(\mathcal{P'}_{n m})$

$\hspace{1.3cm}\min \quad\quad u $
\begin{eqnarray*}
 s.t. & & \left| \sum_{i:\,x_i\in S_n} (Lf)(x_i)\lambda_i\right| \leq u,\, \mbox{ for all } f\in \mathcal{F}_m \hspace{3.0cm} [\approx\mbox{BAR}]\\
 & & \sum_{i:\,x_i\in S_n} \lambda_i = 1 \hspace{6.3cm} [\mbox{Normalization}]\\
 & & \sum_{i:\,x_i\in S_n} g(x_i)\lambda_i \leq M \hspace{6.0cm} [\mbox{Tightness}]\\
 & & \lambda \geq 0, \hspace{7.4cm} [\mbox{Nonnegativity}]
 \end{eqnarray*}
for a suitably chosen $M< \infty$. In this formulation, $\set{\lambda_i\,:\,i \mbox{ s.t. }x_i\in S_n}$ is a discrete distribution with $S_n$ as its
support, that approximates $\pi$. The function $g:\Real^2\rightarrow\Real$ is chosen such that $g(x)\leq g(y)$ for $x\leq y$ componentwise, and such
that $g(x)\rightarrow\infty$ as $\norm{x}\rightarrow\infty$, and $M$ is set such that $\mathbb{E}_{\pi}(g(X))\leq M$. Convergence to $\pi$ along a
subsequence of solutions to these LP's can be proved following the same steps as in proof of Theorem \ref{theo:convergence}.

\subsection{Practical Considerations}
In order to implement and use the algorithm described above, we need to specify the sequences $\set{S_n}$ and $\set{\mathcal{F}_m}$, and a feasible
specification of the tightness condition, via choice of $M$ and $g$.

\textbf{Tightness.} For illustration purposes we will assume $\gamma_2>0$ to facilitate the explicit construction of a simple Lyapunov function.
Along the lines of $\S$\ref{sec:issues}, consider the function $f(x)= x_1^2+x_2^2$. For this function we have
\begin{eqnarray*}
(Lf)(x) &  =  &-\mu_1\gamma_1x_1 -\mu_1x_1^2 - (\mu_2\gamma_2x_2 +\mu_2x_2^2)\textbf{1}\set{x_1+x_2\leq 0}\\
& &  \hspace{2.8cm}- (\mu_2\gamma_2x_2 -\mu_2x_2x_1)\textbf{1}\set{x_1+x_2> 0} + \sigma_1^2 + \sigma_2^2,\\
\end{eqnarray*}
where $\textbf{1}\set{B}$ is the indicator of the set $B$. It then follows that $(Lf)(x) \leq -c_1 \|x\|_1 + c_2$ for all $x$ such that
$\|x\|^2_2\geq r$, for some $r>0$, and for some constants $c_1, c_2>0$. Therefore, we conclude that $ (Lf)(x) \leq -c_1 \|x\|_1 + c_3$ for all
$x\in\Real^2$, for a suitable constant $c_3\geq c_2$. It follows from Corollary 2 of Glynn and Zeevi \cite{Glynn-Zeevi:07} that $
\mathbb{E}_{\pi}\|X\|_1 \leq c_3/c_1,$ therefore we can take $g(x) = \|x\|_1$, and $M = c_3/c_1$.

\textbf{Distribution and function space approximation.} The first coordinate of the diffusion process $X$ has a normal distribution with mean
$\gamma_1$ and variance $\kappa_1$ as its stationary distribution. Based on this, we set the grid to have a ``Gaussian'' spacing on the first
coordinate, and in absence of prior information on the distribution of $X_2(\infty)$, we take the grid to be equispaced on the second coordinate:
\begin{eqnarray}\nonumber
S_n & = & \{ x\in S \,:\, \{ x_1 = (\Phi^{-1}(1/2+k/(2n))+\gamma_1)/\sqrt{\kappa_1} \quad k=0,\pm1\hdots,\pm n-1\}\\\label{eq:grid2dim}
 & & \times \{x_2
=\alpha j/n \quad j=-n,\hdots,n\}\},
\end{eqnarray}
for some constant $\alpha>0$. Regarding the approximation of the function space $C^2_b$, we take the same sequence of finite dimensional spaces
specified in $\S$\ref{sec:issues}.

\subsection{Numerical Results}
In the 2-class queue described in $\S$\ref{sec:2-class-diffapprox}, take $\kappa_1 =\kappa_2 = 0.5$ and $\gamma_1 =\gamma_2 =0.25$ (these parameters
correspond to the example studied in Maglaras and Zeevi \cite{Maglaras-Zeevi:04}). We take $S_n$ as in (\ref{eq:grid2dim}) with $n=100$ and degree of
polynomials $m=6$. Table \ref{table:two-class-2} shows that the algorithm provides reasonable moment approximations for each coordinate. There,
benchmark moments for the second coordinate were computed via Monte Carlo simulation, using an Euler scheme to approximate (\ref{eq:SDE-two-class})
and arrive at a density estimate from which the moments were derived.
\begin{table}
  \centering
    \begin{tabular}{c||c|c|c|c||c|c|c|c}
       & $\mathbb{E}[X_1^1]$ & $\mathbb{E}[X_1^2]$ & $\mathbb{E}[X_1^3]$ & $\mathbb{E}[X_1^4]$& $\mathbb{E}[X_2^1]$ & $\mathbb{E}[X_2^2]$ & $\mathbb{E}[X_2^3]$ & $\mathbb{E}[X_2^4]$ \\\hline \hline
      True/Simulated value     & -0.25 & 0.56 & -0.39 & 0.94 & 1.19 & 8.31 & 70.67 & 854 \\
      LP-based estimate        & -0.25& 0.56 & -0.39 & 0.94 & 1.26 & 7.79 & 69.68 & 812 \\
    \end{tabular}
\caption{\textbf{Moment approximations for first and second coordinates of the 2-dimensional diffusion}. The algorithm uses grid size $n=100$ and
$m=6$ degree polynomials.}\label{table:two-class-2}
\end{table}
Figure \ref{figure:two-class} compares the approximation of the marginal distributions against the underlying marginal distributions when using grid
size $n=50$ and $m=4$ degree polynomials.
\begin{figure}
\centering
\includegraphics[width=18cm]{N-figure4.eps}
\caption{\textbf{Marginal distribution estimates for the two-dimensional diffusion}. The solid line is computed via Monte Carlo simulation, and the
dotted line represents the algorithm estimates based on grid size $n=50$ and $m=4$ degree polynomials.}\label{figure:two-class}
\end{figure}
It is evident from Figure \ref{figure:two-class} that the LP-based estimate of the stationary distribution, while giving rise to reasonable moment
estimates, does not exhibit the smoothness that the true distribution is anticipated to posses (see the solid lines in Figure
\ref{figure:two-class}). To try to remedy this, we incorporate ``smoothness'' constrains into our LP formulation as follows. Consider a
\emph{neighborhood parameter} $r>0$, and a \emph{smoothness parameter} $\delta$. The following Lipschitz continuity constraint can then be added to
($\mathcal{P}_{nm}$):
\begin{equation}\label{eq:smoothness}
\abs{\pi^n_m(x) - \pi^n_m(y)}\leq \delta \quad \mbox{ for all } \, (x,y)\in S_n \quad \mbox{ s.t. }\, \|x-y\|_2<r.
\end{equation}
In practice we introduce these constraints on each marginal distribution separately. We choose $r$ so that the $r$-ball contains only one neighbor in
each marginal, and $\delta$ sufficiently small to exclude ``peaks'' as those observed in Figure \ref{figure:two-class}. In particular, setting
$\delta_1 =0.0026$ (the first coordinate) and  $\delta_2 =0.0011$ (the second coordinate) we obtain the distribution estimates reported in Figure
\ref{figure:two-class-2}.
\begin{figure}
\centering
\includegraphics[width=18cm]{N-figure3.eps}
  \caption{\textbf{Smoothed marginal distribution estimates for the two-dimensional diffusion}. The solid line is computed via Monte Carlo simulation, and the
  dotted line represents the algorithm estimates based on grid size $n=50$ and $m=4$ degree polynomials, incorporating smoothness
constraints.}\label{figure:two-class-2}
\end{figure}
Despite the smoothness imposed through the LP formulation, moment estimates do not differ significantly compared to those reported in Table
\ref{table:two-class-2}.

\begin{remark}
{\rm Assuming that $\pi$ possesses a density $\xi$ with bounded derivative, one can implement constraint (\ref{eq:smoothness}) by deriving an upper
bound on the Lipschitz constant. This bound can be obtained as follows. For any measurable set $A\subset\Real^d$, we have that $\pi(A) =
\int_{\Real^d} K(y,A) \pi(dy)$, where $K(A,y)$ denotes the transition kernel of the diffusion process. Then, for the density $\xi$ we have that
$\xi(x) = \int_{\Real^d} K(y,x) \xi(y) dy $ and hence one expects that
\begin{equation*}
\|\nabla\xi(x)\|_1 \leq  \int_{\Real^d} \|\nabla K(y,x)\|_1 \xi(y) dy = \int_{\Real^d} \norm{\nabla K(y,x)} \pi(dy).
\end{equation*}
Suppose we can find a twice continuously differentiable function $g:\Real^d\rightarrow \Real_+$ and a constant $c$ for which
$$ (Lg)(x) \leq -\norm{\nabla K(y,x)} + c \,, \mbox{ for all }x\in\Real^d,$$
uniformly for $y\in \Real^d$. Then, by Corollary 2 of Glynn and Zeevi \cite{Glynn-Zeevi:07} we have that
$$ \|\nabla\xi(x)\|_1 \leq \int_{\Real^d} \norm{\nabla K(y,x)} \pi(dy) \leq c.$$
}
\end{remark}

\section{Extensions: On More General Approximating Spaces}\label{sec:prob-meas-choice}

Consider the LP formulation presented on \S \ref{sec:LPformulation}. There, for fixed values of $n$ and $m$, a feasible solution $(\lambda,\gamma)$
to $(\mathcal{P}_{n m})$ is used to approximate $\pi$ and $\nu:=(\nu_1,\hdots,\nu_d)$ through a probability measure, and finite measures, assigning
mass to a finite and discrete set of points on $S$ and $F_i$ $i=1,\hdots,d$, respectively. More specifically,
\begin{equation*}
\pi^{n}_{\lambda}(x) := \sum_{i:\,x_i\in S_n} \lambda_i \delta_{x_i}(x), \quad\quad\nu^{n}_{\gamma j}(x) := \sum_{i:\,x_i\in F_j^n} \gamma_{ji}
\delta_{x_i}(x)\,\,\quad j=1,\hdots,d,
\end{equation*}
where $\delta_y(x)=1$ if $x=y$ and $0$ otherwise. This point-mass-based approximation can be generalized as follows. Let $\M$ denote the space of
probability measures on $S$, and $\M_j$ the space of finite measures on $F_j \,,\, j=1\hdots,d$. For arbitrary sets $\widetilde{\M}\subset\M$ and
$\widetilde{\M}_j\subset\M_j$, $j=1,\hdots, d$ of finite cardinality $\abs{\widetilde{M}}$, $\abs{\widetilde{M}_j}\leq\infty$, and vectors
$\lambda\in \Real^{|\widetilde{\M}|}$ and $\gamma:=(\gamma_1,\hdots,\gamma_d)$, $\gamma_j \in \Real^{|\widetilde{\M}_j|}$, $j=1,\hdots, d$, define:

\begin{equation}\label{eq:new-architecture}
\pi(\widetilde{\M},\lambda) := \sum_{i:\,\phi_i \in \mathcal{\widetilde{M}}} \lambda_i \phi_i, \quad\quad\nu_j(\widetilde{\M}_j,\gamma_j) :=
\sum_{i:\,\varphi_i\in \mathcal{\widetilde{M}}_j} \gamma_{ji} \varphi_i\,\,\quad j=1,\hdots,d.
\end{equation}

Where $\phi_i$ and $\varphi_i$ are arbitrary elements of $\widetilde{M}$ and $\widetilde{M}_j$, respectively. Consider a sequence $\set{\M_n}$ such
that $\M_n\subseteq \M$, $|\M_n|< \infty$, for all $n\geq1$, and: i) for all $\phi\in\M$, there exist a sequence of vectors $\set{\lambda_n}$ such
that $\int_S g(x)\pi(\M_n,\lambda_n)(dx)\leq \int_S g(x)\phi(dx)$, and $\sum_i\lambda_{ni}=1$ for all $n\geq1$, where $g:S\rightarrow\Real_+$ is
specified in \S\ref{sec:LPformulation}; and ii) $\pi(\M_n,\lambda_n)\Rightarrow \phi$. Analogously, consider sequences $\set{\M^n_j}$ $j=1,\hdots,d$,
such that $\M_j^n\subseteq \M_j$, $|\M_j^n|< \infty$, $j=1,\hdots,d$, for all $n\geq1$, and:  i) for all $\varphi\in\M_j$ there exist a sequence of
positive (finite) vectors $\set{\gamma^n_j}$ such that $\int_{F_j} g_j(x)\nu_j(\M_j^n,\gamma^n_j)(dx)\leq \int_{F_j} g_j(x)\varphi(dx)$ for all
$n\geq 1$, where $g_j:F_j\rightarrow\Real_+$ are specified in \S\ref{sec:LPformulation}; and ii) $\nu_j(\M_j^n,\gamma^n_j)\Rightarrow \varphi$, for
$j=1,\hdots,d$.

For any continuous function $f$ and $\phi \in \M$, define $ B(\phi,f) :=\int_S f(x)\, \phi(dx)$, and for $\varphi\in\M_j$ define $B_j(\varphi,f):=
\int_{F_j} f(x)\, \varphi_j(dx)$. Fixing $n\geq 1$ and $m\geq1$, we can restate $\mathcal{P}_{n m}$ as follows:

$\hspace{-0.6cm}(\widetilde{\mathcal{P}_n}_m)$

$\hspace{0.1cm}\min \quad\quad u $
\begin{eqnarray}\label{eq:new-Discrete-Bar}
 \mbox{s.t.} & & \left| \lambda_i B(\pi(\M_n,\lambda),Gf) + \sum_{j=1}^d
 B_j(\nu_j(\M_j^n,\gamma^n_j),D_jf)   \right|\leq u\,\, ,\,\, \mbox{for all } f\in \mathcal{F}_m \hspace{0.4cm} [\approx \mbox{BAR}]\\\label{eq:new-prob-dist}
 & & \sum_{i: \, x_i\in S_n} \lambda_i = 1 \hspace{8.75cm} [\mbox{Normalization}]\\\label{eq:new-tightness-dist}
 & & B(\pi(\M_n,\lambda_n),g) \leq M \hspace{5.25cm} [\mbox{Interior measure tightness}]\\\label{eq:new-tightness-boundary}
 & & B_j(\nu_j(\M_j^n,\gamma^n_j),g_j) \leq M_j \,, \mbox{ for all } j\in\set{1\hdots,d} \hspace{0.65cm} [\mbox{Boundary measure tightness}]\\\label{eq:new-nonnegative}
 & & \lambda_i,\, \gamma_{ij} \geq 0 \,\,,\, \mbox{ for all }\set{i:\,x_i\in S_n} \mbox{ and } j\in\set{1,\hdots,d}.\hspace{2.0cm} [\mbox{Nonnegativity}]
 \end{eqnarray}
Where the constants $M$, $M_j$, $j=1,\hdots,d$, are chosen according to the logic described in \S\ref{sec:LPformulation}. Note that
$\widetilde{\mathcal{P}_n}_m$ remains a finite dimensional linear program. The following corollary follows straightforwardly using the same arguments
used to prove Theorem \ref{theo:convergence}.
\begin{corollary}\label{theo:new-convergence}
Consider an SRBM $Z$ as in Definition \ref{def:SRBM}. Let $\pi$ be the stationary distribution for $Z$, with boundary measures $\nu_j\,,
\,j=1,\hdots,d$. Then, there exist sequences $\set{n_k}$ and $\set{m_k}$, $k=1,2,\hdots$, of positive integers that diverge to infinity as
$k\rightarrow\infty$, such that $\pi_{\lambda^{k}}^{n_km_k}\Rightarrow \pi$ and $\nu^{n_km_k}_{\underline{\gamma}^{k}j}\Rightarrow \nu_j$ as
$k\rightarrow \infty$, for $j=1,\hdots,d$, where $(\lambda^k, \underline{\gamma}^k)$ are optimal solutions to the LP
($\widetilde{\,\mathcal{P}\,}_{n_km_k}$).
\end{corollary}

The complexity of formulating $\widetilde{\,\mathcal{P}\,}_{n_km_k}$ depends directly on the ability to efficiently compute $B(\phi,f)$ for all pairs
$(\phi,f)\in \M_n\times\mathcal{F}_m$, and $B_j(\varphi,f)$ for all pairs $(\varphi,f)\in \M_j\times\mathcal{F}_m$, $j=1,\hdots,d$. For the
implementations presented in \S\ref{sec:app-SRBM} and \S\ref{sec:app-twoclassqueue}, the selection of the probability measure sequences as point
masses ``decouples'' it from the choice of $\mathcal{F}_m$, since computing $B$ and $B_j$, $j=1,\hdots,d$, amounts to evaluating the functions in
$\mathcal{F}_m$ at points in $S$ and $F_j$, $j=1,\hdots,d$, respectively. The following examples indicate other computationally tractable choices.

\noindent \textbf{Example 1.} Consider the algorithm defined in \S\ref{sec:app-SRBM}. Let $\M_n$ be a finite cardinality set of truncated Gaussian
distributions with support in the positive orthant. Let $\set{\mathcal{F}_m}$ be the same space of polynomials of degree less than $m$. Then,
$B(\phi,Gf)$ is nothing but a linear combination of moments of a truncated normal random variable, so it can be computed efficiently. Similarly, one
has a computationally efficient approximation of the boundary measures.

\noindent \textbf{Example 2. } Consider any affine diffusion in $\Real^d$, such that BAR is necessary and sufficient to characterize its unique
stationary distribution. Here one can take $\set{\M_n}$ such that $\M_n$ contains a finite number of Gaussian distributions, for all $n\geq 1$. If
$\set{\mathcal{F}_m}$ is chosen to be the space of polynomials of degree less than $m$, for all $m$, then $B(\phi,Gf)$ is, again, a linear
combination of moments of normal random variables, and can be computed efficiently.


%===================================================================
\bibliographystyle{plain}
\bibliography{SRBM}
%===================================================================

\appendix
\section{Numerical Results}\label{sec:numerical}
This appendix  collects further and more detailed numerical results, a summary of which was given in the main body of the paper. In particular,
tables \ref{Suresh-Whitt-Case1-Table} to \ref{Suresh-Whitt-Case7-Table} give detailed expected steady-state waiting time estimates for the experiment
described on \S\ref{sec:app-SRBM}, Example 4. There, simulation estimates computed by Suresh and Whitt \cite{Suresh-Whitt:90} are compared against:
(i) estimates computed using a software package described by Whitt \cite{Whitt:83}, dubbed the QNA method; (ii) estimates computed using the
algorithm proposed in Harrison and Dai \cite{Dai-Harrison:92}, which they refer to as the QNET method; and (iii) estimates computed using our
LP-based algorithm.
\input{Appendix}
\clearpage


\end{document}
