
Consider a positive recurrent continuous-time Markov chain (CTMC) with infinite state space $S$ and
infinitesimal generator matrix $Q$. A vector $\pi$ is a stationary distribution for this CTMC if
and only if $\pi^t Q=0$ and $e^t\pi=1$, where $e$ (abusing notation) is a column vector with all
entries being 1. One way of finding $\pi$ is to solve the following infinite-dimensional
\emph{linear program} (LP),

\begin{equation}\label{eq:LP-CTMC}
\set{\min u\,\,| \pi^tQ \leq e^t u \,;\, -\pi^tQ \leq e^t u \,;\, e^t\pi = 1 \,;\, \pi\geq 0}
\end{equation}
Solving (\ref{eq:LP-CTMC}) exactly is virtually impossible. One possible approach is trying to ''approximate'' $\pi$ by solving a
finite-dimensional LP that approximate in some sense the original LP.

When $S$ is a locally compact separable metric space, the Banach space $\mathcal{C}_0(S)$ contains a {\it countable} dense subspace (is
{\it separable})
$$ H \equiv \set{h_1,h_2,\hdots}\subset \mathcal{C}_0(S)$$
By the denseness of $H$ in $\mathcal{C}_0(S)$, for any two probability measures $\mu$, $\nu$ we
have
\begin{eqnarray*}
\mu =\nu &\iff & \crs{\mu,h} =\crs{\nu,h} \quad \forall \,h\in \mathcal{C}_0(S)\\
&\iff & \crs{\mu,h} =\crs{\nu,h} \quad \forall \,h\in H\\
\end{eqnarray*}
Hence, solving (\ref{eq:LP-CTMC}) is equivalent to solving
\begin{equation}\label{eq:LP-CTMC-2}
\set{\min u\,\,| \crs{\pi^tQ,h} \leq e^t u \,\forall\,h\in H \,;\, -\crs{\pi^tQ,h} \leq e^t u
\,\forall\,h\in H \,;\, e^t\pi = 1 \,;\, \pi\geq 0}
\end{equation}
Notice that (\ref{eq:LP-CTMC-2}) has a {\it countable} set of constraints. The condition
$\crs{\pi^tQ,h}=0 \forall \,h\in\mathcal{C}_0(S) ;\quad e^t\pi=1$ corresponds to the {\it basic
adjoint relationship} (BAR), introduced by Harrison and Williams \cite{Harrison-Williams:87a} to
characterize the stationary distribution of a reflected Brownian motion in the nonnegative orthant.

Under the same assumptions, there exists a countable subset $\hat{S}$ dense in $S$.
$$\hat{S} \equiv \set{x_1,x_2,\hdots}$$
The key feature here is that, for any probability measure $\pi$ on $S$, there exist a sequence of
measures with finite support, that converges weakly to $\pi$. We can try to approximate $\pi$ by
solving a sequence of finite dimensional LPs.

\begin{equation}\label{eq:LP-CTMC-3}
\Prob_{nm}\equiv \set{\min u \,| \quad\left|\sum_{i=1}^n \lambda_i [ Q(h(x_i))]\right|\leq u \,
\forall \, h\,\in H_m \,;\,e^t\lambda=1\,;\,\lambda\geq0}
\end{equation}
where $H_m = \set{h_1,\hdots,h_m}$. The intuition here is that, as $n$ grows large we are able to
give better a approximation of an arbitrary probability distribution in $S$. On the other hand, as
$k$ grows large, we are giving a better approximation to the BAR condition. Hernandez and Lasserre
\cite{Hernandez-Lasserre} proved the validity of the approach when approximating the moment of an
inf-compact nonnegative function under the stationary distribution. The approach remains valid when
estimating $\pi$, provided that a tightness condition is imposed on each LP.

In what follows, we will adapt this idea to approximate stationary distributions for diffusions for
which the equivalent BAR condition is known to be necessary and sufficient, and will prove the
validity of the approach under certain assumptions.

