A survey on singular perturbations of Markov chains and decision processes

Konstantin E. Avrachenkov^{∗}, Jerzy Filar^{†}, and Moshe Haviv^{‡}

May 28, 2000

Abstract

In this survey we present a unified treatment of both singular and regular perturbations in finite Markov chains and decision processes. The treatment is based on the analysis of series expansions of various important entities such as the perturbed stationary distribution matrix, the deviation matrix, the mean-passage times matrix and others.

1 Background and Motivation

Finite state Markov Chains (MC’s) are among the most widely used probabilistic models of discrete event stochastic phenomena. Named after A.A. Markov, a famous Russian mathematician, they capture the essence of the existentialist “here and now” philosophy in the so-called “Markov property” which, roughly speaking, states that probability transitions to a subsequent state depend only on the current state and time. This property is less restrictive than might appear at first because there is a great deal of flexibility in the choice of what constitutes the “current state”. Because of their ubiquitous nature Markov Chains are, nowadays, taught in many undergraduate and graduate courses ranging from mathematics, through engineering to business administration and finance.

∗Department of Mathematics, The University of South Australia, The Levels, South Australia 5095, Australia.

†Department of Mathematics, The University of South Australia, The Levels, South Australia 5095, Australia. This research was supported in part by the ARC grant #A49906132.

‡Department of Statistics, The Hebrew University, 91905 Jerusalem, Israel and Department of Econometrics, The University of Sydney, Sydney, NSW 2006, Australia.

Whereas a MC often forms a good description of some discrete event stochastic process, it is not automatically equipped with a capability to model such a process in the situation where there may be a “controller” or a “decision-maker” who – by a judicious choice of actions – can influence the trajectory of the process. This innovation was not introduced until the seminal works of Howard ^{\underline{How}} [44] and Blackwell ^{\underline{Black}} [16] that are generally regarded as the starting point of the modern theory of Markov Decision Processes (or MDP’s for short). Since then, MDP’s have evolved rapidly to the point that there is a now a fairly complete existence theory, and a number of good algorithms for computing optimal policies with respect to criteria such as maximisation of limiting average expected reward, or the discounted expected reward.

The bulk of the, now vast, literature on both Markov Chains and Markov Decision Models deals with the “perfect information” situations where all the model parameters – in particular probability transitions – are assumed to be known precisely. However, in most applications this assumption will be violated. For instance, a typical parameter, ρ, would normally be replaced by an estimate

ρ̂ = ρ + ε(n)

where the error term, ε(n), comes from a statistical procedure used to estimate ρ and n is the number of observations used in that estimation. In most of the valid statistical procedures |ε(n)| ↓ 0 as n ↑ ∞, in an appropriate sense. Thus, from a perturbation analysis point of view, it is reasonable to suppress the argument n and simply concern ourselves with the effects of ε → 0.

Roughly speaking, the subject of perturbation analysis of MC’s and MDP’s divides naturally into the study of “regular” and “singular” perturbations. Intuitively, regular perturbations are “good” in the sense that the effect of the perturbation dissipates harmlessly as ε → 0, whereas singular pertubations are “bad” in the sense that small changes of ε (in a neighbourhood of 0) can induce “large” effects. Mathematically, it can be shown that singular perturbations are associated with a change of the rank of a suitably selected matrix. This can be easily seen from the now classical example due to Schweitzer ^{\underline{Sc3}} [76] where the perturbed probability transition matrix

P (ε) = ( \begin{matrix}1 − \frac{ε}{2} & \frac{ε}{2} \\ \frac{ε}{2} & 1 − \frac{ε}{2} \end{matrix} ) ⟶^{ε↓0} P (0)

but the stationary distribution matrix

P^{∗}(ε) = ( \begin{matrix}\frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & \frac{1}{2} \end{matrix} ) ↛ P^{∗}(0) = ( \begin{matrix}1 & 0 \\ 0 & 1 \end{matrix} ) .

Indeed, the rank of P^{∗}(ε) is 1 for all ε > 0 and near 0, but it increases to 2 at ε = 0; despite the fact that P (ε) → P (0). Thus we see that singular perturbations can occur in MC’s in a very natural and essential way. The latter point can be underscored by observing the behaviour of (I − λP)^{−1} as λ → 1, where P is a Markov (probability transition) matrix. It is wellknown (e.g. see ^{\underline{Black}} [16] or ^{\underline{MV}} [65]) that this inverse can be expanded as a Laurent series (with a pole of order 1) in the powers ε := 1 − λ. Indeed, much of the theory devoted to the connections between the discounted and limiting average MDP’s exploits the asymptotic properties of this expansion, as ε ↓ 0. Of course, the rank of (I − λP) changes when λ = 1(ε = 0).

It is not surprising, therefore, that the literature devoted to singularly perturbed MC’s and MDP’s has been growing steadily in recent years. In fact there have been quite a few developments since the 1995 survey by Abbad and Filar ^{\underline{AbFlBl}} [3].

The purpose of this survey paper is to present an up to date outline of a unified treatment of both singular and regular perturbations in MC’s and MDP’s that is based on series expansions of various important entities such as the perturbed stationary distribution matrix, the deviation matrix, the mean-passage times matrix and the resolvent-like matrix (I − λP)^{−1}. ¿From this series expansion perspective, the regular perturbations are simply the cases where Laurent series reduce to power series. Consequently, the capability to characterise and/or compute the coefficients of these expansions and the order of the pole (if any, at ε = 0) becomes of paramount importance.

This survey covers only the results on discrete time MC’s and MDP’s. For a parallel development of the continuous time models we refer an interested reader to the comprehensive book of Yin and Zhang ^{\underline{YZ}} [81].

The logical structure of the survey is as follows: in Section 2, perturbations of (uncontrolled) Markov chains are discussed from the series expansions perspective; in Section 3 the consequences of these results are discussed in the context of optimisation problems arising naturally in the (controlled) MDP case; and, finally, in Section 4 applications of perturbed MDP’s to the Hamiltonian Cycle Problem are outlined. This last section demonstrates that the theory of perturbed MDP’s has applications outside of its own domain.

2 Uncontrolled perturbed Markov chains

2.1 Introduction and preliminaries

Let P ∈ R^{n×n} be a transition stochastic matrix representing transition probabilities in a Markov chain. Suppose that the structure of the underlying Markov chain is aperiodic. Let P^{∗} = lim_{t→∞} P^{t} which is well-known to exist for aperiodic processes. In the case when the process is also ergodic, P^{∗} has identical rows, each of which is the stationary distribution of P, denoted by π. Let Y be the deviation matrix of P which is defined by Y = (I − P + P^{∗})^{−1} − P^{∗}. It is well known (e.g. see ^{\underline{KS}} [55]) that Y exists and it is the unique matrix satisfying Y (I − P) = I − P^{∗} = (I − P)Y and P^{∗}Y = 0 = Y 1 (where 1 here is a matrix full of 1’s) making it the group inverse of I − P. Finally, Y = lim_{T →∞} Σ_{t=0}^{T}(P^{t} − P^{∗}). Let M_{ij} be the mean passage time from state-i into state-j. It is also known that when the corresponding random variable is proper, then M_{ij} is finite. Of course, the matrix M is well-defined if and only if the Markov chain is ergodic. In this case, we have M_{ij} = (δ_{ij} + Y_{jj} − Y_{ij})/π_{j} and, in particular, M_{ii} = 1/π_{i}. The above mentioned results can be found in many sources, for instance, see Meyer ^{\underline{Me}} [64].

We consider (linear) perturbations of the matrix P and their impact on the structure of the process and on various essential matrices such as the stationary distribution matrix, the deviation matrix and the mean passage time matrix. Specifically, for a scalar ε, 0 < ε < ε_{max}, and for some zero rowsum matrix C, we look at the set of perturbed stochastic matrices P (ε) = P +εC which are assumed to be ergodic for any ε in the above mentioned region. Note that ergodicity is not assumed with regard to P = P (0). Actually, the case where P (0) contains some unrelated chains (with or without transient states) is our main focus. Our goal here is to survey the existing literature on series expansions for π(ε), P^{∗}(ε), Y (ε) and M(ε), which denote the stationary distribution, the limit matrix, the deviation matrix and the mean passage time matrix, repsectively, of P (ε), for 0 < ε < ε_{max} and consider their relationship to the corresponding entities in the unperturbed MC for P = P (0).

The rest of this section is organized as follows. In the next subsection we discuss the regular case, namely the case in which the unperturbed system is ergodic. Then, in Subsection 2.3, we look at the nearly completely decomposable (NCD) case, namely the case in which the state space under the unperturbed process is decomposed into a number of ergodic classes. This number is assumed here to be at least two and no transient states are allowed. This assumption is removed in Subsection 2.4, where we allow the unperturbed system to have two or more ergodic classes plus a number of transient states. It will be seen that the presence of transient states induces some interesting phenomena. In Subsections 2.2 through 2.4, we assume that the perturbed process in ergodic. Subsection 2.5, we comment on some issues involving the removal of this assumption.

2.2 The regular case

In this subsection we assume that the unperturbed Markov chain is ergodic. This leads to the case of regular perturbations. The following results appear in a seminal paper by Schweitzer ^{\underline{Sc1}} [74].

Theorem 1 Assume that the unperturbed Markov chain is ergodic. Then,

(i) The matrix functions P^{∗}(ε), Y (ε) and M(ε) are analytic in some (undeleted) neighborhood of zero. In particular, they all admit Maclaurin series expansions:

P^{∗}(ε) = ∑_{m=0}^{∞} ε^{m}P^{(∗m)} , Y (ε) = ∑_{m=0}^{∞} ε^{m}Y^{(m)} and M(ε) = ∑_{m=0}^{∞} ε^{m}M^{(m)}

with some coefficient sequences {P^{(∗m)}}_{m=0}^{∞}, {Y^{(m)}}_{m=0}^{∞} and {M^{(m)}}_{m=0}^{∞}.

(ii) The limit matrices P^{∗}(ε) and the deviation matrix of the perturbed Markov chain admit the following updating formulae

P^{∗}(ε) = P^{∗}(0)[I − εU]^{−1}

and

Y (ε) = [I − P^{∗}(ε)]Y (0)[I − εU]^{−1}
= Y (0)[I − εU]^{−1} − P^{∗}(0)[I − εU]^{−1}Y (0)[I − εU]^{−1},

where U := CY (0).

(iii) These updating formulae yield the following expressions for the power series coefficients.

P^{(∗0)} = P^{∗}(0), Y^{(0)} = Y (0), M^{(0)} = M(0)
P^{(∗m)} = P^{(∗0)}U^{m}, m ≥ 0

Y^{(m)} = Y (0)U^{m} − P^{∗}(0) ∑_{j=1}^{m} U^{j}Y (0)U^{m−j}, m ≥ 0

M_{ij}^{(m)} = \frac{1}{π_{j}^{(0)}} (Y_{jj}^{(m)} − Y_{ij}^{(m)}) − \frac{1}{π_{j}^{(0)}} ∑_{l=1}^{m} π_{j}^{(l)}M_{ij}^{(m−l)}, m ≥ 0

(iv) The validity of any of the above series expansion holds for any ε, 0 ≤ ε < min{ε_{max}, ρ^{−1}(U)} where ρ(U) is the spectral radius of U.

Example 1. For 0 ≤ ε < .25 let

P (ε) = P (0) + εC = ( \begin{matrix}.5 & .5 \\ .5 & .5 \end{matrix} ) + ε ( \begin{matrix}2 & −2 \\ −1 & 1 \end{matrix} ) .

Clearly,

P^{∗}(0) = P^{(∗0)} = ( \begin{matrix}.5 & .5 \\ .5 & .5 \end{matrix} )

Also,

Y (0) = Y^{(0)} = ( \begin{matrix}.5 & −.5 \\ −.5 & .5 \end{matrix} ) and M(0) = M^{(0)} = ( \begin{matrix}2 & 2 \\ 2 & 2 \end{matrix} ) .

Hence,

U = CY (0) = ( \begin{matrix}2 & −2 \\ −1 & 1 \end{matrix} ) .

It is easy to see that for m ≥ 1, U^{m} = 3^{m−1}U and hence for m ≥ 1,

P^{(∗m)} = P^{∗}(0)U^{m} = 3^{m−1} ( \begin{matrix}.5 & −.5 \\ .5 & −.5 \end{matrix} ) .

Also, for m ≥ 1,

