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

% Assaf definitions
\input{mydefinitions.tex}

% My definitions
\newcommand{\Field}{\mathcal{F}}
\newcommand{\Real}{\mathbb{R}}
\newcommand{\RDPlus}{\Real^d_{+}}
\newcommand{\RD}{\Real^d}
\newcommand{\Prob}{\mathbb{P}}
\newcommand{\set}[1]{\left\{#1\right\}}
\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 Heavy Traffic Approximations}
\author{{\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.zeevi@columbia.edu}}
\\Columbia University
\and {\sf Denis  Saure}\thanks{Columbia Graduate School of Business, e-mail: {\tt
dsaure05@gsb.columbia.edu}}
\\Columbia University}

\date{January  2008}
\begin{document}

\maketitle

\begin{abstract}
This paper proposes a linear programming algorithm for computing the stationary distribution of
multidimensional diffusions that arise as approximate models of queueing networks. Our algorithm is
based on the necessity and sufficiency of a Basic Adjoint Relationship (BAR), which characterizes
the stationary distribution. To use the BAR characterization , we approximate the state space with
a finite grid of points, and use a finite set of ``test'' functions. With these approximations, the
BAR reduces to a set of linear equations that can be solved using standard linear programming
techniques. We show that the sequence of stationary distributions which arise as solutions to the
algorithm converges, in a suitable sense, to the stationary distribution of the original diffusion
processes. We present applications to diffusions arising under the standard heavy-traffic regime (
Semi-martingale Reflected Brownian Motions) , and to diffusions arising under the Halfin-Whitt
heavy-traffic regime. Extensive computational experiences are reported. The algorithm is shown to
produce good estimates of the stationary moments, as well as the entire distribution.
\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:} a LP Algorithm for S-S Distribution of Heavy Traffic Approximations\\
%\vspace*{1em}

\noi {\bf Keywords:} queueing, diffusion approximations, reflected Brownian motion,
 Halfin-Whitt regime, numerical methods.\\
%\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}

Multiclass queueing networks are widely used as mathematical models of many real-world complex
systems, e.g., communications networks, manufacturing system, and service operations. That said,
most relevant queueing network models are rather intractable for purposes of exact performance
analysis. Consequently, considerable research efforts have been placed on deriving performance
bounds for these networks, the typical goal being to bound moments of the steady-state queue
lengths or workload [see, e.g., Kumar and Kumar \cite{Kumar-Kumar:94}, Sigman and Yao
\cite{Sigman-Yao:93}, Bertsimas, Gamarnik and Tsitsiklis \cite{Betsimas-etal:01}, Gamarnik and
Zeevi \cite{Gamarnik-Zeevi:06}, as well as references therein].

An alternative path is to approximate the queueing network model with a more simplified structure.
The idea then is to use the stationary distribution of the approximating model, whose dynamics are
``close'' to the ones of the original network as a proxy for the original model.Perhaps the most
prevalent approach in this context has focused on the study of diffusion approximations that arise
under the so-called heavy-traffic operating regimes. The main objective of this paper is to provide
a method to approximate the steady-state distribution of these diffusion approximations.

Consider a  diffusion process $X= \set{X(t),t\geq 0}$ with state-space $S\subset \Real^k$ that
arises as a heavy-traffic limit for some queueing network, and assume that $X$ admits a unique
stationary distribution $\pi$. In most cases, $\pi$ it is impossible to compute $\pi$ in closed
form, and a plausible alternative is to compute it numerically. To this end, let $\mathcal{A}$
denote the generator of the process $X$, and let  $\mathcal{D}_{\mathcal{A}}$ denote its domain.
Under suitable technical conditions on  $\mathcal{A}$ and $S$, we have that
\begin{equation}\label{eq:bar-intro}
\int_S \left(\mathcal{A}f\right) \pi(dx) = 0 \quad  \mbox{for all} f\in \mathcal{D}_{\mathcal{A}}
\end{equation}
is necessary and sufficient to characterize $\pi$. The equation above is known as the {\it basic
adjoint relationship} (BAR). In the context of queueing and heavy-traffic theory, it was first used
by Harrison and Williams \cite{Harrison-Williams:87a} to characterize the stationary distribution
of a reflected Brownian motion in the nonnegative orthant. Since (\ref{eq:bar-intro}) characterizes
$\pi$ as the solution to an infinite dimensional system of linear equations, it is impractical to
try to solve this directly. An alternative involves solving approximating $S$ and
$\mathcal{D}_{\mathcal{A}}$ by suitable chosen sequences $\set{S_n}$ and $\set{\mathcal{D}_m}$ of
finite subsets of $S$ and finite subspaces of $\mathcal{D}_{\mathcal{A}}$ respectively, such that
$S_n\uparrow S$ and $\mathcal{D}_m\uparrow\mathcal{D} = \mathcal{D}_{\mathcal{A}}$ in a precise
sense. For fixed values of $n$ and $m$, (\ref{eq:bar-intro}) can be stated as a finite-dimensional
linear program (LP), that can be solve efficiently. The hope is then, that $\pi$ arises as a limit
of a suitable sequence of optimal solutions to the aforementioned LP's. The main contribution of
this papers is to articulate mathematical conditions under which this holds true, for the classes
of diffusions processes arising in common heavy-traffic approximations, and to investigate the
performance of the algorithm numerically, contrasting it with existing methods.

\textbf{Related literature.} A related algorithm was proposed by Dai and Harrison \cite{Dai:90} for
solving stationary distribution of Semimartingale Reflected Brownian Motions (SRBM). In there, the
authors view the BAR as an orthogonality condition between a infinite dimensional functional space
and the stationary distribution. By consider an increasing sequence of finite dimensional
approximations of the infinite dimensional functional space, they obtain a sequence of
approximating densities by means of orthogonal projections into a sequence of finite dimensional
spaces. They prove that this sequence of densities converges to the stationary density of the SRBM,
but their proof relies on a certain conjecture that concerns the behavior of solutions to the BAR.
Later, Shen et al. \cite{Shen-Chen:02} considered a variant of the Dai-Harrison algorithm for the
case of a SRBM on a hypercube using a finite element method (or piecewise polynomials) to form the
finite dimensional approximations of the functional space of ``test'' functions. Shen and Chen
\cite{Shen-Chen:03} extend this algorithm for SRBM on the positive orthan. Our work differs from
these previous algorithms along three dimensions: First, it allows for one to impose further
constraints on the structure of the underlying stationary distribution, which results in reasonable
estimates not only for stationary moments but also for the distribution itself. Since the tails  of
the queue-length/work-load distribution play important roles as performance indicators, obtaining
good approximations to the distribution itself is of obvious value. Second, the convergence of our
algorithm does not rely on the structural conjecture that appears in the Dai-Harrison and follow-up
papers. Finally, our method can be easily extended to derive approximations to the stationary
distribution under parameter uncertainty using robust optimization (see section \ref{sec:robust}).

%Reiman's seminar paper \cite{Reiman:84} prove that, in the context of open-single-class
%queueing-networks, as the traffic intensity approaches one, the normalized vector queue length
%processes converges in distribution to a Reflected Brownian motion (RBM) on the positive orthant,
%$S = \mathbb{R}^d_+$. This diffusion was first constructed in earlier work of Harrison and Reiman
%\cite{Harrison-Reiman:81}.  Similar heavy-traffic limit results had been proven for various types
%of open d-station networks with possibly multiple customer types [see, e.g., Harrison and Williams
%\cite{Harrison-Williams:87a}, Harrison and Reiman \cite{Harrison-Reiman:81}].

%Another type of heavy-traffic regime was introduced by Halfin and Whitt \cite{Halfin-Whitt:81} in
%the context of a $M/M/C$ queue when letting the $C$ grow large and concurrently letting the traffic
%intensity approach 1 at an appropriate rate. Under this regime, the number of users in the system
%can be expressed roughly as $C+\sqrt{C}X(t)$ for $C$ large enough, where $X(t)$ is a simple
%dimensional diffusion process. There exists an extensive literature providing similar results for
%many-servers multiple customer types systems [see Mandelbaum, Reiman and Massey
%\cite{Reiman-Mandelbaum-Massey:98}, as well as references therein].

