Let $X$ be the constrained random walk on ${\mathbb Z}_+^d$ representing the queue lengths of a stable Jackson network and $x$ its initial position. Let $τ_n$ be the first time the sum of the components of $X$ equals $n$. $p_n \doteq P_x(τ_n < τ_0)$ is a key performance measure for the queueing system represented by $X$, stability implies $p_n\rightarrow 0$ exponentially. Currently the only analytic method available to approximate $p_n$ is large deviations analysis, which gives the exponential decay rate of $p_n$. Finer results are available via rare event simulation. The present article develops a new method to approximate $p_n$ and related expectations. The method has two steps: 1) with an affine transformation, move the origin onto the exit boundary of $τ_n$, take limits to remove some of the constraints on the dynamics, this yields a limit unstable constrained walk $Y$ 2) Construct a basis of harmonic functions of $Y$ and use them to apply the classical superposition principle of linear analysis. The basis functions are linear combinations of $\log$-linear functions and come from solutions of "harmonic systems," which are graphs whose vertices represent points on the "characteristic surface" of $Y$, the edges between the vertices represent conjugacy relations between the points, the loops represent membership in "the boundary characteristic surfaces." Using our method we derive explicit, simple and almost exact formulas for $P_x(τ_n < τ_0)$ for $d$-tandem queues, similar to the product form formulas for the stationary distribution of $X$. The same method allows us to approximate the Balayage operator mapping $f$ to $x \rightarrow {\mathbb E}_x \left[ f(X_{τ_n}) 1_{\{τ_n < τ_0\}} \right]$ for a range of stable constrained random walks in $2$ dimensions. We indicate how the ideas of the paper relate to more general processes and exit boundaries.