Suppose the k-dimensional process $\set{X_t\in S \subset \Real^k :t\geq 0}$ arises as a
heavy-traffic limit for a stochastic process describing some queueing system. In generality $X_t$
could be decomposed as the sum of a ``free'' time-homogeneous diffusion process plus some
finite-variational process, reflecting the behavior of $X$ on the boundary of $S$, $\partial S$. We
will assume that $X_t$ is a diffusion process, ignoring the boundary behavior, although our results
are easily extended to the general case, as we will illustrate in section \ref{sec:app-SRBM}, when
dealing with the semimartingale reflected Brownian motion case.

Assume $\set{X_t\in S \subset \Real^k :t\geq 0}$ solves the following stochastic differential
equation
\begin{equation}\label{eq:SDE}
dX_t = b(X_t)dt + \sigma(x_t)dB_t \quad t\geq 0, \quad X_0=x
\end{equation}
where $b(x) \in \Real^k$, $\sigma(x) \in \Real^{k\times l}$ are continuous functions , and  $B_t$
is a l-dimensional Brownian motion. The infinitesimal generator $\mathcal{A}$ of $X_t$ is defined
by
\begin{equation}
\mathcal{A}f(x) = \lim_{t\downarrow 0} \frac{E^x[f(X_t)]-f(x)}{t} \quad x\in S
\end{equation}
The set of functions $f:S\rightarrow \Real$ such that the limit exits at $x$ is denoted by
$\mathcal{D}_{\mathcal{A}}(x)$, while $\mathcal{D}_{\mathcal{A}}$ denotes the set of functions for
which the limit exists for all $x\in S$. We will assume there exits a further set of functions $H
\subset C^2(S)$ such that $f\in \mathcal{D}_{\mathcal{A}}$ and

\begin{equation}
\mathcal{A}f(x) = \sum_{i}b_i(x)\frac{\partial f}{\partial x_i} +
\frac{1}{2}\sum_{i,j}(\sigma\sigma^t)_{i,j}(x) \frac{\partial^2 f}{\partial x_i \partial x_j}
\end{equation}
For example, if $\sigma$ is bounded, then we can take $H \equiv C^2_0$. Suppose that $X_t$ is a
positive recurrent process and $\pi$ is its unique stationary distribution, then
\begin{equation}\label{eq:BAR}
\int_{S} \mathcal{A}f(x) \pi(dx) = 0 \quad \forall \, f\in H
\end{equation}
Unfortunately, for the diffusion case the BAR condition (\ref{eq:BAR}) is not always both necessary
and sufficient to characterize the stationary distribution $\pi$, and further assumptions are
required. For the purpose of this paper, we will just assume that the BAR condition is both
necessary and sufficient, and we will elaborate from it.

Suppose that $\set{X_t\in S \subset \Real^k :t\geq 0}$ solves (\ref{eq:SDE}), that exists a
countable subset $\hat{S}$ dense in $S$, and that $H\subset C^2_0(S)$ contains a {\it countable}
dense subspace $ \hat{H} \equiv \set{h_1,h_2,\hdots}\subset H$. Using the
aggregation-relaxation-inner approximation procedure described by Hernandez and Lasserre
\cite{Hernandez-Lasserre} we can set the following LP
\begin{equation}\label{eq:LPnm}
\Prob_{nm} \equiv \left\{\min u |\sum_{i=1}^n \lambda_i \mathcal{A}h(x_i)\leq u \, \forall \,
h\,\in \hat{H}_m \,;\, -\sum_{i=1}^n \lambda_i \mathcal{A}h(x_i)\leq u \, \forall \, h\,\in
\hat{H}_m\,;\, \lambda\in \Lambda_n\right\}
\end{equation}
where $\hat{H}_m= \set{h_1,\hdots,h_m}$, and  $\Lambda_n \equiv \set{\lambda \in\Real^n|
\sum_{i=1}^n \lambda_i g(x_i)\leq M \,;\, e^t\lambda =1\,;\,\lambda\geq0}$ with
$g:S\rightarrow\Real$ nonnegative and continuous and $M>0$ such that $\exists\, K>0$, $r>0$ for
which $g(x)\geq K\,\forall\,|x|>r$, $x\in S$, and $E_{\pi}[g(x)]< M$. The idea here is to ensure
tightness for a sequence of solutions to $\Prob_{nm}$. Let $\Theta_{nm} $ denote the set of optimal
solutions to $\Prob_{nm}$.