Harrison and Nguyen [CITE] were among the first to propose the use of the stationary distribution
of the diffusion model as a mean to approximate the one corresponding to the original network,
mostly on the context of moment calculation. It should be noted, however that even if the
approximating diffusion model is arrived at rigourously (i.e., as a formal heavy-traffic limit),
its stationary distribution may not provide a rigorous approximation of that of the underlying
model. To date, the validity of this interchange-of-limits has not been established in full
generality (see Gamarnik and Zeevi \cite{Gamarnik-Zeevi:06} and Gurvich and Zeevi
\cite{Gurvich-Zeevi} to some analysis supporting this interchange argument).

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 space Markov processes. Hernandez and
Lasserre \cite{Hernandez-Lasserre} extend this to analyze the convergence of linear-programming
approximations for discrete time controlled Markov processes in metric spaces. In
\cite{Hernandez-Lasserre}, the authors approximate a discrete time analogous of the BAR using a
discrete probability distribution and a finite subspace of test functions. The main objective there
is not to characterize the steady-state distribution but rather to minimize a steady-state cost
function. Mendiondo and Stockbridge \cite{Mendiondo-stockbrige:98} extends this work to the
continuous time setting in the context of long-term average and discounted control problems, when
the state space and control space are assumed to be compact. The key features that distinguish our
approach relative to these papers are: (i) Our focus is on computing the steady-state distribution
rather than a long-run average cost;(ii) Since our focus is on diffusions arising as heavy-traffic
limits, the state space is not necessarily compact, which creates a potential problem with  the
sequence of approximating measures; (iii) we explicitly show our algorithm's application in the
case of diffusions with reflecting boundaries, and (iv) feasibility of each linear program in our
sequence of approximating problems is guaranteed, while the work cited above assumes feasibility in
the respective measure-control space.

\textbf{The remainder of the paper.}  Section \ref{sec:algorithm} formulates our algorithm in
 general form, and its convergence is established under general conditions. Section
\ref{sec:app-SRBM} presents the application of the algorithm to the SRBM, which arises as an
approximating diffusion model under the classical heavy-traffic scaling, and Section
\ref{sec:app-twoclassqueue} illustrates its application to a Ornstein-Uhlenbeck type diffusion that
arises in the so-called many server heavy-traffic regime. Extensive computational experiences are
reported for both applications. Finally, Section \ref{sec:robust} extends the algorithm for the
case of uncertain parameters on the underlying queueing system, by formulating the ``robust''
counterpart of the LP approximation. Computational experiences are reported for the SRBM case.


\section{Proposed Algorithm}\label{sec:algorithm}
\input{BAR.tex}
\section{Applications: Semimartingale Reflected Brownian Motions}\label{sec:app-SRBM}


This section illustrates the application of the algorithm to a class of diffusion processes that
plays a central role in queuing theory: The Semimartingale Reflected Brownian Motions (abbreviated
as SRBM). These processes have been shown to arise as diffusion limits of open multiclass queueing
networks operating under ``conventional'' Heavy-Traffic conditions: Consider a sequence of queueing
networks indexed by $n$. The key idea of conventional Heavy-Traffic theory is to scale amounts to
express time in multiples of $n$ and space in multiples of $n^{1/2}$ so that functional central
limit theorem applies to a properly normalized queue length process while the sequence of traffic
intensities converges to 1 as $n$ increases. A formal for the SRBM definition is the following:

\begin{definition} Let $S\equiv\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$, a 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, SRBM behaves like Brownian motion in the interior of $S$, and is confined to the
orthant by instantaneous ``reflection'' at the boundary faces, where the direction of  reflection
is dictated by the matrix $R$. The most general condition currently known to ensure existence and
uniqueness (in law) of SRBM in the orthant is the matrix $R$ to be \emph{completely S}.

\begin{definition}
A $d\times d$ matrix R is said to be \emph{S} if there exits a d-dimensional vector $u\geq 0$ such
that $R u>0$, and to be a \emph{completely S} matrix if each of its principal submatrices is an
\emph{S} matrix.
\end{definition}

This \emph{completely S} condition is in fact necessary (Reiman and Williams
\cite{Reiman-Williams:88}) and sufficient (Taylor and Williams \cite{Taylor-Williams:93}).
Regarding its stationary distribution it has been shown (Dupuis and Williams
\cite{Dupuis-Williams:94}) that a sufficient condition for its existence and uniqueness is that all
solutions of an associated deterministic Skorohod problem are attracted to the origin in finite
time. Define $\gamma = R^{-1}\mu$. A more tractable condition, $\gamma<0$ is known to be necessary
when $R^{-1}\geq 0$ (Harrison and Williams \cite{Harrison-Williams:87a}, Dai \cite{Dai:90})  but
not sufficient (Dai and Harrison \cite{Dai-Harrison:07}).

If we want to apply our algorithm to this class of processes we need to check that (i) The BAR
condition is sufficient and necessary, and that (ii) we can provide a tightness bound for the
interior distribution and tightness and finiteness bounds for boundary measures.

\subsection{Bar Condition}
For $f\in C^2_b(S)$, define
\begin{equation}
Lf  \eqq \frac{1}{2} \sum_{i=1}^d\sum_{j=1}^d \sigma_{ij} \frac{\partial^2 f}{\partial x_i\partial
x_j} + \sum_{i=1}^d\mu_i \frac{\partial f}{\partial x_i}
\end{equation}
\begin{equation}
D_i f \eqq R_{\cdot i} \nabla f\, \mbox{  for }x\in F_i\, (i=1,\hdots,d)
\end{equation}

For these processes, BAR takes the following form:
\begin{equation}
\int_{\pi}L f d\pi + \sum_{i=1}^d \int_{F_i} D_i f d\nu_i=0 \quad \forall \,f \in
C^2_b(S)\label{eq:bar}
\end{equation}
where $\pi$ is the stationary distribution for $Z$ associated with boundary measures $\nu_i\,
(i=1,\hdots,d)$.

Necessity of BAR was first derived by Harrison and Williams \cite{Harrison-Williams:87a} when the
matrix $R$ is Minkowski ($I-R\geq 0$ and $I-R$ is transient), and later by Dai \cite{Dai:90} for
the \emph{completely S} case. Sufficiency was proven by Dai and Kurtz \cite{Dai-Kurtz:Unpublished}
through the following theorem:

\begin{theorem} Assume  R is a \emph{completely S} matrix. Suppose that $\pi_0$ is a
probability measure on S with support in the interior of S, and $\pi_1,\hdots,\pi_d$ are positive
finite measures with supports on $F_1,\hdots,F_d$ respectively. If they jointly satisfies
(\ref{eq:bar}), then $\pi_0$ is the stationary distribution for a $(\sigma,\mu,R)$-SRBM Z.
\end{theorem}

\subsection{Tightness condition}

The use of Lyapunov functions to bound expectations of Markov processes is a widely used technique
(for a clear exposition see Glynn and Zeevi \cite{Glynn-Zeevi:07} ). Given a function $f$ for which
one is interested in computing stationary moments, the key idea is to find a Lyapunov-like function
$g$ that satisfies a certain ``mean'' drift inequality with respect to $f$, which in turn leads to
bounds on the stationary expectation of the later. This idea has been used to prove the existence
of exponential moments of the interior stationary distribution for the SRBM: Glynn and Zeevi
\cite{Glynn-Zeevi:07} provided an explicit construction for the case of $R$ being symmetric and
positive definite. Budhiraja and Lee \cite{Budhiraja-Lee:06} have proven the existence of such
function for the \emph{completely S} case. We prove that same result holds for the boundary
measures. In what follows we will use the following proposition, proven by Dai \cite{Dai:90}
\begin{proposition}\label{prop:boundary}
If $\pi$ is the stationary distribution, associated with boundary measures $\nu_i\,
(i=1,\hdots,d)$, then
\begin{equation}
E_{\pi}\left[ \int_0^t f(Z(s)) dY_i(s)\right] = t \int_{F_i}f \nu_i, \quad (i=1,\hdots,d)
\end{equation}
where $\sigma_i$ denote the (d$-$1)-dimensional Lebesgue measure on the face $F_i$.
\end{proposition}