Y^{(m)} = 3^{m−1}Y (0)U − 3^{m−2}(m − 1)P^{∗}(0)UY (0)U − 3^{m−1}P^{∗}(0)UY (0)

= 3^{m−1} ( \begin{matrix}1.5 & −1.5 \\ −1.5 & 1.5 \end{matrix} )−3^{m−2}(m−1) ( \begin{matrix}1.5 & −1.5 \\ 1.5 & −1.5 \end{matrix} )−3^{m−1} ( \begin{matrix}.5 & −.5 \\ .5 & −.5 \end{matrix} ) .

Finally,

M_{12}(ε) = 2 + 8ε + ... and M_{21}(ε) = 2 + 4ε + ...

2.3 The nearly completely decomposable case

Let P (0) ∈ R^{n×n} be a stochastic matrix representing transition probabilities in a completely decomposable Markov chain. By the latter we mean that there exists a partition Ω of the state space into p, p ≥ 2, subsets Ω = {I_{1}, …, I_{p}} each of which being an ergodic class. We assume that the order of the rows and of the columns of P is compatible with Ω, i.e., for p stochastic matrices, P_{I_{1}}, …, P_{I_{p}},

P (0) = ( \begin{matrix}P_{I_{1}} & 0 & ⋯ & 0 \\ 0 & P_{I_{2}} & ⋯ & 0 \\ ⋮ & ⋮ & ... & 0 \\ 0 & 0 & ⋯ & P_{I_{p}} \end{matrix} )

Note that we assume above that none of the states is transient. Let C ∈ R^{n×n} be a zero rowsum matrix such that for some ε_{max} > 0, the matrix P + εC is stochastic for ε ∈ (0, ε_{max}) representing transition probabilities in an ergodic Markov chain. For small values of ε, P (ε) is called nearly completely decomposable (NCD) or sometimes nearly uncoupled. Clearly, C_{ij} ≥ 0 for any pair of states i and j belonging to different subsets.

Probably, the first motivation to study the singular perturbed Markov chains was given by Simon and Ando ^{\underline{SA}} [77]. They demonstrated that several problems in econometrics lead to the mathematical model based on singularly perturbed Markov chains. The first rigorous theoretical developments of the singularly perturbed Markov chains have been carried out by Pervozvanski and Smirnov ^{\underline{PS}} [68] and Gaitsgori and Pervozvanskii ^{\underline{GP}} [32]. In particular, they have shown that the limiting probability distribution π^{(0)} = lim_{ε→0} π(ε) can be expressed in terms of the invariant probability distributions of the ergodic classes I_{k}, k = 1, ..., p and of the stationary distribution of an aggregated chain. Similar ideas were also developed in works of Courtois and his co-authors ^{\underline{Co,CL,CS}} [19, 20, 21] and in the work of Haviv and his co-authors ^{\underline{Ha2,Ha3,Ha1, HR1,HR2,HV}} [35, 36, 37, 39, 40, 45]. Schweitzer ^{\underline{Sc2}} [75] showed that π(ε) is analytic in some deleted neighbourhood of zero. That is, π(ε) = ∑_{m=0}^{∞} π^{(m)}ε^{m}, where π^{(0)} = lim_{ε→0} π(ε) and where π^{(m)}, m ≥ 1, are zerosum vectors. Note that π_{i}^{(0)} > 0 for all i. Moreover, {π^{(m)}}_{m=0}^{∞} is a geometric series, that is, for some matrix U ∈ R^{n×n}, π^{(m)} = π^{(0)}U^{m}, 0 ≤ m < ∞. See (^{\underline{UNCD}} 4) below for an explicit expression for U. Finally, the series expansion holds for 0 < ε < max{ε_{max}, ρ^{−1}(U)}.

For any subset I ∈ Ω, let

k_{I} = Σ_{i∈I}π_{i}^{(0)}  (1) k

and let γ_{I} be the subvector of π^{(0)} corresponding to subset I rescaled so as its entry-sum is now one. Then, γ_{I} is the unique stationary distribution of P_{I}. Note that computing γ_{I} is relatively easy as only the knowledge of P_{I} is needed.

Next define the matrix Q ∈ R^{p×p} which is

usually referred to as the aggregate transition matrix. Each row, and likewise each column in Q corresponds to a subset in Ω. Then, for subsets I and J, I ≠ J, let

Q_{IJ} = ∑_{i∈I}(γ_{I})_{i} ∑_{j∈J} C_{ij}  (2) QIJ

and let

Q_{II} = 1 + ∑_{i∈I}(γ_{I})_{i} ∑_{j∈I} C_{ij} = 1 − ∑_{J≠I} Q_{IJ} .  (3) QII

Without loss of generality, assume that Q_{II} is non-negative for all subsets I and hence Q is easily seen to be a stochastic matrix.^{1} Moreover, Q is irreducible and the vector k ∈ R^{p} (see (^{\underline{k}}1)) is easily checked to be its unique stationary distribution. Often it is convenient to express the aggregated transition matrix Q in matrix terms. Specificaly, let V ∈ R^{p×n} be such that its i-th row is full of zeros except for γ_{I_{i}} at the entries corresponding to subset I_{i}, and where W ∈ R^{n×p} is such that its j-th column is full of zeros except for 1’s in the entries corresponding to subset I_{j}.^{2} Then, we can write

Q = I + V CW.

The aggregate stochastic matrix Q represents transition probabilities between subsets which in this context are sometimes referred to as macrostates. However, although the original process among states is Markovian, this is not necessarily the case with the process among macro-states (and indeed typically it is not).^{3} Yet, as the following indicates, much can be learned on the original process from the analysis of the aggregate matrix.

1Note that the matrix C can be divided by any constant and ε can be multiplied by this constant leading to the same n × n transtion matrices. Taking this constant small enough guarantees the stochasticity of Q and hence this is assumed without loss of generality. In particular, the stationary distribution of Q is invariant with respect to the choice of this constant. Alternatively, one can define Q_{II} by −Σ_{J≠I}Q_{IJ} and consider Q as the generator of the aggregated process, i.e., the process among subsets (and hence no need to assume anything further with regard to the size of the entries of the matrix C.)

2Note that V W ∈ R^{p×p} is the identity matrix. Moreover, V and W correspond to orthonormal sets of eigenvectors of P (0) belonging to the eigenvalue 1, V as left eigenvectors and W as right eigenvectors.

3The process among macro-states is an example of a partially observable Markov process.

Theorem 2 Let the perturbed Markov chain be nearly completely decomposable. The stationary distribution π(ε) admits a Maclaurin series expansion in a deleted neighborhood of zero. Specifically, for some vectors {π^{(m)}}_{m=0}^{∞} with π^{(0)} being a probability vector positive in all its entries and satisfying π^{(0)} = π^{(0)}P (0), and for some zerosum vectors π^{(m)}, m ≥ 1, π(ε) = Σ_{m=0}^{∞}π^{(m)}ε^{m}. Moreover, for I ∈ Ω, π_{I}^{(0)} = k_{I}γ_{I}.^{4} Also, the series {π^{(m)}}_{m=0}^{∞} is geometric, i.e., for some square matrix U, π^{(m)} = π^{(0)}U^{m} for any m ≥ 0. Actually,

U = CY (0)(I + CW DV),  (4) UNCD

where D is the deviation matrix of the aggregated transition matrix Q. Alternatively,

U = CY^{(0)}  (5) Alt

where Y^{(0)} is the first regular term of the Laurent series expansion for Y (ε) (see also the next theorem). Finally, the validity of the series expansion holds for any ε, 0 ≤ ε < min{ε_{max}, ρ^{−1}(U)} where ρ(U) is the spactral radius of U.

We note that a series expression for π(ε) appeared orginally in ^{\underline{Sc2}} [75]. The expression for U in (^{\underline{UNCD}} 4) taken from ^{\underline{HR1}} [39] sheds more light on the role of the aggregate process played in the original process than its original expression given in ^{\underline{Sc2}} [75] (see Equation 3-1 there). Also we note that U = CY^{(0)} is an elegant generalization ^{\underline{Sc3}} [76] of the regular case where U = CY (0).

For ε, 0 ≤ ε < ε_{max}, let Y (ε) be the deviation matrix of P (ε). This matrix is uniquely defined and the case ε = 0 is no exception. Yet, as we see shortly, there is no continuity of Y (ε) at ε = 0. In particular, Y (0) has the same shape as P has, namely

Y (0) = ( \begin{matrix}Y_{I_{1}} & 0 & ⋯ & 0 \\ 0 & Y_{I_{2}} & ⋯ & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & ⋯ & Y_{I_{p}} \end{matrix} )  (6) Y(0)

where Y_{I_{i}} is the deviation matrix of P_{I_{i}}, 1 ≤ i ≤ p.

Theorem 3 In the case of NCD Markov chains, the matrix Y (ε) admits a Laurent series expansion in a deleted neighborhood of zero with the order of the pole being exactly one. Specifically, for some matrices {Y^{(m)}}_{m=−1}^{∞} with Y^{(−1)} ≠ 0,

4Recall that γ_{I} is the stationary distribution of P_{I} and that k in the stationary distribution of the aggreagted matrix Q as defined in (\underline{QIJ} 2) and (\underline{QII} 3).

Y (ε) = \frac{1}{ε}Y^{(−1)} + Y^{(0)} + εY^{(1)} + ε^{2}Y^{(2)} + ⋯  (7) Yexp

for 0 < ε < ε_{max}. Moreover,

Y^{(−1)} = W DV,

or in a component form,

Y_{ij}^{(−1)} = D_{IJ}(γ_{J})_{j} , i ∈ I , j ∈ J.  (8) Y-1

A series expression for Y (ε) was originally developed in ^{\underline{Sc2}} [75]. Yet, the expression for Y^{(−1)} given above in (^{\underline{Y-1}} 8) is taken from ^{\underline{HR1}} [39]. We find (^{\underline{Y-1}} 8) appealing as it explicitly shows the role played by the aggregate matrix (and hence by the process among macro-states) in the original process. This point was already mentioned in ^{\underline{Ha1}} [37] (see Equation 3 there).

We now focus our attention on M(ε). Note that as opposed to Y (0), M(0) is not well-defined as the corresponding mean value (when ε = 0 and states i and j belong to two different subsets) does not exist. Let E ∈ R^{p×p} be the mean passage time matrix associated with the aggregated process. I.e., for any pair of subsets I and J (I = J included), E_{IJ} is the mean passage time from the macro-state I into the macro- state J when transition probabilities are governed by the stochastic matrix Q.

The following theorem appears in ^{\underline{AH}} [10].

Theorem 4 The matrix M(ε) admits a Laurent series expansion in a deleted neighborhood of zero with the order of the pole being exactly one. Specifically, for some matrices {M^{(m)}}_{m=−1}^{∞} with M^{(−1)} ≠ 0,

M(ε) = \frac{1}{ε}M^{(−1)} + M^{(0)} + εM^{(1)} + ε^{2}M^{(2)} + ⋯  (9) Mseries

for 0 < ε < ε_{max}. Moreover, for i ∈ I and j ∈ J,