\begin{theorem}
Let $\set{X_t\in S \subset \Real^k :t\geq 0}$ be a positive recurrent process that solves
(\ref{eq:SDE}) and let $\pi$ be its unique stationary distribution. Suppose that (\ref{eq:BAR}) is
both necessary and sufficient, and that both $\hat{S}$ and $\hat{H}$ exists. Then, there exist
sequences of integers $n(i)$ and $m(i)$ such that $\pi_{n(i)m(i)}\rightarrow \pi$ as
$i\uparrow\infty$, with $\pi_{mn}\in \Theta_{nm}$.
\end{theorem}

\noindent\proof We know there exists a sequence $\set{\mu_n}$ of distribution functions on $S$ such
that
\begin{itemize}
\item[(a)]  For every $n=1,\hdots$, $\mu_n$ has finite support $\set{x_1,\hdots,x_n}$,
that is, $\mu_n$ is of the form $\mu_n=\sum_{i=1}^n \beta^n_i\delta_{x_i}$, with $\beta^n_i\geq
0\,\forall \, i=1\hdots,n$, and $\sum_{i=1}^n \beta^n_i =1$.
\item[(b)] The sequence $\set{\mu_n}$ converges weakly to $\pi$.
\end{itemize}
From the definition of weak convergence (plus a truncation argument) we know that $E_{\mu_n}[g(x)]
\rightarrow E_{\pi}[g(x)]$, therefore there exists a $n^*$ for which $\beta_n\in \Lambda_n$
$\forall \,n\geq n^*$. Fix $m>0$. Consider $\pi_{mn}\equiv \set{\lambda_{mni} \,i=1\hdots n}$. For
$h\in\hat{H}_m$ we have that
\begin{equation}
\lim_{n\uparrow\infty}\abs{\sum_{i=1}^n \lambda_{mni}\mathcal{A}h(x_i)} \leq \lim_{n\uparrow\infty}
\abs{\sum_{i=1}^n \beta^n_i\mathcal{A}h(x_i)} = \lim_{n\uparrow\infty}\abs{\int_S
\mathcal{A}h(x)\mu_n(dx)}\rightarrow \abs{\int_S \mathcal{A}h(x)\pi(dx)}=0
\end{equation}
The first inequality comes from the optimality of $\pi_{mn}$ and the feasibility of $\mu_n$ for
$n\geq n^*$. The limit holds since $\mathcal{A}h(x)$ is bounded on the support of $h$ by its
continuity. Due to the tightness condition we have that there exists a subsequence $\set{n^m(i)}$
of integers such that $\pi_{nm}\Rightarrow\pi_m$, with $\pi_m$ a proper distribution function. From
the reasons above we can also conclude that
\begin{equation}
\lim_{n\uparrow\infty}\abs{\sum_{i=1}^n \lambda_{mni}\mathcal{A}h(x_i)} = \int_{S}
\mathcal{A}h(x)\pi_m(dx) =0
\end{equation}
Now consider the sequence $\set{\pi_m}$. Again, this sequence is tight, and therefore
$\pi_m\Rightarrow\hat{\pi}$ along a further subsequence. Finally we have that, for an arbitrary
function $h\in\hat{H}$ $\exists\, m^*$ such that $h \in \hat{H}_m$ $\forall \,m\geq m^*$ and we
have that $\abs{\int_S \mathcal{A}h(x)\pi_m(dx)} = 0$. Therefore
\begin{equation}
\abs{\int_S \mathcal{A}h(x)\hat{\pi(dx)}} = \lim_{m\uparrow \infty} \abs{\int_S
\mathcal{A}h(x)\pi_m(dx)} = 0
\end{equation}
Since $h$ was arbitrary, we conclude that $\hat{\pi}=\pi$ by the necessity and sufficiency of
(\ref{eq:BAR}), and the uniqueness of the stationary distribution. \qed