\begin{theorem}\label{thm:finite-boundary}
Assume R is a \emph{completely S} matrix, and suppose $\pi$ is the stationary distribution for $Z$
associated with boundary measures $\nu_i\, (i=1,\hdots,d)$ . There exits a vector $\upsilon>0$ such
that
\begin{equation*}
\int_{F_i} e^{\upsilon^tx}\, \nu_i(dx) <\infty \quad , i=1,\hdots,d
\end{equation*}
\end{theorem}
These results indicates that moment bounds exists, but we still need practical bounds to be used in
the setting of the algorithm, for tightness and finiteness of the boundary measures.

\subsection{Practical Bounds}


\subsubsection*{Finiteness of boundary measures} We will use a Lyapunov-function argument. Consider a
d-dimensional vector $v>0$ such that $R^tv>0$ and set $f(x) = x^tv$. Proceeding as in proof of
Theorem \ref{thm:finite-boundary} (see appendix) we have:

$$E_x[f(Z(t))]- E_x[f(x)] + \alpha  =  E_x[\sum_{i=1}^d \beta_i \int_0^t dY_i(s)]$$
where $\alpha \equiv v^t\mu>0$ and $\beta_i\eqq v^tR_{\cdot i}>0\quad i=1,\hdots,d$. Taking
expectations with respect to $\pi$, and using Proposition \ref{prop:boundary} we have
$$ \int_{F_i} \nu_i(ds) \leq \frac{\alpha}{\beta_i}\quad \,i=1\hdots,d$$

\subsubsection*{Tightness of stationary distribution} Here we will assume that $R$ is symmetric to
facilitate the explicit construction of a simple Lyapunov function. This can be relaxed by means of
a more clever choice of such function. Consider $f(x) = x^tR^{-1}x$. Proceeding as in the proof of
Theorem \ref{thm:finite-boundary}, we have

$$-f(x)+  E_x[\int_0^t Z(s)^t \gamma ds] \leq C + E_x[\sum_{i=1}^d \int_0^t Z_i(s) dY_i(s)]$$
where $C = \sum_{i j} \sigma_{ij}R_{ij}$. Notice that $E_x[\sum_{i=1}^d \int_0^t Z_i(s) dY_i(s)] =
0$ by definition of the SRBM. Taking expectations with respect to $\pi$ we have
$$\int_S \langle e, s \rangle \pi(ds) \leq \frac{C}{\min_j\set{\gamma_j}}$$

For the Minkowski case, a bound can be derived directly from lemma 8.4 in Harrison-Williams
\cite{Harrison-Williams:87a}.

\subsubsection*{Tightness of boundary measures} Take a d-dimensional vector $v>0$ such that $R^tv>0$
and set the $d \times d$ symmetric matrix $V$ such that  $V_{ij}=v_j v_i>0$. Notice that $V>0$ and
$VR>0$. Let $f(x) = x^tVx$. Proceeding as in the proof of Theorem \ref{thm:finite-boundary}, we
have
$$E_x[f(Z(t))]+  E_x[\int_0^t Z(s)^t V R \gamma ds] \geq  E_x[\sum_{i=1}^d \int_0^t Z(s)^t V R_{\cdot i} dY_i(s)]$$

Taking expectations with respect to $\pi$, and using Proposition \ref{prop:boundary} we have
$$\int_{F_i} \langle e, s \rangle \nu_i(ds) \leq \frac{\alpha}{\beta_i}\int_S \langle e, s \rangle \pi(ds) ,\, i=1,\hdots,d$$
where $\alpha \eqq \max_j\set{(V R \gamma)_j}>0$ and $\beta_i \eqq \min_j\set{(V R_{\cdot
i})_j}>0,\, i=1,\hdots,d$. We have implicitly used the finiteness of second moments for the
stationary distribution.

\subsection{Algorithm setting}

To apply the algorithm we still need to specify (i) $\{S_n\}$, the sequence of grid for which we
will approximate (\ref{eq:bar}), and (ii) $H_m$, the sequence of finite dimensional subspaces
approximating $C^2_b(S)$.

\subsubsection*{Grid choice} Harrison and Williams \cite{Harrison-Williams:87b} proved that the stationary distribution for a ``standard'' SRBM has a
separable density function (in the usual Cartesian coordinates) if and only if the covariance
matrix $\sigma$ satisfies the ``skew symmetric'' condition
\begin{equation}
2 \sigma_{ij} = R_{ij}+R_{ji} \quad \mbox{for} \quad i\neq j\label{eq:skew-symmetric}
\end{equation}
In this case, the marginal distribution for coordinate $i$ is exponential with rate $2\gamma_i$.
More over, the boundary measures are the restriction of the join distribution to the corresponding
faces of $S$. In practice the matrix $R$ does not fulfill this condition, but the closer it is to
fulfill it, the closer the marginal distributions should be to be exponentials.

Recently, Budhiraja and Lee \cite{Budhiraja-Lee:06} established the finiteness of the moment
generating function of the steady state distribution in a neighborhood of zero, proving the
exponential decay of it.

With this in mind we choose the grid to have an ``exponential spacing on the marginal'':
\begin{equation}
S_n = \set{ x\in S | x_i \in\set{[\log(n)-\log(j)]/\lambda_i \quad j=1,\hdots,n} \quad
i=1,\hdots,d}
\end{equation}
For the boundary, we just project this grind on the corresponding face
\begin{equation}
F^n_i = \set{ x\in S_n | x_i =0} \quad i=1,\hdots,n
\end{equation}
This is somehow related to the choice of a reference density on Harrison and Dai's algorithm
\cite{Dai-Harrison:92}.

\textbf{Remark.} In practice, the values for $\lambda$ were selected as follows: First we perform a
first run with the algorithm using low values for $\lambda$. Then these values were adjusted
according to the first moments on the marginal distributions.

\subsubsection*{Subspace choice}
For each $m\geq 1$, we choose $H_m = \set{f \mbox{ polynomial of degree }\leq m}$. The convergence
of $H_m\rightarrow C^2_b(S)$ is a well known result (Proper citation).

\subsection{Numerical Results}
In this section we will compare numerical results from our algorithm with some known analytical
results of particular instances of SRBM. Also, we will compare our algorithm with previous
algorithms using instances for which either simulation or good approximation results are known.

\subsubsection*{A Two-Dimensional SRBM}
Consider a two-dimensional SRBM 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 $\gamma<0$ reduces to $\mu_1<0$. Harrison \cite{Harrison:78}
computed a closed form solution for the stationary distribution density in polar coordinates:
\begin{equation}
p(x,y) = (2|\mu|)^{3/2}/(\pi^{1/2}) r^{-1/2}\exp(\mu_1r(1+\cos(\theta)))\cos(\theta/2)
\end{equation}
where $(x,y) =(r\cos(\theta),r\sin(\theta))$

Without loss of generality consider $\mu_1 =-1$. It can be shown (Greenberg \cite{Greenberg:86})
that
\begin{equation}
E_{\pi}[x] = 0.5 \quad \quad  E_{\pi}[y] = 0.75
\end{equation}
Taking $n=100$ we used our algorithm to get estimates for these moments. The results are shown on
Table \ref{table:two-dim}.
\begin{table}
  \centering
    \begin{tabular}{|c||c|c|c|c|c|c|c|c|} \hline
        $m$ & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \\\hline \hline
        $E_{\pi}[x]$ & 0.5 & 0.5 & 0.5 & 0.5 & 0.5 & 0.5 & 0.5 & 0.5 \\\hline
        $E_{\pi}[y]$ & 0.793 & 0.821 & 0.764 & 0.768 & 0.762 &0.758 & 0.760 & 0.750 \\\hline
    \end{tabular}

  \caption{Moment estimates Two-dimensional SRBM ($n=100$)}\label{table:two-dim}