M_{ij}^{(−1)} = \{ \begin{matrix}0 & if J = I \\ E_{IJ} & ifJ ≠ I \end{matrix}  (10) M-1

M_{ij}^{(m)} = \frac{1}{π_{j}^{(0)}} (Y_{jj}^{(m)} − Y_{ij}^{(m)}) − \frac{1}{π_{j}^{(0)}} ∑_{l=1}^{m+1} π_{j}^{(l)}M_{ij}^{(m−l)} , m ≥ −1

The following example is taken from ^{\underline{Sc3}} [76].

Example 2 Let

P (0) = ( \begin{matrix}1 & 0 & 0 \\ 0 & .5 & .5 \\ 0 & .5 & .5 \end{matrix} ) and C = \frac{2}{7} ( \begin{matrix}−2 & 1 & 1 \\ 3 & −1 & −2 \\ 4 & −3 & −1 \end{matrix} ) .

The number of subset equals 2 with γ_{I_{1}} = 1 and γ_{I_{2}} = (0.5, 0.5). First, we construct the following matrices.

V = ( \begin{matrix}1 & 0 & 0 \\ 0 & 0.5 & 0.5 \end{matrix} ) W = ( \begin{matrix}1 & 0 \\ 0 & 1 \\ 0 & 1 \end{matrix} ) Y (0) = ( \begin{matrix}0 & 0 & 0 \\ 0 & 0.5 & −0.5 \\ 0 & −0.5 & 0.5 \end{matrix} )

The aggregated transition matrix is given by

Q = I + V CW = ( \begin{matrix}3/7 & 4/7 \\ 1 & 0 \end{matrix} ) .

Hence, k = (7/11, 4/11) and π^{(0)} = (7/11, 2/11, 2/11). Next, we calculate D, the deviation matrix of Q,

D = (I − Q + Q^{∗})^{−1} − Q^{∗} = \frac{1}{121} ( \begin{matrix}28 & −28 \\ −49 & 49 \end{matrix} )

and hence, using (^{\underline{Y-1}} 8),

Y^{(−1)} = W DV = \frac{1}{242} ( \begin{matrix}56 & −28 & −28 \\ −98 & 49 & 49 \\ −98 & 49 & 49 \end{matrix} )

The matrix E, which is the mean passage time matrix for the aggregated process, equals

E = ( \begin{matrix}11/7 & 7/4 \\ 1 & 11/4 \end{matrix} )

and hence, using (^{\underline{M-1}} 10),

M^{(−1)} = ( \begin{matrix}0 & 7/4 & 7/4 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \end{matrix} )

Finally, from (^{\underline{UNCD}} 4) we get that

U = CY (0)(I + CW DV) = \frac{1}{77} ( \begin{matrix}0 & 0 & 0 \\ −2 & 12 & −10 \\ 4 & −24 & 20 \end{matrix} )

2.4 The general case

In this section we impose no assumptions on the ergodic structure of the unperturbed chain. That is, the unperturbed chain may now consist of several ergodic classes plus a set of transient states. Probably, Delebecque and Quadrat ^{\underline{DqQt}} [24] were the first who studied the model in this general setting. However, they have analysed only the case of the first order singularity, namely the case where the Laurent series expansion for the perturbed deviation matrix has a simple pole. This analysis is not much different from the analysis of the NCD case. Later Delebecque ^{\underline{D}}[25], using the reduction technique of Kato ^{\underline{Ka}} [52], removed the first order singularity assumption and showed how to compute the asymptotic expansion for π(ε) in the general case. However, the reduction technique of ^{\underline{Ka}} [52] does not provide an efficient computational scheme. Other related papers devoted to the case where the unperturbed chain has a transient set are ^{\underline{CW,CdWl,BlSt,KT,RW,Vat}} [17, 18, 15, 56, 72, 78].

The following typical example of the general situation was considered extensively in the past (see ^{\underline{HR1, HH, HR2, AL}} [39, 34, 40, 12]):

Example 3.

P (ε) = P (0) + εC = ( \begin{matrix}0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \end{matrix} ) + ε ( \begin{matrix}0 & −1 & 0 & 1 \\ 1 & −1 & 0 & 0 \\ 0 & 1 & 0 & −1 \\ 0 & 0 & 1 & −1 \end{matrix} ) .

In this example the unperturbed chain contains two ergodic classes (states 2 and 4) and two transient states (states 1 and 3). They all are coupled to a single chain when ε > 0. Moreover, states 2 and 4 (i.e., the ergodic chains in the unperturbed process) communicate under the perturbed case only via states 1 and 3 (i.e., transient states in the unperturbed case). This phonomenon makes this case more involved than the NCD case. For example, the reader is invited to check that the order of magnitude of M_{24}(ε) is ε^{−2}.

Theorem 5

(i) The stationary distribution π(ε) admits a Maclaurin series expansion in a deleted neighbourhood of zero, that is, for some sequence {π^{(m)}}_{m=0}^{∞} with π^{(0)} being a probability vector and π^{(m)} for m ≥ 1, being zero sum vectors,

π(ε) = ∑_{m=0}^{∞} ε^{m}π^{(m)} .

Moreover, π^{(0)}P (0) = π^{(0)} and for some matrix U,

π^{(m)} = π^{(0)}U^{m} .  (11) pi_geom

(ii) The deviation matrix Y (ε) admits a Laurent series expansion in a deleted neighborhood of zero with a non-essential pole, that is, for some nonnegative integer s and for some sequence of matrices {Y^{(m)}}_{m=−s}^{∞} with Y^{(−s)} ≠ 0,

Y (ε) = \frac{Y^{(−s)}}{ε^{s}} + \frac{Y^{(−s+1)}}{ε^{s−1}} + ⋯ + Y^{(0)} + Y^{(1)}ε + · · · .

Moreover, the regular part Y^{R}(ε) = Y^{(0)} + εY^{(1)} + ⋯ can be expressed by the updating formula

Y^{R}(ε) = [I − P^{∗}(ε)]Y^{(0)}[I − εU]^{−1} − P^{∗}(ε) ∑_{i=1}^{s} U^{i}Y^{(−i)} . (12) up_Y

The latter provides a computationally efficient recursive formula for the coefficients of the regular part:

Y^{(m)} = Y^{(0)}U^{m} − P_{0}^{∗} ∑_{j=0}^{m} U^{j}Y^{(0)}U^{m−j}

+P^{∗(0)}U^{m}[I − ∑_{i=1}^{s} U^{i}Y^{(−i)}], m ≥ 1.  (13) recur_Y

Furthermore, the quotient matrix U can be expressed as

U = CY^{(0)}  (14) Y0G

(iii) The mean passage time matrix M(ε) admits a Laurent series expansion in a deleted neighborhood of zero with a non-essential pole, that is, for some nonnegative integer t and for some sequence of matrices {M^{(m)}}_{m=−t}^{∞} with M^{(−t)} ≠ 0,

M(ε) = \frac{M^{(−t)}}{ε^{t}} + \frac{M^{(−t+1)}}{ε^{t−1}} + ⋯ + M^{(0)} + M^{(1)}ε + · · · .

Finally, t ≥ s.

We would like to emphasize that once one has found the value of s and computed the terms Y^{(−s)}, ..., Y^{(0)} and π^{(0)}, all the other terms in the series expansion for Y (ε) and π(ε) can be computed by (^{\underline{recur_Y}} 13). We refer below to a procedure for computing s. The term π^{(0)} (which we believe to be the most important coefficient of the expansion) can be computed either by a reduction process ^{\underline{D}}[25] or via the determination of the null space of an auxiliary augmented matrix. The latter method will be described in some detail below. Finally, one can find methods for computing Y^{(−s)}, ..., Y^{(0)} (and hence U) (once s and π^{(0)} are in hand) by methods outlined in ^{\underline{AL}} [12] and ^{\underline{HR1}} [39].

Remark 1 Note that the formulae (^{\underline{pi_geom}} 11) and (^{\underline{up_Y}} 12) are the natural generalization of the regular case. Formula (^{\underline{pi_geom}} 11) and (^{\underline{Y0G}} 14) appeared originally in \underline{Sc3} [76]. Finally, (^{\underline{up_Y}} 12) was derived in ^{\underline{AL}} [12].

The next theorem shows that the order of poles of Y (ε) and M(ε) (denoted above by s and t, respectively), can be computed by a combinatorial algorithm.

Theorem 6 For 1 ≤ i, j ≤ n, let u_{ij} be the order of the pole of M_{ij}(ε) at zero. Then, there exists an O(n^{4}) algorithm for computing these orders whose input is the addresses of the nonzero entries in P (0) and in C. Likewise, for 1 ≤ i, j ≤ n, let v_{ij} be the order of the pole of Y_{ij}(ε) at zero. Then,

t ≥ max_{ij}(u_{ij} − u_{jj}) = max_{ij} v_{ij} = s .

The above theorem and a detailed algorithm are given in ^{\underline{HH}} [34].

Example 3 (cont.). Running the algorithm developed in ^{\underline{HH}} [34] for computing u_{ij}, 1 ≤ i, j ≤ 4, results in

(u_{ij}) = ( \begin{matrix}1 & 1 & 2 & 2 \\ 1 & 0 & 2 & 2 \\ 2 & 2 & 1 & 1 \\ 2 & 2 & 1 & 0 \end{matrix} )

In particular, s = t = 2.

For m ≥ 0, let B_{m} ∈ R^{n(m+1)×n(m+1)} be the matrix such that each of its n × n block diagonal matrices equals I − P (0) and each of its n × n matrices above the diagonal equals −C. Namely,

B_{m} = ( \begin{matrix}I − P (0) & −C & 0 & ⋯ & 0 & 0 \\ 0 & I − P (0) & −C & ⋯ & 0 & 0 \\ 0 & 0 & I − P (0) & ⋯ & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ... & ⋮ & ⋮ \\ 0 & 0 & 0 & ⋯ & I − P (0) & −C \\ 0 & 0 & 0 & ⋯ & 0 & I − P (0) \end{matrix} ) .

where I −P (0) appears m+1 times while −C appears m times. By equating coefficients of the same order in the identity π(ε) = π(ε)(I − P (0) − εC), it is easy to conclude that for any m ≥ 0, (π^{(0)}, π^{(1)}, ⋯, π^{(m)}) is a left eigenvector of B_{m} belonging to the eigenvalue zero. Of course, it is not an unique eigenvector. Also, the dimension of this eigenspace of B_{m} is nonincreasing in m. Let V_{m} be the subspace of vectors of length n such that the first n entires of vectors in the above mentioned eigenspace of B_{m} induce. Again, the dimension of V_{m} is non-increasing in m. Finally, these dimensions are always greater than or equal to 1. The following theorem is taken from \underline{HR2} [40].

Theorem 7 (The general case: computing π^{(0)}).

t + 1 = min{m; dimV_{m} = 1} .

The above two theorems suuggest a way to compute π^{(0)}. Specifically, one can first find the value of t by the combinatorial algorithm suggested in \underline{HH} [34] and then the left eigenspace of B_{t} belonging to the eigenvalue zero can be explored. For the latter task some procedures exist. A few of them are based on reduction steps. In ^{\underline{HR1}} [39] one reduction step is carried out while ^{\underline{D}}[25] and ^{\underline{AL}} [12] suggest performing a number of reduction steps. The latter number is bounded by t as at the t-th step one gets a non-singular system. Actually, the procedure was defined in the previous subsection and takes care of NCD matrices with one (and terminal) reduction step.

We again note that once π^{(0)} and Y^{(0)} are calculated, the other coefficients of the series expansion for π(ε) can be computed by the recurssive formula π^{(m)} = π^{(m−1)}U, m ≥ 1, where U = CY^{(0)}.

Example 3 (cont.). We have pointed out above that t = 2. Hence, solving for the left null space of

B_{2} = ( \begin{matrix}I − P (0) & −C & 0 \\ 0 & I − P (0) & −C \\ 0 & 0 & I − P (0) \end{matrix} )

= ( \begin{matrix}1 & −1 & 0 & 0 & 0 & 1 & 0 & −1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & −1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & −1 & 0 & −1 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & −1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & −1 & 0 & 0 & 0 & 1 & 0 & −1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & −1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & −1 & 0 & −1 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & −1 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & −1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & −1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{matrix} )

results in (0, α, 0, α, ∗, ∗, ∗, ∗, ∗, ∗, ∗, ∗) for any α.^{5} Hence, π^{(0)} = (0, .5, 0, .5). In ^{\underline{HR1}} [39] it is shown that the same set of equations which is needed in order to solve for π^{(0)} is also all required in order to compute U and Y^{(−s)}. For this numerical example, using the methods of either ^{\underline{HR1}} [39] or ^{\underline{AL}} [12], one obtains

Y^{(0)} = ( \begin{matrix}.5 & −.5 & .5 & −.5 \\ 0 & 0 & 0 & 0 \\ .5 & −.5 & .5 & −.5 \\ 0 & 0 & 0 & 0 \end{matrix} )

and hence

U = CY^{(0)} = ( \begin{matrix}0 & 0 & 0 & 0 \\ 0.5 & −0.5 & 0.5 & −0.5 \\ 0 & 0 & 0 & 0 \\ 0.5 & −0.5 & 0.5 & −0.5 \end{matrix} )

The terms of the singular part of the Laurent expansion for the deviation matrix are:

Y^{(−2)} = ( \begin{matrix}0 & .25 & 0 & −.25 \\ 0 & .25 & 0 & −.25 \\ 0 & −.25 & 0 & .25 \\ 0 & −.25 & 0 & .25 \end{matrix} ) and Y^{(−1)} = ( \begin{matrix}.25 & −.5 & −.25 & .5 \\ .25 & 0 & −.25 & 0 \\ −.25 & .5 & .25 & −.5 \\ −.25 & 0 & .25 & 0 \end{matrix} )

Finally,

M_{21}(ε) = \frac{Y_{11}(ε) − Y_{21}(ε)}{π_{1}(ε)} = [(\frac{0.25}{ε}+0.5+...)−(\frac{0.25}{ε}+0+...)]/[0.5ε−0.5ε^{2}+...] = \frac{1}{ε}+...,

5The sign ‘*’ corresponds to values not computed (and not relevant for our current purposes.)

M_{31}(ε) = \frac{1}{ε^{2}} + \frac{1}{ε} + ..., M_{41}(ε) = \frac{1}{ε^{2}} + \frac{2}{ε} + ...,

M_{12}(ε) = \frac{1}{ε} + ..., M_{32}(ε) = \frac{1}{ε^{2}} + \frac{0}{ε} + ..., M_{42}(ε) = \frac{1}{ε^{2}} + \frac{1}{ε} + ... .

We like to point out that the expression for Y^{(−1)} for this example given in \underline{HR1} [39] had a numerical error that was corrected in ^{\underline{AL}} [12].

2.5 The case of non ergodic perturbed Markov chains

In this section we briefly discuss the singularity phenomena that can appear if the perturbed Markov chain itself is not ergodic. The reader can find a more detailed exposition in ^{\underline{A}}[9].

In case where the perturbed Markov chain has several ergodic classes plus a non- empty set of transient states, then the perturbed transition matrix can be written in the following canonical form

P (ε) = [ \begin{matrix}P_{1}(ε) & ⋯ & 0 & 0 \\ ⋮ & ... & ⋮ & ⋮ \\ 0 & ⋯ & P_{ℓ}(ε) & 0 \\ R_{1}(ε) & ⋯ & R_{ℓ}(ε) & S(ε) \end{matrix} ] \} Ω_{1} ; ⋮ ; \} Ω_{ℓ} ; \} Ω_{T}  (15) can_pert

where P_{i}(ε) = P_{i}(0) + εC_{i}, R_{i}(ε) = R_{i}(0) + εC_{Ri}, S(ε) = S(0) + εC_{S}. Now note that all invariant measures m_{i}(ε) of the perturbed Markov chain can be immediately constructed from the invariant measures of the ergodic classes associated with stochastic matrices P_{i}(ε), i = 1, ..., ℓ. Namely, m_{i}(ε) = [0 ⋯ 0 π_{i}(ε) 0 ⋯ 0], where π_{i}(ε) is uniquely determined by the system

\{ \begin{matrix}π_{i}(ε)P_{i}(ε) = π_{i}(ε), \\ π_{i}(ε)\underline{1} = 1. \end{matrix}

This is exactly the perturbation problem under the assumption that there are several ergodic classes plus possibly a set of transient states coupled by the perturbation into a single ergodic class. This issue was covered in the previous subsections.

The challenging task that we face in this case of multi-chain perturbed process is the computation of power series for the right eigenvectors of the perturbed Markov chains. The pertubed right eigenvectors contain the probabilities of being absorbed in each of the ergodic classes under the perturbed processes. These quantities exhibit the multi-chain ergodic properties of the perturbed process. For instance, nothing conceptually new would emerge had the set of transient states in the perturbed process been empty.

The right 0-eigenvectors of the perturbed chain can be written in the following form ^{\underline{KS}} [55]

q_{i}(ε) = [ 0 ; \underline{1} ; 0 ; φ_{i}(ε) ] \begin{matrix}\} Ω_{i} \\ \} Ω_{T} \end{matrix}  (16) reig_pert

where subvector φ_{i}(ε) is given by

φ_{i}(ε) = (I − S(ε))^{−1}(ε)R_{i}(ε)\underline{1},  (17) phi_pert

Note that if some ergodic classes become transient after the perturbation, then the matrix valued function (I − S(ε))^{−1} has a singularity at ε = 0. To further elaborate on this phenomenon, let us consider the structure of the substochastic matrix S(0).

S(0) = [ \begin{matrix}P̃_{1} & ⋯ & 0 & 0 \\ ⋮ & ... & ⋮ & ⋮ \\ 0 & ⋯ & P̃_{m} & 0 \\ R̃_{1} & ⋯ & R̃_{m} & S̃ \end{matrix} ] .

Blocks P̃_{1}, ..., P̃_{m} represent the ergodic classes of the original Markov chain which merge with the transient set after the perturbation. Of course, m = 0 is possible. Since each of P̃_{1}, ..., P̃_{m} is a stochastic matrix, we conclude that matrix I − S(0) has zero as an eigenvalue with multiplicity of at least m. Of course, I − S(0) is not invertible. However, the matrix (I − S(ε))^{−1} exists for small positive (but not zero) values of ε. ¿From the results of ^{\underline{AHH}} [11], it follows that one can expand (I − S(ε))^{−1} as a Laurent series at ε = 0

(I − S(ε))^{−1} = \frac{1}{ε^{p}} U^{(−p)} + ... + \frac{1}{ε}U^{(−1)} + U^{(0)} + εU^{(1)} + ... (18) laur_T

One can use the methods of ^{\underline{AHH}} [11] to calculate the coefficients of the above series. Substituting the expression R_{i}(ε) = R_{i}(0) + εC_{Ri} and the Laurent series (^{\underline{laur_T}} 18) into formula (^{\underline{phi_pert}} 17), we obtain the asymptotic expansion

φ_{i}(ε) = φ_{i}^{(0)} + εφ_{i}^{(1)} + ε^{2}φ_{i}^{(2)} + ... ,  (19) pow_phi

where

φ_{i}^{(m)} = U^{(m)}R_{i}^{(0)} + U^{(m−1)}C_{Ri}, m ≥ 0.  (20) phi_coef

The above formulae are valid in the most general setting. Some interesting particular cases are discussed in ^{\underline{A}}[9].

3 Singularly Perturbed Markov Decision Processes

In this section we review some results on singular perturbations of Markov decision processes (MDPs) and introduce a unified approach to singular perturbation, Blackwell optimality and branching Markov decision processes. First we briefly intoduce notation and define various optimality criteria. The reader can find a more detailed study on MDPs in the books and surveys \underline{Der,HL,How,Kal,Kal2,Put,Vt2} [26, 42, 44, 53, 54, 70, 80] and the references therein.

We consider a discrete-time MDP with a finite state space 𝕏 = {1, ..., N} and a finite action space 𝔸(i) = {1, ..., m_{i}} for each state i ∈ 𝕏. At any time point t the system is in one of the states i ∈ 𝕏 and the controller or “decision-maker” chooses an action a ∈ 𝔸(i); as a result the following occur:

(a) the controller gains an immediate reward r(i, a), and (b) the process moves to a state j ∈ 𝕏 with probability p(j|i, a), where p(j|i, a) ≥ 0 and ∑_{j∈𝕏} p(j|i, a) = 1.

A decision rule π_{t} at time t is a function which assigns a probability to the event that any particular action a is taken at time t. In general, π_{t} may depend on history h_{t} = (i_{0}, a_{0}, i_{1}, a_{1}, ..., a_{t−1}, i_{t}) up to time t. The distribution π_{t}(·|h_{t}) defines the probability of selecting any action a at time t given the history h_{t}.

A policy (or strategy) is a sequence of decision rules π = (π_{0}, π_{1}, ..., π_{t}, ...). A policy π is called Markov if π_{t}(·|h_{t}) = π_{t}(·|i_{t}). If π_{t}(·|i) = π_{t′}(·|i) for all t, t′ ∈ ℕ then the Markov policy π is called stationary. Furthermore, a deterministic policy π is a stationary policy whose single decision rule is nonrandomized. It can be defined by the function f(i) = a, a ∈ 𝔸(i).

Let Π, Π^{S} and Π^{D} denote the sets of all policies, of all stationary policies and of all deterministic policies, respectively. For the stationary policy π ∈ Π^{S} define the corresponding transition matrix P (π) = {p_{ij}(π)}_{i,j=1}^{N} and the reward vector r(π) = {r_{i}(π)}_{i=1}^{N}

p_{ij}(π) := ∑_{a∈𝔸(i)}p(j|i, a)π(a|i), r_{i}(π) := ∑_{a∈𝔸(i)}r(i, a)π(a|i).

Let the limit matrix P^{∗}(π) and the deviation matrix (or reduced resolvent as it is sometimes refer to in this setting) Y (π) associated with policy π be defined by

P^{∗}(π) := lim_{T →∞} \frac{1}{T} ∑_{t=1}^{T} P^{t−1}(π)

and

Y (π) := (I − P (π) + P^{∗}(π))^{−1} − P^{∗}(π).

The expected average reward g_{i}(π) and the expected discounted reward v_{i}^{λ}(π) are defined as follows:

g_{i}(π) := lim_{T →∞} \frac{1}{T} ∑_{t=1}^{T} [P^{t−1}(π)r(π)]_{i}

and

v_{i}^{λ}(π) := ∑_{t=1}^{∞} λ^{t−1} [P^{t−1}(π)r(π)]_{i} = [(I − λP (π))^{−1}r(π)]_{i} ,

where i ∈ 𝕏 is an initial state and λ ∈ (0, 1) is a discount factor. One can also use the interest rate ρ = (1 − λ)/λ ∈ (0, ∞] instead of the discount factor λ. Note that λ ↑ 1 is equivalent to ρ ↓ 0. Then, the expected discount reward v^{ρ}(π) := v^{λ(ρ)}(π) can be expanded as a Laurent series ^{\underline{MV,Vt1,LdPut}} [65, 79, 59] for sufficiently small ρ

v^{ρ}(π) = (1 + ρ) [ρ^{−1}y_{−1}(π) + ∑_{m=0}^{∞} ρ^{m}y_{m}(π)] ,  (21) Laur

where y_{−1}(π) = P^{∗}(π)r(π) and and y_{m} = (−1)^{m}Y (π)^{m+1}r(π) for m ≥ 0. Note that the above series can be considered as a particular case of resolvent expansion ^{\underline{Ka,Rt}} [52, 73], which in turn can be viewed as a particular case of the inversion of singularly perturbed matrices. Note that in particular y_{−1}(π) = g(π) and y_{0} = h(π), where g(π) is the expected average reward vector defined above and h(π) is the bias vector.

Definition 1 The stationary policy π_{∗} is called a discount optimal policy for the discount factor λ ∈ (0, 1) if v_{i}^{λ}(π_{∗}) ≥ v_{i}^{λ}(π) for all i ∈ 𝕏 and all π ∈ Π^{S}.

Definition 2 The stationary policy π_{∗} is called a gain optimal policy if g_{i}(π_{∗}) ≥ g_{i}(π) for each i ∈ 𝕏 and all π ∈ Π^{S}.

It is well known that there exists a deterministic discount policy and a gain optimal policy, each of which can be computed by a number of efficient algorithms ^{\underline{Der,Kal,Put}} [26, 53, 70]. The Laurent expansion (^{\underline{Laur}} 21) allows us to define more selective optimality criteria ^{\underline{MV,Vt1,Vt2}} [65, 79, 80].

m-discount Definition 3 A policy π^{∗} ∈ Π^{S} is called an m-discount optimal policy for

some integer m ≥ 0 if

[y_{m}(π^{∗})]_{i} ≥ [y_{m}(π)]_{i}

for all i ∈ 𝕏 and all policies π that are (m − 1)-discount optimal. By convention, the gain optimal policy is (-1)-discount optimal policy.

In particular, 0-discount optimal policy is called bias optimal.

Moreover, it is known that there exists a stationary (and even deterministic) policy that is n-discount optimal for all n ≥ −1 ^{\underline{MV,Vt1,Vt2}} [65, 79, 80]. It is referred to as the Blackwell optimal policy ^{\underline{Black}} [16] and is defined equivalently in the following way.

Definition 4 A policy π_{∗} is said to be Blackwell optimal if there exists some ρ_{0} > 0 such that for all ρ ∈ (0, ρ_{0}] v^{ρ}(π^{∗}) ≥ v^{ρ}(π) for all π ∈ Π^{S}.

In other words, a Blackwell optimal policy is a policy which is discount optimal policy for any discount factor sufficiently close to one.

The results of singular perturbation theory have several applications to Markov decision processes. The first example is the Laurent series (^{\underline{Laur}} 21) for the expected discount reward. Another important application is singularly perturbed MDP. For the clarity of exposition we restrict ourselves to the case of linear perturbation, namely, we assume that the transition probabilities of the perturbed MDP are given by

p^{ε}(j|i, a) = p(j|i, a) + εd(j|i, a) ≥ 0,  (22) pert_mark