\end{table}
We see that our algorithm provides good estimates even for low values of $m$. Each one of the
instances shown on Table \ref{table:two-dim} ran in less than $3$ minutes (on a regular Desktop PC)
on a MATLAB implementation of the algorithm. This type of SRBM is the only one (beside the skew
symmetric type) for which the actual distribution is know, so we can use this knowledge to test our
algorithm. Figure \ref{figure:two-dim} shows our marginal distribution estimates for $n=100$ and
$m=6$, compared with the actual marginal distributions (doted line).
\begin{figure}
\includegraphics[width=18cm]{figure1b.eps}
  \caption{Marginal distribution estimates Two-dimensional SRBM ($n=100$, $m=6$)}\label{figure:two-dim}
\end{figure}
It is proven (Harrison \cite{Harrison:78}) that the one dimensional marginal distribution for $x$
is exponential, however the same result does not hold for $y$. We see that our algorithm provides
good estimates for the marginal distributions.

\subsubsection*{Symmetric SRBM's}

A standard SRBM is said to be 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$. Positiveness of $\sigma$ implies $-1/(d-1)<\rho<1$ and the \emph{completely 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, Dai \cite{Dai:90} showed that
\begin{equation}
m_1=\hdots=m_d = \frac{1-(d-2)r+(d-1)r\rho}{2(r+1)}
\end{equation}
where $m_i = E_{\pi}[x_i]$ $i=1,\hdots,d$. We use our algorithm to compute estimates for these
moments for the case $d=2$. Imitating Dai's work, we let $\rho$ range through $\set{-0.9,
-0.5,0.0,0.5,0.9}$ and $r$ range through $\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$. We see that in each case the relative error is lower that $1\%$.

\begin{table}
  \centering
    \begin{tabular}{|c||c|c|c|c|c|} \hline
        $r/\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  \\\hline
    \end{tabular}

  \caption{Relative errors for Symmetric SRBM}\label{table:symmetric}
\end{table}


\subsubsection*{Skew-symmetric SRBM}
Harrison and Williams \cite{Harrison-Williams:87b} proved that a standard SRBM has a product form
stationary distribution if and only if $\gamma<0$ and condition \ref{eq:skew-symmetric} holds. In
this, case we know that the marginal distribution for $x_i$ is exponential with mean
$1/(2\gamma_i)$. Consider a two-dimensional SRBM associated to 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 $\gamma^t =(-1,-0.25)$. This
implies that $E_{\pi}[x_1] =0.5$ and  $E_{\pi}[x_2] =2$. Taking $n=100$ we used our algorithm to
get estimates for these moments. The results are shown on Table \ref{table:skew}.
\begin{table}
  \centering
    \begin{tabular}{|c||c|c|c|c|c|c|c|c|} \hline
        $m$ & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \\\hline \hline
        $E_{\pi}[x_1]$ & 0.4942 & 0.4943 & 0.4968 & 0.5014 & 0.5005 & 0.5011 & 0.4996 & 0.5002 \\\hline
        $E_{\pi}[x_2]$ & 2.0095 & 2.0094 & 2.0053 & 1.9976 & 1.9990 &1.9981 & 2.0005 & 1.9996 \\\hline
    \end{tabular}

  \caption{Moment estimates Skew-symmetric SRBM}\label{table:skew}
\end{table}
We see that our algorithm provides good estimates even for low values of $m$. Each one of the
instances shown on Table \ref{table:skew} ran in less than $3$ minutes (on a regular desktop PC) on
a MATLAB implementation of the algorithm. Since we know that each marginal distribution is
exponential, we can compare our estimated marginal distribution against the actual marginal
distribution for each coordinate. Figure \ref{figure:skew} shows our marginal distribution estimate
for $n=100$ and $m=6$, compared with the actual marginal distributions (doted line).
\begin{figure}
\includegraphics[width=18cm]{figure1b.eps}
  \caption{Marginal distribution estimates Skew-symmetric SRBM}\label{figure:skew}
\end{figure}
We see that our algorithm provides good estimates for the marginal distributions.

\subsubsection*{Suresh and Whitt's Experiments}

This section follows closely the analysis presented in Chapter 4 of Dai's dissertation
\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 associated with the data $(\sigma,R,\mu)$. When $d=2$
we have

$$ 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^2_2(C^2_{s_1}+C^2_{s_2})& -rho_1\rho_2C^2_{s_1}\\\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_i$ represents the traffic intensity at station
$i$.

Suresh and Whitt \cite{Suresh-Whitt:90} studied a system of two queues in tandem with the purpose
of finding the best order for the queues to minimize the average total waiting time. They consider
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 triples, namely, (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 them as Case 1 to Case 5, respectively. For
$C^2_{s_1}= C^2_{s_2}$ they consider two variability triples, namely, (1.0,0.5,0.5) and
(1.0,4.0,4.0). We will refer them as Case 6 and Case 7 respectively. When  $C^2_{s_1}\neq
C^2_{s_2}$, for each queue 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.

Extensive simulation experiments were conducted in order to obtain estimates for the expected
steady-state waiting times. When $C^2 = 0.5$, the $E_2$ distribution was used. When $C^2 = 1.0$,
the exponential distribution was used. When $C^2 > 1.0$, the $H_2$ distribution with balanced means
was used.

They also compare their simulation results with QNA approximations. Dai \cite{Dai:90} use this set
of experiments to test the QNET method. There, simulation and approximation results for the
expected waiting time at the second station are compared. We used these results to compare our
algorithm against the QNA and QNET methods. Table X-Y (appendix) give simulation estimates, QNET
estimates, QNA estimates, and our ``LP'' estimates. Table \ref{Suresh-Whitt-Table} summarizes all
balanced heavy traffic cases and gives an overall comparison of QNA estimates, QNET estimates and
LP estimates under heavy traffic. There, quantities shown represent the minimum between absolute
relative difference and absolute difference, with respect to simulation estimates.

\begin{table}
  \centering
    \begin{tabular}{||c|c|c||c|c|c||c|c|c||} \hline
      \multicolumn{3}{||c||}{} & \multicolumn{3}{|c||}{Order I} & \multicolumn{3}{|c||}{Order II}\\
      \hline
      Case & $\rho_1$ & $\rho_2$ & QNET-D & QNA-D & LP-D & QNET-D & QNA-D & LP-D\\ \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
      \multicolumn{3}{||c||}{Average} & 0.05 & 0.09 & 0.04 & 0.05 & 0.09 & 0.04\\ \hline
    \end{tabular}
  \caption{Overall Comparisons with QNA and QNET Approximations in Heavy Traffic}\label{Suresh-Whitt-Table}
\end{table}

We used $n=100$ and $m=6$ for all experiments. Each run took less than 3 minutes on regular Desktop
PC. As expected, QNET and LP methods give ``better'' estimates than the ones from QNA method, under
balanced heavy traffic conditions. Also, LP method performs as good as the QNET method, most of the
time better.

\section{Applications: Two Class Queue}\label{sec:app-twoclassqueue}

Consider a single queue system with $m$ classes, with class 1 denoting a guaranteed-rate or
high-priority class, and classes $2,\hdots,m$ being different best-effort classes that are labeled
according to their priority level i.e., class $i$ has higher priority than all classes $j>i$. The
guaranteed-rate class, assumed to arrive according to a  Poisson process with rate $\lambda_1$,
engage {\it one unit of capacity} each for i.i.d. exponentially distributed amounts of time with
rate $\mu_1$, provided  that the total number of guaranteed users connected is less than the system
capacity $C$; otherwise they are denied service.

The {\it best-effort users of class }$i$ ($i=2,\hdots,m$), assumed to arrive according to a Poisson
process with rate $\lambda_i$ (independent of everything else), are always admitted into the
system, requiring a total processing time exponentially distributed with rate $\mu_i$. When there
is enough capacity not used up by users from classes $j<i$, class $i$ users are allocated a nominal
processing rate corresponding to one unit of capacity, and when capacity is not sufficient, they
share the available capacity in an {\it egalitarian manner}:
$$\mbox{Class } i \mbox{ service rate at time }t = \frac{(C-\sum_{j<i}Q_j(t))^+}{Q_i(t)}\wedge1$$
where $Q_i(t)$ denotes the number of class $i$ users in the system.

We see that $Q_1$ evolves according to the state of an M/M/C/C system, while $Q_i\, (i=2,\hdots,m)$
evolve as a M/M/C(t), where  the available capacity, $C(t)$, is stochastically modulated by the
number of higher priority users present in the system.  Despite its simple structure, exact
analysis of this multiclass system is not straightforward and relies either on simulation or on
numerical methods that offer little insight as to its structural behavior.

Maglaras and Zeevi \cite{Maglaras-Zeevi:04} studied this system and derived a diffusion
approximation under the Halfin and Whitt \cite{Halfin-Whitt:81} asymptotic heavy traffic regime.
This regime is defined by letting capacity grow large and concurrently letting the system
utilization approach 1 at an appropriate rate. Specifically, Maglaras and Zeevi considered a
sequence of systems with capacity $C^n=n$ 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,\hdots,m$ such that
$\sum_i \kappa_i =1$. There $\kappa$ denotes the vector of {\it relative workload contributions}

\begin{equation}\label{eq:kappa}
\kappa_i = \lim_{n\rightarrow\infty}\frac{\lambda^n_i/\mu_i}{\sum_{j=1^m}\lambda^n_j/\mu_j}\quad
i=1,\hdots,m
\end{equation}

Defining the normalized state processes
$$ X_i^n(t) \equiv \frac{Q_i^n(t)-\kappa_i n}{\sqrt{n}}\quad i=1,\hdots,m$$
they proved the following theorem:

\begin{theorem}
Suppose that for some $\xi\,\in\,\Real^m$, $Q^n_i(0) = \lfloor n \kappa_i + \sqrt{n}\xi_i\rfloor$
for $i=1,\hdots,m$. Then, $X^n\rightarrow X$ in $D^m[0,\infty)$ as $n\rightarrow\infty$, where $X$
is a diffusion process. Specifically, $X$ is the unique strong solution of the following stochastic
differential equation:

\begin{equation}\label{eq:SDE-two-class}
dX(t) = b(X(t)dt+\Sigma dW(t), \quad \quad X(0) =\xi,
\end{equation}
where $W =(W(t):t\geq 0)$ is a standard Brownian motion in $\Real^m$, the infinitesimal drift
function $b_i(\cdot)$ for the $i$th component is
\begin{eqnarray*}
b_i(x) & = & -\mu_i \gamma_i - \mu_i x_i \quad \quad i=1,\hdots,m-1\\
b_m(x) & = & \left\{\begin{array}{lr} -\mu_m \gamma_m -\mu_m x_m & \sum_{i=1}^m x_i\leq 0\\
-mu_m\gamma_m + \sum_{i=1}^{m-1}x_i & \sum_{i=1}^{m}x_i >0\\ \end{array} \right.\\
\end{eqnarray*}
and $\Sigma \equiv $ diag$(\sigma_1, \sigma_2,\hdots,\sigma_m$), with $\sigma_i^2 =
2\mu_i\kappa_i$.
\end{theorem}
The first $m-1$ components of $X$ are simple O-U processes and clearly admits a steady-state,
however the last component has more complicated structure, and thus is not clear {\it a priori}
under what conditions a steady-state will exist for $X$. Maglaras and Zeevi proved $X$ is positive
recurrent and that it admits a unique stationary distribution if and only if $\gamma =
\sum_{i=1}^m\gamma_i >0$, which is also the steady-state distribution.

Our objective in this section is to apply our algorithm to approximate the stationary distribution
of these type of processes. For this we need first to check the sufficiency and necessity of the
BAR condition, and second, to provide tightness bounds for the steady-state distribution.

\subsection{Bar Condition}
Necessity of the BAR condition comes from the following proposition.
\begin{proposition}
Suppose $\pi$ is the stationary distribution for $X$. Then for each $f\in C^2_b(\Real^m)$,
\begin{equation}
\int_{\Real^m} \left(\sum_{i=1}^m 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 \label{eq:bar-twoserver}
\end{equation}
\end{proposition}
\proof First we need to prove that $P\left[\int_0^t |b_i(x)|dx <\infty \quad \forall \, t\geq 0
\right]=1$. If this is the case, then we can apply Ito's lemma, and take expectation $E_x$ to
obtain
\begin{equation*}
E_x[f(X(t))] - f(x) = E_x\left[\int_0^t \left(\sum_{i=1}^m b_i(X(s)) \frac{\partial f(x)}{\partial
x_i} + \frac{\sigma^2_i}{2} \frac{\partial^2 f(x)}{\partial x_i^2} \right)ds)\right] +
E_x\left[\int_0^t \crs{ \nabla f(X(s))\Sigma, dW(s)}\right]
\end{equation*}
Since $f\in C^2_b(\Real^m)$ the last term is a martingale, so its expectation is 0. Integrating
both sides with respect to the stationary distribution $\pi$, we obtain (\ref{eq:bar-twoserver})
when applying Fubini's lemma.\qed

Sufficiency of the BAR condition comes from the following proposition (see appendix for a proof).
\begin{proposition}\label{prop:bar-sufficiency}
Suppose that $\pi$ is a probability measure on $\Real^m$. If $\pi$ satisfies
(\ref{eq:bar-twoserver}), then  $\pi$ is  the stationary distribution for $X$.
\end{proposition}

\subsection{Tightness condition}
Since $X$ lives on $\Real^m$, this application does not involve any boundary measure associated to
the stationary distribution. This leaves us with just one ''tightness'' constraint, the one for the
''interior'' distribution. In exchange the simple sum of first moments of each coordinate under the
stationary distribution does not imply tightness.

For $i=1,\hdots,m-1$ consider $f_i(x) = \mu_1^{-1}(x_i-\gamma_i)$. Applying Ito, taking expectation
w.r.t to $\pi$, dividing by $t$, and taking $t\rightarrow \infty$ we have that $E_{\pi}[x_i^2] =
\gamma_i^2 + \kappa_i$ ($i=1,\hdots,m-1$). Taking $f_i(x) = x_i$ and repeating the same reasoning
we have that $E_{\pi}[x_i] = -\gamma_i$ ($i=1,\hdots,m-1$). Therefore we can conclude that
$X_i(\infty) \sim N(-\gamma_i,\kappa_i)$. In order to state the tightness constraint, we still need
to analyze the last component of $X$. Consider:
$$f(x)= \sum_{i<m} (1/2) C_{1,i} x_i^2 + \sum_{i<m} (1/2) C_{2,i} x_i  + C_3\sqrt{1+x_m^2} +C_4$$
where $C_4$ is set so that $f$ is nonnegative, $c>0$ is a constant that ensures that
$\mathcal{A}f(x)<-\eps$ for $x\in K^c$, where $K$ is a compact set, and
$$C_3=\frac{c+ \sum_{i<m}\gamma^2_i}{\gamma \mu_m}\quad\quad C_{2,i}=(C_3\mu_m -\gamma_i)/\mu_i\quad\quad
C_{1,i}=1/\mu_i$$

Maglaras and Zeevi \cite{Maglaras-Zeevi:04} shown that this function is a valid Lyapunov function
for $X$.  This function can be used, for example, to establish the finiteness of the moment
generating function of the steady state distribution. In particular it can be shown (see Glynn and
Zeevi \cite{Glynn-Zeevi:07}) that
$$E[\sum_{i=1}^m |x_i|]\leq \max_{x \in K} \set{|\mathcal{A}f(x)|}$$
which the tightness constraint that we will use.

\subsection{Algorithm setting}

To apply the algorithm we still need to specify (i) $\{S_n\}$, the sequence of grid for which we
will approximate (\ref{eq:bar}), and (ii) $H_m$, the sequence of finite dimensional subspaces
approximating $C^2_b(S)$.

\subsubsection*{Grid choice}

We know that $X_i(\infty)\sim N(-\gamma_i,\kappa_i) (i=1,\hdots,m-1)$, but we don't have any prior
information about the distribution $X_m(\infty)$, but that it has an exponential decay. With this
in mind we choose the grid to have an ``Normal'' spacing on the marginal distribution for $X_i$
$i=1,\hdots,m-1$, and an ''Uniform'' spacing on the marginal distribution for $X_m$:
\begin{eqnarray}
S_n & = & \{ x\in S |\{ x_i = (\Phi^{-1}(1/2+j/(2n))+\gamma_i)/\sqrt{\kappa_i} \quad
j=0,\pm1\hdots,\pm n-1\\
& & i=1,\hdots,m-1\} \times \{x_m =\alpha j/n \quad j=-n,\hdots,n\}\}
\end{eqnarray}

\subsubsection*{Subspace choice}
For each $m\geq 1$, we choose $H_m = \set{f \mbox{ polynomial of degree }\leq m}$.

\subsection{Numerical Results}

The following example was studied by Maglaras and Zeevi \cite{Maglaras-Zeevi:04}, and restrict
attention to a two-class system with a single guaranteed and best-effort class. Consider a system
with capacity $C= 100$, $\mu_1 =1$, $\mu_2 =2$, $\lambda_1=47.5$ and $\lambda_2=95$ (here
$\rho=0.95)$. The idea is to approximate this system using the two dimensional diffusion described
in this section. For that purpose we can solve for $\kappa_i$ and $\gamma_i$ $(i=1,\hdots,m)$ using
(\ref{eq:kappa}). In this example:

$$\kappa_i = \frac{\lambda_i/\mu_i}{\lambda_1/\mu_1 + \lambda_2/\mu_2} \Rightarrow \kappa_1 =
47.5/(47.5+95/2)=0.5 \quad \mbox{and}\quad \kappa_1 =1-\kappa_2 = 0.5$$ and
$$\gamma_i =\frac{\kappa_i\mu_iC-\lambda_i}{\mu_i\sqrt{C}}\Rightarrow \gamma_1 =\frac{0.5\cdot 1\cdot
100-47.5}{\sqrt{100}}=0.25 \quad \mbox{and}\quad \gamma_2 =\frac{0.5\cdot 2\cdot
100-95}{2\cdot\sqrt{100}}=0.25$$

Queue length steady state distribution can be approximated as follows
$$Q^C_i(\infty) = \kappa_i C +\sqrt{C}X_i(\infty)\quad i=1,\hdots,m$$

Figure \ref{figure:two-class} compares the approximated marginal distributions resulting from the
algorithm (in what follows we will use $m=100$, $n=4$) with the actual marginal distribution
computed via Monte Carlo simulation (doted line).

\begin{figure}
\centering
\includegraphics[width=13cm]{figure3a.eps}
  \caption{Marginal distribution estimates Two Class Queue ($n=100$, $m=6$)}\label{figure:two-class}
\end{figure}

This does not look like good approximation to the marginal distributions. However, when comparing
moments for the second coordinate, we have that the approximation does a good job. See table \ref{table:two-class}.
\begin{table}
  \centering
    \begin{tabular}{||c||c|c|c|c||}
      \hline
       & $E[x^1_2]$ & $E[x^2_2]$ & $E[x^3_2]$ & $E[x^4_2]$ \\\hline \hline
      Simulation     & 1.19 & 8.31 & 70.67 & 854 \\
      Algorithm      & 1.22 & 9.14 & 68.2 & 592 \\
      Relative error & 0.03 & 0.10 & 0.03 & 0.30 \\
      \hline
    \end{tabular}
  \caption{Moment approximation for second coordinate Two Class Queue}\label{table:two-class}
\end{table}

We would expect for the actual steady state distribution to have ''smooth'' marginal distributions,
but this condition is not explicitly imposed on the LP formulation. It turns out that the
approximated BAR condition can not rule out such erratic behavior on the marginal distributions, at
least not for such low values of $m$ and $n$. Also, we would expect for a discretization of the
stationary distribution to be both feasible and to achieve an objective vale close to the optimal
one, and to have an ''smooth'' behavior. With this in mind we decided to impose ''smoothness'' on
the algorithm's output. Consider a neighborhood parameter $r>0$, and a ''smoothness'' parameter
$\delta$. We will impose that the probability assigned to two points on $S^n$ whose distance
(euclidian norm) from each other is less or equal to $r$ should not be greater in absolute value to
$\delta$

\begin{equation}
\abs{\pi^n_m(x) - \pi^n_m(y)}\leq \delta \quad \forall\, (x,y)\in S^n \quad s.t. \, ||x-y||<r
\end{equation}

\textbf{Remark} The choice of $r$ and $\delta$ will be case dependent. The greater the value of $r$
is, the greater is the number of constraints added. This will affect the complexity of the LP. On
the other hand, low values for $\delta$ will reduce the space of feasible solutions, and at some
point it could even discard a true steady-state distribution discretization.

In practice we introduce these constraints on each marginal distribution separately. We choose $r$
in order to include just one neighbor in each marginal, and $\delta$ low enough to rule out peaks
as in figure \ref{figure:two-class}. For $\delta_1 =0.0026$ ( smoothness parameter for the first
coordinate marginal) and  $\delta_2 =0.0011$ ( smoothness parameter for the first coordinate
marginal) we obtain the moment approximation shown on Table \ref{table:two-class-2}. Figure
\ref{figure:two-class-2} shows the marginal obtained including the smoothing constraints.

\begin{figure}
\centering
\includegraphics[width=13cm]{figure3b.eps}
  \caption{Smooth marginal distribution estimates Two Class Queue ($n=100$, $m=6$)}\label{figure:two-class-2}
\end{figure}

\begin{table}
  \centering
    \begin{tabular}{||c||c|c|c|c||}
      \hline
       & $E[x^1_2]$ & $E[x^2_2]$ & $E[x^3_2]$ & $E[x^4_2]$ \\\hline \hline
      Simulation     & 1.19 & 8.31 & 70.67 & 854 \\
      Algorithm      & 1.22 & 9.14 & 68.2 & 592 \\
      Smooth Alg.    & 1.23 & 8.92 & 68.7 & 618 \\\hline\hline
      Relative error & 0.03 & 0.07 & 0.03 & 0.28 \\
      \hline
    \end{tabular}
  \caption{Smooth moment approximation for second coordinate Two Class Queue}\label{table:two-class-2}
\end{table}

\section{Robust Steady State Estimates (?)}\label{sec:robust}

\newpage
\appendix
\section{Proofs}\label{sec:proofs}
\subsection{Proof of Theorem \ref{thm:finite-boundary}}
Define $\gamma = - R^{-1}\mu > 0$. Now, take a d-dimensional vector $v>0$ such that $R^tv>0$ (the
\emph{completely S} is close under the transpose and inverse operations). Fix $a>0$ and let $f(x) =
\exp(a x^tv)$. Now let $f_n \in C^2_b(S)$ be a positive function that agrees with $f$ for $x\in
B_S(n)\eqq \set{x\in S, \abs{x}\leq n}$. Applying It\^{o}'s formula, we have that
$$f_n(Z(t))-f_n(z(0)) -\int_0^t L f_n(Z(s))ds - \sum_{i=1}^d \int_0^t D_i f_n(Z(s)) dY_i(s)$$
is a $P_x$ martingale. Let $T_n \eqq \inf\set{t\geq 0, \abs{Z(t)}\geq n}$. Since $f_n>0$ we have
$$E_x[f(Z(t\wedge T_n))]\geq E_x[\int_0^{t\wedge T_N} L f(Z(s))ds] + E_x[\sum_{i=1}^d \int_0^{t\wedge T_N} D_i f(Z(s)) dY_i(s)]$$
For this particular function we have $L f(x) = f(x)\left[-a v^t R \cdot \gamma +
\frac{a^2}{2}\sum_{i=1}^d \sum_{j=1}^d \Gamma_{ij} v_i v_j \right]$ and $D_i f(x) = a f(x) [v^t
R_{\cdot i}]$. Taking $a< \frac{2 v^t R \cdot \gamma}{\sum_{i=1}^d \sum_{j=1}^d \Gamma_{ij} v_i
v_j}$ we have

$$E_x[f(Z(t\wedge T_n))]+ \alpha E_x[\int_0^{t\wedge T_N} f(Z(s))ds] \geq  E_x[\sum_{i=1}^d \beta_i \int_0^{t\wedge T_N} f(Z(s)) dY_i(s)]$$
where $\alpha \eqq a v^t R \cdot \gamma - \frac{a^2}{2}\sum_{i=1}^d \sum_{j=1}^d \Gamma_{ij} v_i
v_j >0$ and $\beta_i\eqq v^tR_{\cdot i}>0\quad i=1,\hdots,d$. Since $f>0$ and taking
$\set{f_n}_{n\in \mathbb{N}}$ increasing, the Monotone convergence theorem implies the following:

$$E_x[f(Z(t))]+ \alpha E_x[\int_0^t f(Z(s))ds] \geq  E_x[\sum_{i=1}^d \beta_i \int_0^t f(Z(s)) dY_i(s)]$$

Taking expectations with respect to $\pi$ (and considering $a>0$ small enough for $E_{\pi}[e^{a
v^tx}] $ to exist) we have that

$$E_{\pi}[f]+ \alpha t E_{\pi}[f] \geq  t \sum_{i=1}^d \beta_i \int_{F_i} f \nu_i(ds)$$
where we have used proposition \ref{prop:boundary}, the definition of the stationary distribution
and Tonelli's theorem. Dividing by $t$ and taking $t\rightarrow \infty$ we conclude the following
$$ \int_{F_i} e^{a v^tx} \nu_i(ds) \leq \frac{\alpha}{\beta_i} E_{\pi}[e^{av^tx}]<\infty , \,i=1\hdots,d$$

\subsection{Proof of Proposition \ref{prop:bar-sufficiency}}
First we will need some side results. For $f\, \in \, C^2(\Real^m)$ define the operator
$\mathcal{A}$ as follows
$$\mathcal{A} f \equiv \sum_{i=1}^m b_i(x) \frac{\partial f(x)}{\partial x_i} + \frac{\sigma^2_i}{2}
\frac{\partial^2 f(x)}{\partial x_i^2}$$
\begin{lemma}\label{lem:maximum}
The operator $(\mathcal{A},C^2_b(\Real^m)$ satisfies the positive maximum principle, i.e., whenever
$f\,\in C^2_b(\Real^m)$, $x_0\in \Real^m$, and $\sup_{x \in \Real^m}f(x) = f(x_0)\geq 0$, we have
$\mathcal{A}f(x_0) \leq 0$.
\end{lemma}

\proof Consider $f\, \in C^2_b(\Real^m)$. Applying Ito
\begin{equation}
f(X(t)) - \int_0^t \mathcal{A}f(X(s))ds =  f(x_0) + \int_0^t \crs{ \nabla f(X(s))\Sigma, dW(s)}
\label{eq:ito-twoserver}
\end{equation}

Now take expectations w.r.t. $P_{x_0}$, we have
\begin{equation}
E_{x_0}[f(X(t))] - f(x_0) = E_{x_0}\left[\int_0^t \mathcal{A}f(X(s))ds\right]
\end{equation}
Dividing by $t$ and taking $t\rightarrow 0$ by the continuity of $\mathcal{A}f$ and the continuity
of $X$, we get $\mathcal{A}f(x_0)\leq 0$, since the left hand side above is non-positive.\qed

\begin{proposition}\label{prop:mart-prob}
For any probability measure $\pi$ on $\Real^m$, the measure $P_{\pi} \equiv \int_{\Real^m} P_{x}
\pi(dx)$ is a solution for the martingale problem for $(\mathcal{A},\pi)$.
\end{proposition}
\proof It is enough to show that for $f\in C^2_b(\Real^m)$ and each $x \in \Real^m$,
\begin{equation}
f(X(t)) - \int_0^t\mathcal{A}f(X(s))ds
\end{equation}
is a $P_x$-martingale. However, this is direct a consequence of (\ref{eq:ito-twoserver}).\qed

\noindent \textbf{Proof of Proposition \ref{prop:bar-sufficiency}.} Clearly $C^2_b(\Real^m)$ is an
algebra, and also it is dense in $C(\Real^m$) (\textbf{CITE}). By Lemma \ref{lem:maximum}, the
operator ($\mathcal{A},C^2_b(\Real^m$) 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 distribution for a solution of the martingale problem for
($\mathcal{A},\pi$). By proposition \ref{prop:mart-prob} we know that $P_{\pi}$ is a solution of
the martingale problem for ($\mathcal{A},\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 results follows from the uniqueness of
the stationary distribution.

\section{Numerical Results}\label{sec:numerical}

 7 Pages of numerical results for the SRBM?????

%\input{Appendix}

\noi\textbf{Acknowledgements.} We are indebted to...

%===================================================================
\begin{thebibliography}{99}
\bibitem{Hernandez-Lasserre} O. Hernandez-Lerma and J. B. Lasserre. Markov Chains and Invariant
Probabilities. {\it Progress in Mathematics}, Vol 211, Birkh$\ddot{a}$user.
\bibitem{Gamarnik-Zeevi:06} D. Gamarnik and A. Zeevi. Validity of heavy traffic steady-state
approximations in generalized Jackson networks, in {\it Annals of Applied Probability}, Vol 16, No.
1, 2006.
\bibitem{Gurvich-Zeevi} I. Gurvich and A. Zeevi. Validity of Heavy-Traffic Steady-State
Approximations in Open Queueing Networks: Sufficient conditions involving state-space collapse,
forthcoming.
\bibitem{Harrison-Reiman:81} J. M. Harrison and M. Reiman. Reflected Brownian Motion on an Orthan,
in {\it The Annals of Probability}, Vol. 9, No. 2, 1981.
\bibitem{Reiman:84} M. I. Reiman. Open Queueing Networks in Heavy Traffic, in {\it Mathematics of Operations
Research}, Vol. 9, No. 3, 1984.
\bibitem{Dai:90} J. G. Dai. Steady-state analysis of reflected Brownian motions: Characterization,
numerical methods and queueing applications. Ph.D. dissertation, Dept. Math. Stanford University,
1990.
\bibitem{Dai-Harrison:92} J. G. Dai and J. M. Harrison. Reflected Brownian motion in an orthant:
Numerical methods for steady-state analysis, in {\it The Annals of Applied probability}, Vol 2, No.
1, 1992.
\bibitem{Dai-Kurtz:Unpublished} J.G. Dai and T. G. Kurtz. The sufficiency of the basic adjoint
relationship. Unpublished manuscript.
\bibitem{Harrison-Williams:87a} J. M. Harrison and R. J. Williams.  Brownian models of open
queueing networks with homogeneous customer populations, in {\it Stochastics}, 1987.
\bibitem{Harrison-Williams:87b} J. M. Harrison and R. J. Williams. Multidimensional reflected
Brownian motions having exponential stationary distributions, in {\it Annals of Applied
Probability}, 1987.
\bibitem{Harrison-Nguyen:90} J. M. Harrison and V. Nguyen. The QNET method for two-moment analysis
of open queueing networks, in {\it Queueing Systems: Theory and Applications}, 1990.
\bibitem{Taylor-Williams:93} L.M. Taylor and R. J. Williams. Existence and uniqueness of
semimartingale reflecting Brownian motion in an orthant, in {\it Probability Theory and Related
Fields}, 1993.
\bibitem{Greenberg:86} B. S. Greenberg. Queueing systems with returning customers and the order of
tandem queues. Ph.D. dissertation, University of California, Berkeley, 1986.
\bibitem{Harrison:78} J. M. Harrison. The diffusion approximation for tandem queues in heavy
traffic, in {\it Advances in Applied Probability}, 1978.
\bibitem{Budhiraja-Lee:06} A. Budhiraja and C. Lee. Long time asymptotics for constrained
diffusions in polyhedral domains, in {\it Stochastic Processes and their
Applications},doi:10.1016/j.spa.2006.11.007.
\bibitem{Reiman-Williams:88} M. I. Reiman and R. J. Williams.  A boundary property of
semimartingale reflecting  Brownian motions, in {\it Probability Theory and Related Fields}, 1988.
\bibitem{Dupuis-Williams:94} P. Dupuis and R. J. Williams. Lyapunov functions for semimartingale Brownian
motions, in {\it Annals of Probability}, 1994.
\bibitem{Dai-Harrison:07} J. G. Dai and J. M. Harrison. Reflecting Brownian motion in the orthan:
An illuminating example of inestability, forthcoming.
\bibitem{Glynn-Zeevi:07} P. W. Glynn and A. Zeevi. Bounding Stationary Expectations of Markov
Procesess, forthcomming.
\bibitem{Suresh-Whitt:90} S. Suresh and W. Whitt. Arranging queues in series: a simulation
experiment, in {\it Management Sciences}, Vol. 9, No. 3, 1990.
\bibitem{Maglaras-Zeevi:04} C. Maglaras and A. Zeevi. Diffusion Approximations for a Multiclass
Markovian Service System, in {\it Mathematics of Operations Research}, Vol. 29, No. 4, 2004.
\bibitem{Halfin-Whitt:81} S. Halfin and W. Whitt. Heavy-traffic limits for queues with many
exponential servers, in {\it Operations Research}, Vol 29, 1981.
\bibitem{Ehtier-Kurtz} S. N. Either and T. G. Kurtz. Markov Processes, Characterization and
Convergence. {\it Wiley Series in Probability and Mathematical Statistics}.
\bibitem{Shen-Chen:02} X. Shen, H. Chen, J. Dai and W. Dai. The Finite Element method for Computing
the Stationary Distribution of an SRBM in a Hypercube with Applications to Finite Buffer Queueing
Networds, in {\it Queueing Systems: Theory and Applications}, Vol 42, 2002
\bibitem{Shen-Chen:03} X. Shen, H. Chen. Computing the stationary distribution of an SRBM in an orthant with applications to queueing
networks, in {\it Queueing Systems}, Vol 45, 2003
\bibitem{Manne:60} A. Manne. Linear programming and sequencial desicions, in {\it Management Science}, Vol.
6, 1960.
\bibitem{Mendiondo-stockbrige:98} M. Mendiondo, R. Stockbridge. Approximation of Infinite-Dimensional Linear Programming Problems
which Arise in Stochastic Control, in {\it Journal of  Control Optimization}, Vol. 36, 1998.
\bibitem{Betsimas-etal:01} D. Bertsimas, J. Tsitsiklis and D. Gamarnik. Performance of Multiclass Markovian Queueing Networks via Piecewise
Linear LYapunov Functions, in {\it Annals of Applied Probability,} Vol. 11, 2001.
\bibitem{Kumar-Kumar:94} S. Kumar and P. Kumar. Performance bounds for queueing networks and scheduling policies
, in {\it IEEE Trans. on Automatic Control,} Vol. 39, 1994.
\bibitem{Sigman-Yao:93} K. Sigman and D. Yao. Finite Moments for Inventory Processes, in {\it Annals of Applied
Probability}, Vol. 3, 1993.
\bibitem{Reiman-Mandelbaum-Massey:98} M.Reiman, A. Mandelbaum and W. Massey. Strong Approximations for Markovian Service Networks, in {\it Queueing
Systems}, Vol. 30, 1998.
\end{thebibliography}
%===================================================================

\end{document}


Old Theorem 2 proof!

\begin{proof} We know there exist sequences $\set{\phi_n}$ ($\set{\phi_{jn}}$)of distribution (finite)
measures on $S$ such that
\begin{itemize}
\item[(a)] For every $n\geq 1$, $\phi_n$ ($\phi_{jn}$) has finite support $S^n$ ($F_j^n$),
that is, $\phi_n$ ($\phi_{jn}$) is of the form $\phi_n=\sum_{x_i\in S^n} \beta^n_i\delta_{x_i}$, ($\phi_{jn}=\sum_{x_i\in F_j^n}
\beta^n_{ji}\delta_{x_i}$), with $\beta^n_i\geq 0$ $(\beta^n_{ji}\geq 0)$, and $\sum_{x_i\in S^n}^n \beta^n_i =1$.
\item[(b)] The sequences $\set{\phi_n}$ ($\set{\phi_{jn}}$)converges weakly to $\pi$ ($\nu_j$).
\end{itemize}

From the definition of weak convergence (plus a truncation argument) we know that $E_{\phi_n}[g(x)] \rightarrow E_{\pi}[g(x)]$
($E_{\phi_{jn}}[g_j(x)] \rightarrow E_{\nu_j}[g_j(x)]$), therefore there exists a $n^*$ for which $(\beta^n, \beta^n_j)$ is
$\mathcal{P}_{nm}$-feasible $\forall \,m,\,\forall\,n\geq 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}$, $\phi_j\in\mathcal{M}_j \, (j=1\hdots,d)$ and $f\in C^2_b(S)$ let $B:\mathcal{M}\times\Pi_{j=1}^d\mathcal{M}_j\rightarrow
\Real$ denote the following operator:
$$ B(\phi,\phi_j,f) =\int_S \left(Gf\right)(x)\, \phi(dx)  + \sum_{i=1}^d \int_{F_i}
\left(D_if\right)(x)\, \phi_i(dx).$$

Fix $m>0$, and consider $(\lambda^n_{m},\, \gamma^n_{jm}) $ an optimal solution to $\mathcal{P}_{nm}$. For $h\in\hat{H}_m$ we have that

$$ \lim_{n\uparrow\infty} \abs{B(\pi(\lambda^n_m),\nu_j(\gamma^n_{jm}),h)}\leq \lim_{n\uparrow\infty}
\abs{B(\pi(\beta^n_m),\nu_j(\beta^n_{jm}),h)} = B(\pi,\nu_j,h) = 0.$$

The first inequality comes from the optimality of $(\lambda^n_{m},\, \gamma^n_{jm})$ and the feasibility of $(\beta^n, \beta^n_j)$ for $n>n^*$, while
the second inequality comes from the weak convergence result for $\set{\phi_n}$ ($\set{\phi_{jn}}$). Due to the tightness condition, there exists a
subsequence $\set{n^m_k}$ of integers such that $(\pi(\lambda^{n^m_k}_m), \nu_j(\gamma_j^{n^m_k}))\Rightarrow (\pi_m,\nu_{jm})$, proper probability
and finite measures respectively. We also conclude that $B(\pi_m,\nu_m,h) = 0 \,\forall\,h\in H_m$. Now consider the sequences $\set{\pi_m}$ and
$\set{\nu_{jm}}$. These sequences are tight, therefore $(\pi_m, \nu_{jm})\Rightarrow (\hat{\pi},\hat{\nu_j})$ along a further subsequence.

Fix $f\in C^2_b(S)$ and assume that $(Gf)(x)<C$ for all $x\in S$ and $(D_jf)(x)<C_j(j=1,\hdots,d)$ for all $x\in F_i$, for some finite constants $C$
and $C_j (j=1,\hdots,d)$. Now consider $A\subset S$ and $A_j\subset F_i$ such that $\pi_m(A^c)\leq \epsilon/C$, and $\nu_{jm}(A_j^c)\leq \epsilon/C_j
\,(j=1,\hdots,d)$ for all $m>0$. There exists a $m^*$ such that exists a $g_m \in H^m$ such that $\sup_A | (G(g_m-f))(x) |\leq \epsilon $,
$\sup_{A_j} | (D_j(g_m-f))(x) |\leq \epsilon \quad \forall j=\set{1,\hdots,d}$, $(Gg_m)(x)<C$ for all $x\in S$ and $(D_jg_m)(x)<C_j(j=1,\hdots,d)$
for all $x\in F_i$ for all $m>m^*$. We have that for $m>m^*$,
\begin{eqnarray*}
B(\pi_m,\nu_{jm},f) & = & B(\pi_m,\nu_{jm},f-g_m)\\
& \leq &  \int_A (G (f-g))(x)\pi_m(dx) + \sum_{j=1}^d \int_{A_j} (D_j (f-g_j))(x) \nu_{mj}(dx) + 4d\epsilon\\
& \leq & \epsilon + \epsilon \cdot \sum_{i=1}^d K_j + 4 d \epsilon.
\end{eqnarray*}
Since $\epsilon$ is arbitrary  we conclude that, for all $f\in C^2_b(S)$, $ B(\hat{\pi},\hat{\nu_{j}},f) = \lim_{m\uparrow \infty}
B(\pi_m,\nu_{jm},f) = 0.$ The result comes from the uniqueness of the stationary distribution, and a diagonal argument. $\Box$
\end{proof}