where p(j|i, a) are transition probabilities of the original unperturbed chain, ∑_{j} d(j|i, a) = 0, and ε is a “small” perturbation parameter. We are especially interested in the case of singular perturbations, that is when the perturbation changes the ergodic structure of the underlying Markov chain.

As the following example shows, the policies that are optimal for the unperturbed MDP may not coincide with the optimal policies for the perturbed MDP.

Example 2.1 Let us consider an MDP model with 𝕏 = {1, 2}, 𝔸(1) = {a_{1}, b_{1}}, 𝔸(2) = {a_{2}} and

p^{ε}(1|1, a_{1}) = 1, p^{ε}(2|1, a_{1}) = 0;
p^{ε}(1|1, b_{1}) = 1 − ε, p^{ε}(2|1, b_{1}) = ε;
p^{ε}(1|2, a_{2}) = ε, P^{ε}(2|2, a_{2}) = 1 − ε;
r(1, a_{1}) = 1, r(1, b_{1}) = 1.5, r(2, a_{2}) = 0

There are only two deterministic policies, u = [a_{1}, a_{2}] and v = [b_{1}, a_{2}]. For these policies one can easily calculate the average reward vectors. Namely,

g^{ε}(u) = \{ \begin{matrix}[ \begin{matrix}1 \\ 0 \end{matrix} ] & ε = 0 \\ [ \begin{matrix}1 \\ 0 \end{matrix} ] & ε > 0 \end{matrix} and g^{ε}(v) = \{ \begin{matrix}[ \begin{matrix}1.5 \\ 0 \end{matrix} ] & ε = 0 \\ [ \begin{matrix}0.75 \\ 0 \end{matrix} ] & ε > 0 \end{matrix}

Thus, we can see that for ε = 0 the average optimal policy is v, whereas for ε > 0 the average optimal policy is u.

Since often we do not know the exact value of the perturbation parameter ε, we are interested in finding the policy which is “close” to the optimal one for ε small but different from zero. We will call it a suboptimal policy or a limit control optimal policy. A strict definition will be given later in this section.

Most research in this topic was carried out with the assumption that the unperturbed MDP is completely decomposable and that the perturbed process is ergodic. More precisely, one can introduce the following four assumptions:

A1) 𝕏 = ∪_{k=1}^{p}𝕏_{k}, where 𝕏_{k}∩𝕏_{l} = ∅ if k ≠ l, n > 1, and ∑_{k=1}^{p} card(𝕏_{k}) = N

A2) p(j|i, a) = 0 whenever i ∈ 𝕏_{k}, j ∈ 𝕏_{l} and k ≠ l.

A3) For every i = 1, ..., p the unperturbed MDP associated with the subspace 𝕏_{k} is ergodic.

A4) The transition matrix P^{ε}(π) is irreducible for any π ∈ Π^{S} and any ε sufficiently small but different from zero, that is the perturbed MDP is ergodic.

The structure defined by assumptions A1)-A4) has a clear interpretation. The perturbed MDP model can be viewed as a complex system consisting of p “weakly-interacting” subsytems associated with 𝕏_{k}, k = 1, ..., p. Note that perturbation εd(j|i, a), where i and j are the states of different subsystems 𝕏_{k} and 𝕏_{l} respectively, represents the probability of rare transitions between the subsystems, which are independent in the unperturbed chain.

This model of singularly perturbed MDP was first studied by Pervozvanski and Gaitsgory ^{\underline{PG}} [67]. They proposed an aggregation - disaggregation algorithm for the computation of suboptimal policies for the perturbed MDP with average and discount optimality criteria. Furthermore, they have shown that in the case of discount optimality criterion, an optimal policy for the original problem is also suboptimal for the perturbed problem. The latter is not true in general if one uses the gain optimality criterion. In particular, if the perturbation is singular, the suboptimal policy of the perturbed problem can be quite different from the optimal policies of the unperturbed problem (see Example 2.1). The explanation for this phenomenon is that the interaction between weakly-coupled subsytems in the singularly perturbed process makes an effect only after sufficiently long time interval, and, if one uses discounting, the end of the trajectory has no significant contribution to the expected discount reward.

To investigate the above phenomenon, it was proposed in ^{\underline{DqQt}} [24] to consider the discount optimality criterion for the case where the interest rate goes to zero (equivalently, the discount factor goes to one) as the perturbation parameter vanishes. In particular, in ^{\underline{DqQt}} [24] it is assumed that ρ = ρ(ε) = με for some constant μ. The latter allows one to exhibit the two time scale behaviour of the singularly perturbed MDP. In addition, the model in \underline{DqQt} [24] admits a set of transient states. The dynamic programming equation for the perturbed model is solved in ^{\underline{DqQt}} [24] by using an approach based on lexicographical ordering. This concept was first applied to MDPs by Veinott \underline{Vt1} [79] to calculate a Blackwell optimal policy. Moreover, as one can seen there, singular perturbed MDPs have a close relation with Blackwell optimality \underline{Black,MV,Vt1,Put,Vt2} [16, 65, 79, 70, 80] and Markov branching decision chains (for comprehensive study of Markov branching chains the reader is refered to the paper by Huang and Veinott ^{\underline{HV}} [45]). Later, based on ^{\underline{D}}[25], results which appeared first in ^{\underline{DqQt}} [24], were extended in ^{\underline{Q}}[71] to include the possibility of several time scales. The key point in ^{\underline{D}}[25] and ^{\underline{Q}}[71] is the use of the following dependency of the interest rate on the perturbation parameter: ρ = ρ(ε) = με^{l}, where l represents the order of a time scale. In ^{\underline{Al}} [4] Altman and Gaitsgory have generalized the results of ^{\underline{PG}} [67] to constrained Markov decision processes.

Abbad, Bielecki and Filar ^{\underline{AbFl,AbFlBl,BlFl,AFs}} [1, 3, 14, 2] formulated the limit control principle, which provides a formal setting for the concept of suboptimal policies. First, as elaborated in Section 2, we note that for any stationary policy π ∈ Π^{S}, there exists a limiting ergodic projection

P̂^{∗}(π) := lim_{ε→0} P^{∗,ε}(π).

The average reward optimization problem for the perturbed MDP can be written in the form

g_{i}^{opt,ε} = max_{π∈Π^{S}}[P^{∗,ε}(π)r(π)]_{i} (L^{ε}).

The limit control principle says that instead of the above singular program one can consider a well-defined limit Markov control problem:

ĝ_{i}^{opt} = max_{π∈Π^{S}}[P̂^{∗}(π)r(π)]_{i} (L).

It is natural to expect that an optimal strategy to (L), if exists, will approximate well the optimal strategy for the perturbed problem (L^{ε}), in the case where the perturbation parameter is small. In ^{\underline{BlFl}} [14] it was shown that this is indeed the case. Specifically, let π_{∗} be any maximizer in (L), then

lim_{ε→0} max_{i∈𝕏}|g_{i}^{ε}(π_{∗}) − g^{opt,ε}| = 0.

In ^{\underline{AbFl}} [1] it is proved that there exists a deterministic policy that solves the limit control problem (L). Recently, Bielecki and Stettner ^{\underline{BlSt}} [15] have generalized the limit control principle to MDPs with general Borel state spaces and compact action spaces.

Under the assumptions A1) – A4) the limit Markov control problem (L) can be solved by solving the following linear programming problem (P):

maximise ∑_{k=1}^{n} ∑_{i∈𝕏_{k}} ∑_{a∈A(i)} r(i, a)z_{ia}^{k}

subject to:

∑_{i∈𝕏_{k}} ∑_{a∈A(i)}(δ_{ij} − p(j|i, a))z_{ia}^{k} = 0, j ∈ 𝕏_{k}; k = 1, …, n

∑_{k=1}^{n} ∑_{j∈𝕏_{ℓ}} ∑_{i∈𝕏_{k}} ∑_{a∈A(i)} d(j|i, a)z_{ia}^{k} = 0; ℓ = 1, …, n

∑_{k=1}^{n} ∑_{i∈𝕏_{k}} ∑_{a∈A(i)} z_{ia}^{k} = 1

z_{ia}^{k} ≥ 0 k = 1, …, n; i ∈ 𝕏_{k}, a ∈ A(i).

It can be shown (see ^{\underline{AbFlBl}} [3]) that an optimal strategy in the limit Markov control problem (L) can be constructed as follows.

Theorem 8 Let {z_{ia}^{k}|k = 1, …, n; i ∈ 𝕏_{k}; a ∈ A(i)} be an optimal extreme solution to the linear program (P), then the deterministic strategy defined by

f_{∗}(i) = a, i ∈ 𝕏_{k}, k = 1, …, n ⟺ z_{ia}^{k} > 0

is optimal in the limit Markov control problem (L).

The linear program (P) is similar to one given by Gaitsgory and Pervozvanski ^{\underline{GP2}} [33]. However, these authors used techniques different from those in

\underline{AbFlBl} [3].

As a result of the solution of problem (L), one gets limit control optimal policy f_{∗} ∈ Π^{D} such that for any other deterministic policy f ∈ Π^{D} the holds

lim_{ε→0}(g_{i}^{ε}(f_{∗}) − g_{i}^{ε}(f)) ≥ 0, i ∈ 𝕏.

However, as was noted in ^{\underline{AbFl}} [1], a policy that solves (L), in general, is only suboptimal. Interestingly, there exists a policy that is optimal for all sufficiently small ε but different from zero (see ^{\underline{AbFl}} [1]). This policy is called uniform optimal ^{\underline{AlAvFl}} [5] and it satisfies the following inequality

g_{i}^{ε}(f_{∗}) ≥ g_{i}^{ε}(f), i ∈ 𝕏,

for any deterministic policy f and all ε sufficiently small but different from zero. We like to emphasize that a uniform optimal policy is limit control optimal, but the converse need not hold, as the following example illustrates.

Example 2.2 Consider 𝕏 = {1, 2}, 𝔸(1) = {a_{1}, b_{1}}, 𝔸(2) = {a_{2}}; let

unif_examp

p^{ε}(1|1, a_{1}) = 1,  r(1, a_{1}) = 10
p^{ε}(1|1, b_{1}) = 1 − ε, r(1, b_{1}) = 10
p^{ε}(2|1, b_{1}) = ε,
p^{ε}(1|2, a_{2}) = 1,  r(2, a_{2}) = 0

Then the stationary policy u(1) = a_{1}, u(2) = a_{2} is uniformly optimal with expected average reward g_{i}^{ε}(u) = 10. The stationary policy v(1) = b_{1}, v(2) = a_{2} is limit control optimal as lim_{ε→0} g_{i}^{ε}(v) = 10, but for every ε > 0,

g_{i}^{ε}(v) = \frac{10}{1 + ε} < g_{i}^{ε}(u).

A heuristic procedure to determine the uniform optimal policy is proposed in ^{\underline{AbFl}} [1]. In contrast, a well-defined computational procedure was proposed in ^{\underline{AlAvFl}} [5], which is based on asymptotic linear programming. The asymptotic linear programming was first introduced by Jeroslow ^{\underline{Jer1,Jer2}} [50, 51] and later refined by Hordijk, Dekker and Kallenberg ^{\underline{HDK}} [43]. It is designed to solve the following parametric linear program

max c(ε)x  (23) 1.3
A(ε)x = b(ε), x ≥ 0  (24) 1.4

for ε small but different from zero. The elements of A(ε), b(ε) and c(ε) are assumed to be polynomials of ε. The basic idea of this method is to perform the simplex algorithm over the field of rational functions instead of the field of real numbers. The latter results in a solution that is optimal for all ε sufficiently small but different from zero. We would like to note that instead of asymptotic linear programming one can use its modification, the asymptotic simplex method ^{\underline{PG,Lam1,Lam2,AlFlAv}} [67, 57, 58, 7], which operates over the field of Laurent series. Moreover, if the asymptotic linear programming or asymptotic simplex method is used, then one can drop Assumptions A1)-A4) and consider the most general situation.

As can be concluded from the above overview of Markov decision processes, there is an intrinsic relationship between the Blackwell optimality and singular perturbations. Actually, the contributions of Delebecque and Quadrat ^{\underline{DqQt,D,Q}} [24, 25, 71] clearly emphasize this connection. It turns out that the Blackwell optimality, singularly perturbed MDPs and branching Markov decision chains can all be considered in a unified framework based on asymptotic simplex method. Towards this end, let us consider the linear programming formulation for the following four MDP models.

Model I: MDP with Blackwell optimality criterion

A Blackwell optimal policy is a discount optimal policy for all discount factor sufficiently close to one, or, equivalently, the policy which is optimal for all interest rate sufficiently close to zero. A Blackwell optimal policy can be determined by using the asymptotic linear programming approach \underline{AHK,HDK,Lam1,Lam2} [6, 43, 57, 58]. Namely, one needs to find a solution for the following linear program which is optimal for all interest rate ρ > 0 which are sufficiently small:

max ∑_{i,a}(1 + ρ)r(i, a)x_{ia}

∑_{i,a}[(1 + ρ)δ_{ij} − p(j|i, a)]x_{ia} = 1, x_{ia} ≥ 0.  (25) Bl_model

Note that the above linear program can be immediately written in the form (^{\underline{1.3}} 23),(^{\underline{1.4}} 24) with ε = ρ. In ^{\underline{HDK}} [43] a refined version of Jeroslow’s asymptotic linear programming method ^{\underline{Jer1,Jer2}} [50, 51] for solving (^{\underline{Bl_model}} 25) was developed. Later Lamond ^{\underline{Lam2}} [58] applied his version of the asymtotic simplex method to the above parametric LP. Of course, our algorithm ^{\underline{AlFlAv}} [7] can be applied to this model as well. As pointed out in ^{\underline{Lam2}} [58], in this particular case, the Laurent series for the basis matrix always admits a simple pole. This in turn implies that the complexity of the basis updating is O(n^{2}), which is comparable with the complexity of the basis updating for the ordinary simplex method.

Model II: Markov branching decision chains.

The Markov branching decision chains were introduced in ^{\underline{HV}} [45]. These are MDPs with immediate rewards which dependent on the interest rate. Namely, it is assumed that r(i, a) = r^{ρ}(i, a) is a polynomial in the interest rate ρ. To find a policy which is optimal for all sufficiently small ρ (Strongvalue policy ^{\underline{HV}} [45]), we need to solve only slightly modified version of (^{\underline{Bl_model}} 25), that is

max ∑_{i,a}(1 + ρ)r^{ρ}(i, a)x_{ia}

∑_{i,a}[(1 + ρ)δ_{ij} − p(j|i, a)]x_{ia} = 1, x_{ia} ≥ 0.  (26) Br_model

Again the asymptotic simplex method can be immediately applied. A method of finding the optimal policy proposed by Huang and Veinott ^{\underline{HV}} [45] also employs the asymptotic programming approach. The method is based on the subsequent solution of the augmented LPs, whose objective functions and right hand sides (but not the constraint coefficient matrix) depend on the parameter.

Model III: Singularly perturbed MDP with average criterion

In Section 3.4 we applied the asymptotic linear programming ^{\underline{Jer1,Jer2,HDK}} [50, 51, 43], which is based on the ordering in the field of rational functions, to singularly perturbed MDPs with average criterion. Since one needs to solve the following parametric program

max ∑_{i,a} r(i, a)x_{ia}

∑_{i,a}[δ_{ij} − p(j|i, a) − εd(j|i, a)]x_{ia} = 0,

∑_{a} x_{ja} + ∑_{i,a}[(δ_{ij} − p(j|i, a) − εd(j|i, a)]y_{ia} = β_{j},

x_{ia} ≥ 0,  y_{ia} ≥ 0,  (27) SP_model

the asymptotic simplex method, based on the Laurent series expansions, can be applied as well. Moreover, since the basis updating in the asymptotic linear programming takes O(m^{4}log(m)) flops, and the asymptotic simplex method takes s̄m^{2} flops. The latter method is computationally superior.

Model IV: Singularly perturbed MDP with killing interest rate

In ^{\underline{DqQt,Q}} [24, 71] considered singularly perturbed MDPs with “killing interest rate” ρ(ε) = με^{l}, where l is the order of a time scale were considered. This model exhibits the necessity of different control regimes for the different time scales. The papers ^{\underline{DqQt,Q}} [24, 71] provided a lexicographical policy improvement algorithm for the solution of the perturbed dynamic programming equation.

Alternatively, the extension of our asymptotic simplex method for the polynomial perturbation can be used in this problem. Here the parametric linear program takes the following form

max ∑_{i,a}(1 + με^{l})r(i, a)x_{ia}

∑_{i,a}[(1 + με^{l})δ_{ij} − p(j|i, a) − εd(j|i, a)]x_{ia} = 1, x_{ia} ≥ 0.  (28) QD_model

Generalized Model:

Finally, we would like to note that the Models I,II and IV can be viewed as particular cases of the next unified scheme. Let transition probabilities p^{ε}(j|i, a), immediate rewards r^{ε}(i, a) and interest rate ρ(ε) of an MDP model be polynomials of the parameter ε. Then a policy which is optimal for all sufficiently small values of parameter ε can be found from the next perturbed linear program.

max ∑_{i,a}(1 + ρ(ε))r^{ε}(i, a)x_{ia}

∑_{i,a}[(1 + ρ(ε))δ_{ij} − p^{ε}(j|i, a)]x_{ia} = 1, x_{ia} ≥ 0.  (29) gen_model

The above perturbed LP can be efficiently solved by the generalization of the asymptotic simplex method proposed in ^{\underline{AlFlAv}} [7]. Note that we retrieve Model I with ρ(ε) = ε, r^{ε}(i, a) = r(i, a), p^{ε}(j|i, a) = p(j|i, a); Model II with ρ(ε) = ε, r^{ε}(i, a) = ∑_{k=0}^{p} ε^{k}r_{k}(i, a), p^{ε}(j|i, a) = p(j|i, a); and Model IV with ρ(ε) = με^{l}, r^{ε}(i, a) = r(i, a), p^{ε}(j|i, a) = p(j|i, a) + εd(j|i, a).

4 Singularly Perturbed MDP’s and the Combinatorial Hamiltonian Cycle Problem

In this section we demonstrate that a famous combinatorial optimisation problem such as the Hamiltonian Cycle Problem (HCP) can be regarded as a singularly perturbed Markov Decision Process. We begin with a brief description of only one version of the HCP.

In graph theoretic terms, the problem is to find a simple cycle of N arcs, that is a Hamiltonian Cycle or a tour, in a directed graph G with N nodes and with arcs (i, j), or to determine that none exist. Recall that a simple cycle is one that passes exactly once through each node comprising the cycle.

It is known that HCP belongs to the NP-complete class of problems and, as such, is considered very difficult from an algorithmic perspective.

In this section we propose the following, unorthodox, perspective of the Hamiltonian cycle problem: Consider a moving object tracing out a directed path on the graph G with its movement “controlled” by a function f mapping the set of nodes 𝕏 = {1, 2, …, N} into the set of arcs A.

Clearly, we can think of this set of nodes as the state space of a Markov decision process Γ where for each state/node i, the action space

A(i) = {a = j|(i, j) ∈ A}

is in one to one correspondence with the set of arcs emanating from that node.

Of course, we shall ignore the trivial case A(i) = ∅, because in such a case, obviously, no Hamiltonian cycle exists. Furthermore, if we restrict the function f above in such a way that f(i) ∈ A(i), for each i ∈ 𝕏, then we see that f can be thought of as a deterministic strategy f in an MDP Γ. Designating node 1 as the “home node” in G, we shall say that f is a Hamiltonian cycle in G if the set of arcs {(1, f(1)), (2, f(2)), …, (N, f(N))} is a Hamiltonian cycle in G. If the above set of arcs contains cycles of length less than N, we shall say that f has subcycles in G.

Note that if P (f) is the transition probability matrix of a Markov chain induced by f that is a Hamiltonian cycle, then P (f) is irreducible and the long-run frequency of visits to any state x_{i}(f) = 1/N. Of course, if f has subcycles in G, then P (f) contains multiple ergodic classes which complicates the analysis of the Markov decision process Γ, in which we have embedded our graph theoretic problem.

A class of limiting average Markov decision processes that retains most of the desirable properties of the irreducible processes is the so-called “unichained” class. Briefly, a Markov decision process is unichained if for every deterministic stationary control f, P (f) contains only a single ergodic class and possibly a nonempty set of transient states. We now perturb the transition probabilities of Γ slightly to create an ε-perturbed process Γ(ε) (for 0 < ε < 1) defined by:

p^{ε}(j|i, a) = \{ \begin{matrix}1 & if i = 1 and a = j \\ 0 & if i = 1 and a ≠ j \\ 1 & if i > 1 and a = j = 1 \\ ε & if i > 1, a ≠ j, and j = 1 \\ 1 − ε & if i > 1, a = j, and j > 1 \\ 0 & if i > 1, a ≠ j, and j > 1. \end{matrix}

Note that with the above perturbation, for each pair of nodes i, j (not equal to 1) corresponding to a “deterministic arc” (i, j) our perturbation replaces that arc by a pair of “stochastic arcs” (i, 1) and (i, j):

Figure 1:

with weights ε and (1 − ε) respectively (ε ∈ (0, 1)). Note that this perturbation changes Γ to an ε-perturbed Markov decision process Γ(ε).

Example 4.1 Consider the following complete graph G on four nodes (with no self-loops):

Figure 2:

and think of the nodes as the states of an MDP, denoted by Γ, and of the arcs emanating from a given node as actions available at that state. The Hamiltonian cycle c_{1} : 1 → 2 → 3 → 4 → 1 corresponds to the deterministic stationary strategy f_{1} : {1, 2, 3, 4} → {2, 3, 4, 1} where f_{1}(2) = 3 corresponds to the controller choosing arc (2,3) in state 2 with probability 1. The Markov chain induced by f_{1} is given by the transition matrix

P (f_{1}) = ( \begin{matrix}0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 \end{matrix} )

which is irreducible. On the other hand, the union of two sub-cycles: 1 → 2 → 1 and 3 → 4 → 3 corresponds to the policy f_{2} : {1, 2, 3, 4} → {2, 1, 4, 3} which identifies the Markov chain transition matrix

P (f_{2}) = ( \begin{matrix}0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{matrix} )

containing two distinct ergodic classes.

As mentioned earlier, the perturbation destroys multiple ergodic classes and induces a unichained, singularly perturbed, Markov decision process Γ(ε). For instance, the policy f_{2} now has the Markov chain matrix

P (f_{2}) = ( \begin{matrix}0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ ε & 0 & 0 & 1 − ε \\ ε & 0 & 1 − ε & 0 \end{matrix} ) .

Remark 4.1 Our perturbation has made the home node/state 1 rather special. In particular, the home state always belongs to the single ergodic class of P (f) for any f ∈ Π^{S}. Of course, some other states could be transient.

We shall undertake the analysis of the Hamiltonian cycle problem in the “frequency space” of the perturbed process Γ(ε). Recall that with every f ∈ Π^{S} we can associate the long-run frequency vector x(f). This is achieved by defining a map M : Π^{S} → X(ε) by

x_{ia}(f) = π_{i}(ε, f)f(a|i); f ∈ Π^{S}

for each i ∈ 𝕏 and a ∈ A(i), where π_{i}(ε, f) is the i-th element of the stationary distribution vector of the perturbed Markov chain transitions matrix P (f), and f(a|i) is the probability of choosing action a in state i.

Now consider the polyhedral set X(ε) defined by the constraints

(i) ∑_{i=1}^{N} ∑_{a∈A(i)} [δ(i, j) − p_{ε}(j|i, a)] x_{ia} = 0; j ∈ 𝕏.

(ii) ∑_{i=1}^{N} ∑_{a∈A(i)}x_{ia} = 1.

(iii) x_{ia} ≥ 0; a ∈ A(i), i ∈ 𝕏.

Next define a map M̂ : X(ε) → Π^{S} by

f_{x}(i, a) = \{ \begin{matrix}\frac{x_{ia}}{x_{i}}; & if x_{i} = ∑_{a∈A(i)} x_{ia} > 0 \\ 1; & if x_{i} = 0 and a = a_{1} \\ 0; & if x_{i} = 0 and a ≠ a_{1}, \end{matrix}

for every a ∈ A(i), i ∈ 𝕏 where a_{1} denotes the first available action in a given state according to some ordering. The following result can be found in ^{\underline{FK}} [28] and ^{\underline{FiVr}} [31].

Lemma 4.1

(i) The set X(ε) = {x(f)|f ∈ Π^{S}} and will henceforth be called the (longrun) “frequency space” of Γ(ε).

(ii) For every x ∈ X(ε),

M (M̂(x)) = x

but the inverse of M need not exist.

(iii) If x is an extreme point of X(ε), then

f_{x} = M̂(x) ∈ Π^{D}.

(iv) If f ∈ Π^{D} is a Hamiltonian cycle, then x(f) is an extreme point of X(ε).

We shall now derive a useful partition of the class Π^{D} of deterministic strategies that is based on the graphs they “trace out” in G. In particular, note that with each f ∈ Π^{D} we can associate a subgraph G_{f} of G defined by

arc (i, j) ∈ G_{f} ⟺ f(i) = j

We shall also denote a simple cycle of length m and beginning at 1 by a set of arcs

c_{m}^{1} = {(i_{1} = 1, i_{2}), (i_{2}, i_{3}), …, (i_{m}, i_{m+1} = 1)}; m = 2, 3, …, N.

Of course, c_{N}^{1} is a Hamiltonian cycle. If G_{f} contains a cycle c_{m}^{1} we write G_{f} ⊃ c_{m}^{1}. For 2 ≤ m ≤ N, let C_{m} := {f ∈ Π^{D}|G_{f} ⊃ c_{m}^{1}},namely, the set of deterministic strategies that trace out a simple cycle of length m, beginning at 1, for each m = 2, 3, …, N. Of course, C_{N} is the set of strategies that correspond to Hamiltonian cycles and any single C_{m} can be empty, depending on the structure of the original graph G. Thus a typical strategy f ∈ C_{3}, for example, traces out a graph G_{f} in G that might look like Figure 3 where the dots indicate the “immaterial” remainder of G_{f} that corresponds to states that are transient in P (f).

Figure 3:

The partition of the deterministic strategies that seems to be most relevant for our purposes is

Π^{D} = [ ⋃_{m=2}^{N} C_{m}] ⋃ B,  (30) eq2.5.1

where B contains^{6} all the deterministic strategies that are not in any of the C_{m}’s. Note that a typical strategy f in B traces out a graph G_{f} in G that might looks as follows:

6It will soon be seen that the strategies in B are in a certain sense “bad” or, more precisely, difficult to analyse, thereby motivating the symbol B.

Figure 4:

where the dots again denote the immaterial part of G_{f}. However, it is important to note that for any ε > 0, the states 1, i_{2}, …, i_{k−1} are not transient in Γ_{α}(ε).

It is, perhaps, interesting to observe that all strategies in a given set in the partition (^{\underline{eq2.5.1}} 30) induce the same long-run frequency x_{1}(f) of visits to the home node 1. This observation is captured in the following proposition which can be found in ^{\underline{MF}} [66] and ^{\underline{FiVr}} [31]

Proposition 4.2 Let ε ∈ (0, 1), f ∈ Π^{D}, and x(f) be its long-run frequency vector (that is, x(f) = M(f)). The long-run frequency of visits to the home state 1 is given by

x_{1}(f) = ∑_{a∈A(1)} x_{1a}(f) = \{ \begin{matrix}\frac{1}{d_{m}(ε)}; & if f ∈ C_{m}, m = 2, 3, …, N \\ \frac{ε}{1 + ε}; & if f ∈ B, \end{matrix}

where d_{m}(ε) = 1 + ∑_{i=2}^{m}(1 − ε)^{i−2} for m = 2, 3, …, N.

The above proposition leads the following characterizations of the Hamiltonian cycles of a directed graph.

Theorem 9

(i) Let f ∈ Π^{D} be a Hamiltonian cycle in the graph G. Then G_{f} = c_{N}^{1}, x(f) is an extreme point of X(ε) and x_{1}(f) = \frac{1}{d_{N} (ε)}.

(ii) Conversely, suppose that x is an extreme point of X(ε) and that x_{1} = ∑_{a∈A(1)} x_{1a} = \frac{1}{d_{N} (ε)}, then f = M̂(x) is a Hamiltonian cycle in G.

(iii) Hamiltonian cycles of the graph G are in 1 : 1 correspondence with those points of X(ε) which satisfy

(a) x_{1} = ∑_{a∈A(1)}x_{1a} = \frac{1}{d_{N}(ε)}.

(b) For every i ∈ 𝕏, x_{i} = ∑_{a∈A(1)}x_{ia} > 0 and \frac{x_{ia}}{x_{i}} ∈ \{0, 1\} for each

a ∈ A(i), i ∈ 𝕏.

Remark 4.2: It is, perhaps, significant to note that for all ε ∈ (0, 1), m = 2, 3, …, N − 1

\frac{1}{d_{m}(ε)} > \frac{1}{d_{m+1}(ε)} > \frac{ε}{1 + ε}.

Thus Theorem 9 demonstrates that the extreme points x of X(ε) can be “ranked” according to their values of the linear function l(x) = ∑_{a∈A(1)} x_{1a}. Unfortunately, the Hamiltonian cycles (if they exist) may attain only the “second lowest” value of l(x), namely, \frac{1}{d_{N} (ε)}. The latter problem is, partially, rectified in the following result that can be found in ^{\underline{FiLi}} [29].

Theorem 10 Let f^{∗} ∈ Π^{D} be a Hamiltonian cycle in the graph G. Then for ε ≥ 0 and sufficiently small f^{∗} is a global minimiser of the following perturbed optimisation problem:

min_{f∈Π^{D}}{[I − P (f) + P^{∗}(f)]_{11}^{−1}},

where H_{11} denotes the (1,1)^{th} entry of matrix H.

The above theorem shows that if Hamiltonian cycles exist, they are the minimisers of the top left element of the fundamental matrix, over f ∈ Π^{D}, as long as ε is sufficiently near 0. However, from the optimisation point of view, elements of the fundamental matrix are not straightforward to analyse. This leads naturally to the open problem: can asymptotic expansions such as the one in Theorem 5, facilitate algorithmic approaches to the optimisation problem in Theorem 10?

To date the best algorithmic results based on this stochastic approach to the HCP are reported in ^{\underline{And}} [8]. These results exploit the properties in (i)

– (iii) of Theorem 9, above. In particular, it can be checked that the most awkward requirement x_{ia}/x_{i} ∈ {0, 1} for all i ∈ 𝕏, a ∈ 𝔸(i) is equivalent min{x_{ia}, x_{ib}} = 0 for all i ∈ 𝕏, a, b ∈ 𝔸(i) and a ≠ b. This observation immediately leads to the following mixed integer programming formulation of the HCP problem:

min ∑_{i} ∑_{a} c_{ia}x_{ia}
s.t.  x ∈ X(ε)
x_{1} = 1/d_{N}(ε)
x_{ia} ≤ My_{ia}  : i ∈ 𝕏, a ∈ 𝔸(i)
y_{ia} + y_{ib} ≤ 1  ; i ∈ 𝕏, a, b ∈ 𝔸(i), a ≠ b
y_{ia} ∈ {0, 1} ; i ∈ 𝕏, a ∈ 𝔸(i).

In the above M ≥ 1/d_{N}(ε) and c_{ia}’s can be experimented with. In preliminary numerical experiments reported in ^{\underline{And}} [8], randomly generated problems with up to 100 nodes and 300 arcs were solved in less than 150 cpu seconds on a Sun Workstation. Recently, Filar and Lasserre ^{\underline{FL}} [30] proposed a non-standard branch and bound algorithm for the formulation

min ∑_{i} ∑_{a} c_{ia}x_{ia}

s.t. (a) x ∈ X(ε)

(b) x_{1} = 1/d_{N}(ε)

(c) x_{ia}/x_{i} ∈ {0, 1}; i ∈ 𝕏, a ∈ 𝔸(i)

which exploits the structural property that ensures that if the relaxation omitting (c) is solved by the simplex method, then (c) can be violated for at most one i and at most one pair of arcs a, b.

In an interesting, related, development Feinberg ^{\underline{Fein}} [27] has recently considered the embedding of a Hamiltonian cycle problem in a discounted Markov decision process. In that paper the perturbation parameter ε is not necessary but, instead, the discount factor λ ∈ [0, 1) plays a crucial role. In particular Feinberg’s embedding can be obtained by setting ε = 0 in p_{1}^{ε}(j|i, a)’s defined earlier and by setting

r(i, a) = \{ \begin{matrix}1 & if i = 1, a ∈ 𝔸(1) \\ 0 & otherwise \end{matrix}

For any π ∈ Π^{S} the expected discounted reward v_{i}^{λ}(π) is now defined as in Section 3. The analysis in ^{\underline{Fein}} [27] is based on the following interesting observation. Let i_{m} denote the state/node visited at stage m, then

v_{1}^{λ}(π) = ∑_{m=0}^{∞} λ^{m}P_{1}^{π}(i_{m} = 1),

where P_{1}^{π}(·) denotes the probability measure induced by π and the initial state i_{0} = 1, and

P_{1}^{π}(i_{m} = 1) = \frac{1}{m!} [\frac{∂^{m}}{∂λ^{m}} (v_{1}^{λ}(π))]_{λ=0} .

The above lead to novel characterisations of Hamiltonian cycles that are summarised below.

Theorem 11 With the embedding in Γ_{λ} described above the following statements are equivalent:

(i) A policy π = f is deterministic and a Hamiltonian cycle in G.

(ii) A policy π is stationary and a Hamiltonian cycle in G.

(iii) A policy f is deterministic and v_{1}^{λ}(f) = (1 − λ^{N})^{−1} for at least λ ∈ [0, 1).

(iv) A policy π is stationary and v_{1}^{λ}(π) = (1 − λ^{N})^{−1} for 2N − 1 distinct discount factors λ_{k} ∈ (0, 1); k = 1, 2, …, 2N − 1.

The above characterisation naturally leads to a number of mathematical programming formulations of both HCP and TSP that are described in \underline{Fein} [27]. There is clearly a need to explore the algorithmic potential of these formulations.

References

AbFl [1] M. Abbad and J.A. Filar, “Perturbation and stability theory for Markov control problems”, IEEE Trans. Auto. Contr., AC-37, no.9, pp.1415-1420, 1992.

AFs [2] M. Abbad and J.A. Filar, “Algorithms for singularly perturbed Markov control problems: A survey”, in Techniques in discrete-time stochastic control systems (ed.) C.T. Leondes, Series: Control and Dynamic Systems, v.73, Academic Press, New York, 1995.

AbFlBl [3] M. Abbad, J.A. Filar and T.R. Bielecki, Algorithms for singularly perturbed limiting average Markov control problems. IEEE Trans. on Automatic Control 37:1421-1425, 1992.

Al [4] E. Altman and V.G. Gaitsgory, “Stability and Singular Perturbations in Constrained Markov Decision Problems”, IEEE Trans. Autom. Control, v.38, pp.971-975, 1993.

AlAvFl [5] E. Altman, K.E. Avrachenkov, and J.A. Filar, “Asymptotic linear programming and policy improvement for singularly perturbed Markov decision processes”, ZOR: Math. Meth. Oper. Res., v.49, pp.97-109, 1999.

AHK [6] E. Altman, A. Hordijk, and L.C.M. Kallenberg, “On the value function in constrained control of Markov chains”, ZOR: Math. Meth. Oper. Res., v.44, pp.387-399, 1996.

AlFlAv [7] J.A. Filar, E. Altman and K.E. Avrachenkov, “An asymptotic simplex method for singularly perturbed linear programs”, submited to Operations Research Letters, 1999.

And [8] M. Andramonov, J. Filar, A. Rubinov and P. Pardalos, “Hamiltonian Cycle Problem via Markov Chains and Min-type Approaches”in Approximation and Complexity in Numerical Optimization: Continuous and Discrete Problems, Ed. P.M. Pardalos, Kluwer Aacademic Publishers, to appear.

A [9] K.E. Avrachenkov, “Analytic perturbation theory and its applications”, PhD Thesis, University of South Australia, 1999.

AH [10] K.E. Avrachenkov and M. Haviv, “The highest singular coefficients in singular perturbation of stochastic marices,” (to appear).

AHH [11] K.E. Avrachenkov, M. Haviv and P.G. Howlett, “Inversion of analytic matrix functions that are singular at the origin”, submitted to SIAM J. Matr. Anal. Appl.

AL [12] K.E. Avrachenkov and J.B. Lasserre, “The fundamental matrix of singularly perturbed Markov chains,” Advances in Applied Probability, v.31, (to appear).

AL2 [13] K.E. Avrachenkov and J.B. Lasserre, “Perturbation analysis of reduced resolvents and generalized inverses”, LAAS Report No. 98520, also submited to Numerische Mathematik, 1998.

BlFl [14] T.R. Bielecki and J.A. Filar, “Singularly perturbed Markov control problem: Limiting average cost”, Annals O.R., v.28, pp.153-168, 1991.

BlSt [15] T.R. Bielecki and L. Stettner, “Ergodic control of singularly perturbed Markov process in discrete time with general state and compact action spaces”, submitted, 1996.

Black [16] D. Blackwell, “Discrete dynamic programming”, Ann. Math. Stat., v.33, pp.719-726, 1962.

CW [17] M. Cordech, A.S. Willsky, S.S. Sastry and D.A. Castanon, Hierarchical aggregation of linear systems with multiple time scales, IEEE Trans. Autom. Contr., AC-28, (1983) 1029-1071.

CdWl [18] M. Cordech, A.S. Willsky, S.S. Sastry, and D.A. Castanon, “Hierarchical aggregation of singularly perturbed finite state Markov processes,” Stochastics. v.8, pp.259-289, 1983.

Co [19] P.J. Courtois, Decomposability: queueing and computer system applications, Academic Press, New York, 1977.

CL [20] P.J. Courtois and G. Louchard, “Approximation of eigencharacteristics in nearly-completely decomposable stochastic systems”, Stoch. Process. Appl., v.4, pp.283-296, 1976.

CS [21] P.J. Courtois and P. Semel, “Bounds for the positive eigenvectors of non-negative matrices and their approximation by decomposition”, JACM, v.31, pp.804-825, 1984.

Dek [22] R. Dekker, Denumerable Markov decision chains: optimal policies for small interest rate, Ph.D. Thesis, Univ. of Leiden.

DH [23] R. Dekker and A. Hordijk, “Average, sensitive and Blackwell optimal policies in denumerable Markov decision chains with unbounded rewards”, Math. Oper. Res., v.13, pp.395-420, 1988.

DqQt [24] F. Delebecque and J.P. Quadrat, “Optimal control of Markov chains admitting strong and weak interactions”, Automatica, v.17, pp.281-296, 1981.

D [25] F. Delebecque, A reduction process for perturbed Markov chain, SIAM J. Appl. Math., 43, (1983) 325-350.

Der [26] C. Derman, Finite state Markovian decision processes, Academic Press, New York, 1970.

Fein [27] E.A. Feinberg, “Constrained Discounted Markov Decision Processes and Hamiltonian Cycles”, MOR, to appear.

FK [28] J.A. Filar and D. Krass, “Hamiltion cycles and Markov chains,” Mathematics of Operations Research, Vol. 19, pp. 223–237, 1994.

FiLi [29] J.A. Filar and Ke Liu, “Hamiltonian Cycle Problem and Singularly Perturbed Decision Process”, in Statistics, Probability and Game Theory: Papers in Honor of David Blackwell, IMS Lecture Notes – Monograph Series, USA, 1996.

FL [30] J.A. Filar and J-B Lasserre, “A Non-Standard Branch and Bound Method for the Hamiltonian Cycle Problem”, LAAS Report 98247, June 1998.

FiVr [31] J.A. Filar and K. Vrieze, Competitive Markov Decision Processes, Springer-Verlag, N.Y., 1996.

GP [32] V.G. Gaitsgori and A.A. Pervozvanskii, “Aggregation of states in a Markov chain with weak interactions”, Cybernetics, v.11, pp.441-450, 1975. (Translation of Russian original in Kibernetika, v.11, pp.91-98, 1975.)

GP2 [33] V.G. Gaitsgori and A.A. Pervozvanskii, Theory of Suboptimal Decisions, Kluwer Academic Publishers, 1988.

HH [34] R. Hassin, and M. Haviv, “Mean passage times and nearly uncoupled Markov chains,” SAIM Journal of Discret Mathematics, Vol. 5, pp. 386–397, 1992.

Ha2 [35] M. Haviv, “An approximation to the stationary distribution of a nearly completely decomposable Markov chain and its error analysis,” SIAM Journal on Algebraic and Discrete Methods, vol. 7, pp. 589–594, 1986.

Ha3 [36] M. Haviv, “Aggregation/disaggregation methods for computing the stationary distribution of a Markov chain,” SIAM Journal on Numerical Analysis, vol. 24, pp. 952–966, 1987.

Ha1 [37] M, Haviv, “More on the Rayliegh-Ritz refinement technique for nearly uncoupled matrices,” SIAM Journal of Matrix Analysis and Application, Vol. 10, pp. 287–293, 1989.

[38] M. Haviv and Y. Ritov, “An approximation to the stationary distribution of a nearly completely decomposable Markov chain and its error bounds,” SIAM Journal on Algebraic and Discrete Methods, vol. 7, pp. 583–588, 1986.

HR1 [39] M. Haviv and Y. Ritov, “Series expansions for stochastic matrices,” unpublished manuscript, 1989.

HR2 [40] M. Haviv and Y. Ritov, “On series expansions for stochastic matrices,” SIAM Journal on Matrix Analysis and Applications, vol. 14, pp. 670– 677, 1993.

HV [41] M. Haviv and L. Van der Heyden, “Perturbation bounds for the stationary probabilities of a finite Markov chain,” Advances in Applied Probability, vol. 16, pp. 804–818, 1984.

HL [42] O. Hernandez-Lerma and J.B. Lasserre, Discrete- time Markov control processes: basic optimality criteria, Springer- Verlag, New York, 1996.

HDK [43] A. Hordijk, R. Dekker, and L.C.M. Kallenberg, “Sensitivity analysis in discounted Markovian decision problems”, Spectrum, v.7, pp.143-151, 1985.

How [44] R.A. Howard, Dynamic programming and Markov processes, Cambridge, MA: MIT Press, 1960.

HV [45] Y. Huang and A.F. Veinott, Jr., “Markov branching decision chains with interest-rate-dependent rewards”, Probability in the Engineering and Information Sciences, v.9, pp.99-121, 1995.

Hu [46] Y. Huang, “A canonical form for pencils of matrices with applications to asymptotic linear programs, Lin. Alg. Appl., v.234, pp.97-123, 1996.

Hu1 [47] J.J. Hunter, “Stationary distributions of perturbed Markov chains”, Lin. Alg. Appl., v.82, pp.201-214, 1986.

Hu2 [48] J.J. Hunter, “The computation of stationary distributions of Markov chains through perturbations”, J. Appl. Math. Stoch. Anal., v.4, pp.29-46, 1991.

Hu3 [49] J.J. Hunter, “A survey of generalized inverses and their use in applied probability”, Math. Chronicle, v.20, pp.13-26, 1991.

Jer1 [50] R.G. Jeroslow, “Asymptotic Linear Programming”, Oper. Res., v.21, pp.1128-1141, 1973.

Jer2 [51] R.G. Jeroslow, “Linear Programs Dependent on a Single Parameter”, Disc. Math., v.6, pp.119-140, 1973.

Ka [52] T. Kato, Perturbation theory for linear operators, Springer-Verlag, Berlin, 1966.

Kal [53] L. C. M. Kallenberg, Linear programming and finite Markovian control problems, Mathematical Centre Tracts 148, Amsterdam, 1983.

Kal2 [54] L. C. M. Kallenberg, “Survey of linear programming for standard and nonstandard Markovian control problems, Part I: Theory”, ZOR - Methods and Models in Operations Research, v.40, pp. 1-42, 1994.

KS [55] J.G. Kemeny and J.L. Snell, Finite Markov Chains, Von Nostrand, New York, 1960.

KT [56] V.S. Korolyuk and A.F. Turbin, Mathematical foundations of the state lumping of large systems, Naukova Dumka, Kiev, 1978, (in Russian), translated by Kluwer Academic Publishers, Dordrecht, The Netherlands, 1996.

Lam1 [57] B.F. Lamond, “A generalized inverse method for asymptotic linear programming”, Math. Programming, v.43, pp.71-86, 1989.

Lam2 [58] B.F. Lamond, “An efficient basis update for asymptotic linear programming”, Lin. Alg. Appl., v.184, pp.83-102, 1993.

LdPut [59] B. F. Lamond and M. L. Puterman, “Generalized inverses in discrete time Markov decision processes”, SIAM J. Matrix Anal. Appl., v.10, pp.118-134, 1989.

Ls [60] J.B. Lasserre, “A formula for singular perturbation of Markov chains”, J. Appl. Prob., v.31, pp.829-833, 1994.

LaLu [61] G. Latouche and G. Louchard, “Return times in nearly completely decomposable stochastic processes”, J. Appl. Prob., v.15, pp.251-267, 1978.

Lat [62] G. Latouche, “First passage times in nearly decomposable Markov chains”, in Numerical solution of Markov chains, Workshop 1990, Pure Prob. Appl., v.8, pp.401-411, 1991.

LuLa [63] G. Louchard and G. Latouche, “Geometric bounds on iterative approximations for nearly completely decomposable Markov chains”, J. Appl. Prob., v.27, pp.521-529, 1990.

Me [64] C.D. Meyer, “The role of the group generelized inverse in the theory of finite Markov chains,” SIAM Review, Vol. 17, pp. 443–464, 1975.

MV [65] B. L. Miller and A. F. Veinott, Jr., “Discrete dynamic programming with a small interest rate”, Ann. Math. Stat. v.40, pp.366-370, 1969.

MF [66] Ming Chen and J.A. Filar (1992),“Hamiltonian Cycles, Quadratic Programming and Ranking of Extreme Points” in Global Optimization, C. Floudas and P. Pardalos, eds. Princeton University Press.

PG [67] A.A. Pervozvanski and V.G. Gaitsgori, Theory of suboptimal decisions, Kluwer Academic Publishers, Dordrecht, The Netherlands, 1988, (Translation from the Russian orginal: Decomposition, aggregation and approximate optimization, Nauka, Moscow, 1979.)

PS [68] A.A. Pervozvanskii and I.N. Smirnov, “Stationary-state evaluation for a complex system with slowly varying couplings”, Cybernetics, v.10, pp.603-611, 1974. (Translation of Russian original in Kibernetika, v.10, pp.45-51, 1974.)

PhKo [69] R.G. Phillips and P.V. Kokotovic, “A singular perturbation approach to modeling and control of Markov chains”, IEEE Trans. Auto. Contr., AC-26, no.5, pp.1087-1094, 1981.

Put [70] M. L. Puterman, Markov decision processes, John Wiley & Sons, New York, 1994.

Q [71] J.P. Quadrat, “Optimal control of perturbed Markov chains: the multitime scale case”, in Singular perturbations in systems and control, ed. M.D. Ardema, CISM Courses and Lectures no.280, Springer-Verlag, Wien - New York, 1983.

RW [72] J.R. Rohlicek and A.S. Willsky, “The reduction of Markov generators: An algorithm exposing the role of transient states”, JACM, v.35, pp. 675-696, 1988.

Rt [73] U. Rothblum, Resolvent expansions of matrices and applications, Linear Algebra Appl., 38, (1981) 33–49.

Sc1 [74] P.J. Schweitzer, “Perturbation theory and finite Markov chains,” J. Appl. Probability, Vol. 5, pp. 401–413, 1968.

Sc2 [75] P.J. Schweitzer, “Perturbation series expansion of nearly completelydecomposable Markov chains,” Working Paper Series No. 8122, The Graduate School of Management, The University of Rochester, 1981.

Sc3 [76] P.J. Schweitzer, “Perturbation series expansion of nearly completelydecomposable Markov chains,” appeared in Teletraffic Analysis and Computer Performance Evaluation, by O.J. Boxma, J.W. Cohen and H.C. Tijms (editors), pp. 319–328, Elsevier Science Publishers B.V. (North Holland), 1986.

SA [77] H.A. Simon and A. Ando, “Aggregation of variables in dynamic systems”, Econometrica, v.29, pp.111-138, 1961.

Vat [78] H. Vantilborgh, “Aggregation with an error of O(ε^{2})”, Journal of the Association for Computing Machinery, v.32, pp.161-190, 1985.

Vt1 [79] A. F. Veinott, Jr., “Discrete dynamic programming with sensitive discount optimality criteria”, Ann. Math. Stat. v.40, pp.1635-1660, 1969.

Vt2 [80] A. F. Veinott, Jr., “Markov decision chains”, in Studies in Optimization, eds. G. B. Dantzig and B. C. Eaves, pp. 124-159, 1974.

YZ [81] G.G. Yin and Q. Zhang, “Continuous-time Markov chains and applications: A singular perturbation approach”, Series: Applications of Mathematics, v.37, Springer-Verlag, New York, 1998.