MARKOV DECISION PROCESSES Models, Methods, Directions, and Open Problems

Sean P. Meyn

Department of Electrical and Computer Engineering and the Coordinated Sciences Laboratory University of Illinois at Urbana-Champaign

Urbana, IL 61801, U.S.A.

May 28, 2000

Chapter 1

Stability, Performance Evaluation, and Optimization

Abstract

In this chapter we discuss various aspects of stability and control for MDP models. These ideas are centered around fluid models as approximations to the MDP model, and stochastic Lyapunov functions for verifying stability and bounding performance.

1.1 Introduction

This chapter considers stability of controlled Markov chains, and the relationship between criteria for stability; current approaches to bounding performance for the chain; and the construction of optimal policies. With stability taken as a central issue, the reader must soon realize that we will be considering large state spaces. In this chapter very little structure will be imposed on 𝕏, or its associated sigma field F. In other ways however our viewpoint will be fairly narrow. This chapter only considers the cost minimization problem, and only the average cost optimality criterion. Moreover, the assumptions that we impose will imply that the average cost is independent of the starting point of the process. By restricting attention in this way we hope that we can make the methodology more transparent.

One sees in several chapters in this volume that the generalization from finite state spaces to countable state spaces can lead to considerable technicalities. In particular, invariant distributions may not exist, and the cost functions of interest may not take on finite values. It would be reasonable to assume that the move from countable state spaces, to MDPs on a general state space should be at least as difficult. This assumption is probably valid if one desires a completely general theory.

However, the MDPs that we typically come across in practice exhibit structure which simplifies analysis, sometimes bringing us to the level of difficulty found in the countable, or even the finite state space case. For example, all of the specific models to be considered in this chapter, and most in this volume, have some degree of spatial homogeneity. The processes found in most applications will also exhibit some level of continuity in the sense that from similar starting points, and similar control sequences, two realizations of the model will have similar statistical descriptions. We do not require strong continuity conditions such as the strong Feller property. An assumption of ψ-irreducibility, to be described and developed below, allows one to lift most results in the discrete state space setting to models on a completely general, non-countable state space. This is an exceptionally mild assumption on the model and, without this assumption, the theory of MDPs on a general state space is currently extremely weak.

This surprising fact, that one does not lose much in moving from countable state space to ψ-irreducible chains on a general state space, can be regarded as a consequence of the Nummelin-Athreya-Ney construction of an ‘atom’. This construction involves a time homogeneous Markov chain, where the atom θ is a single point in an enlarged state space. The ψ irreducibility assumption makes this construction possible, and moreover ensures that the atom is reachable with positive probability from any initial condition. In this setting it is clear how to generalize many of the well known concepts in a finite state space setting. An invariant measure will be given by

μ{A} = 𝔼_{θ}[∑_{t=1}^{τ_{θ}} 1_{A}(x_{t})], A ∈ F,

where 1_{A} is the indicator function of the set A. In words, the quantity μ{A} expresses the mean number of times that the chain visits the set A before returning to θ. This expression assumes that the return time τ_{θ} is almost surely finite. If the mean return time 𝔼_{θ}[τ_{θ}] is finite then in fact the measure μ is finite, and it can then be normalized to give an invariant probability measure. Finiteness of the mean return time to some desirable state is the standard stability condition used for Markov chains, and for MDPs in which one is interested in the average cost optimality criterion.

Unfortunately, the split chain construction is cumbersome when developing a theory for controlled Markov chains. The sample path interpretation given above for the invariant probability μ is appealing, but it will be more convenient to work within an operator-theoretic framework. To motivate this, suppose first that we remain in the previous setting with an uncontrolled Markov chain, and suppose that do have an atom satisfying ℙ_{θ}{τ_{θ} < ∞} = 1. Denote by s the function which is equal to one at θ, and zero elsewhere: That is, s = 1_{θ}. We let ν denote the probability measure on 𝕏 given by ν(A) = P (θ, A), A ∈ F, and define the ‘outer product’ of s and ν by

s ⊗ ν (x, A) =^{Δ} s(x)ν(A).

We can then write in operator theoretic notation,

ℙ_{θ}{τ_{θ} ≥ n, x_{n} ∈ A} = ν(P − s ⊗ ν)^{n−1}1_{A}, n ≥ 1.

For example, in the finite state space case the measure ν can be interpreted as a row vector, the functions s and 1_{A} as column vectors, the kernel s ⊗ ν is the standard (outer) product of these two vectors, and any kernel such as P or s ⊗ ν can be interpreted as an N × N matrix, where N is the number of states.

Hence the invariant measure μ can be expressed in this notation as

μ(A) = ∑_{n=1}^{∞} ν(P − s ⊗ ν)^{n−1}1_{A},  A ∈ F.  (1.1) e:inv

It is this algebraic description of μ that will be generalized and exploited in this chapter.

How can we mimic this algebraic structure without constructing an atom? First, we require a function s: 𝕏 → ℝ_{+} and a measure ν on F satisfying the minorization condition,

P (x, A) ≥ s(x)ν(A), x ∈ 𝕏, A ∈ F.

In operator theoretic notation this is written P ≥ s⊗ν, and in the countable state space case this means that the transition matrix P dominates an outer product of two vectors with non-negative entries. For this bound to be useful we usually require s to be strictly positive on some suitably “large” set, and we require ν(𝕏) > 0. Unfortunately, this assumption excludes a large class of models, even the simple linear models to be considered as examples below. One can however move to the resolvent kernel defined by

K(x, A) = (1 − β) ∑_{t=0}^{∞} β^{t}P^{t} (x, A),  x ∈ 𝕏, A ∈ F,  (1.2) e:K

where β ∈]0, 1[ is some fixed constant. For a ψ-irreducible chain the required minorization always holds for the resolvent K. Much of the analysis then will involve the potential kernel, defined via

H(x, A) =^{Δ} ∑_{t=0}^{∞}(K − s ⊗ ν)^{t}(x, A),  x ∈ 𝕏, A ∈ F.  (1.3) e:potentialA

The move to the resolvent is useful since almost any object of interest can be mapped between the resolvent chain, and the original Markov chain. In particular, the invariant measures for P and K coincide.

The next question is stability. The ψ-irreducibility assumption will imply a reachability condition on s, analogous to the irreducibility condition that ℙ_{x}{τ_{θ} < ∞} > 0 for all x ∈ 𝕏. Generalizations of the more useful stability conditions such as L_{p} bounds on τ_{θ} are less straightforward. We will summarize known results in Section ^{\underline{s:lyapunov}} 1.2. A key observation is that the most natural stability assumption is equivalent to the existence of a Lyapunov function, whose form is very similar to the Poisson equation found in the average cost optimality equations (or ACOE). This connection will be exploited in our development of dynamic optimization theory below.

We conclude with an outline of the topics to follow. In the next section we review a bit of the general theory of ψ-irreducible chains, and develop some stochastic Lyapunov theory for such chains following ^{\underline{MT}} [39]. Following this, in Section ^{\underline{s:fish}} 1.3 we develop in some detail the computation of the average cost through the Poisson equation, and the construction of bounds on the average cost. All of these results are developed for time homogeneous chains without control. In Section ^{\underline{s:opt}} 1.4 the stability theory is applied to the analysis of the value iteration and policy iteration algorithms. Sections ^{\underline{s:linearEx}} 1.5 and ^{\underline{s:netEx}} 1.6 illustrate the theory with a detailed application to linear models, and to network scheduling. The chapter is concluded with a discussion of some open problems.

1.2 Stability

s:lyapunov

In this section we consider a Markov chain x with uncontrolled transition function P. The state space 𝕏 is assumed to be a locally compact, separable metric space, and we let F denote the (countably generated) Borel σ-field on 𝕏. Unless other references are given, all of the results described here together with their derivations can be found in ^{\underline{MT}} [39].

1.2.1 ψ-irreducibility

Throughout this chapter we assume that ψ is a σ-finite measure on F. The chain is called ψ-irreducible if the resolvent kernel defined in (^{\underline{e:K}} 1.2) satisfies

K(x, A) > 0, x ∈ 𝕏 ⟺ ψ(A) > 0.

We then call ψ a (maximal) irreducibility measure. We let F^{+} denote the set of all measurable h: 𝕏 → ℝ_{+} satisfying ψ(h) > 0. For A ∈ F we write A ∈ F^{+} provided 1_{A} ∈ F^{+}. That is, ψ(A) > 0. If the chain is ψ-irreducible, then from any initial condition x, the process has a chance of entering any set in F^{+} in the sense that ℙ_{x}{τ_{A} < ∞} > 0, where τ_{A} is the first return time,

τ_{A} = min{t ≥ 1 : x_{t} ∈ A}.

A function s: 𝕏 → ℝ_{+} and a non-trivial measure ν on F are called petite if the resolvent K satisfies the minorization condition,

K(x, A) ≥ s(x)ν(A),  x ∈ 𝕏, A ∈ F.  (1.4) e:accessible

One can show that a Markov chain is ψ-irreducible if and only if there is a function s: 𝕏 → (0, 1) and a probability measure ν satisfying (^{\underline{e:accessible}} 1.4). This bound is the most powerful consequence of the ψ-irreducibility assumption since it allows the construction of the potential kernel (^{\underline{e:potentialA}} 1.3).

A set C ⊆ 𝕏 is called petite if its indicator function 1_{C} is a petite function. Equivalently, for some probability measure ν and δ > 0,

K(x, A) ≥ δν(A), x ∈ C, A ∈ F.

The superlevel set {x : s(x) ≥ η} is always petite, for any η > 0, if the function s is petite. Consequently, for a ψ-irreducible chain, there always exists a countable covering of the state space by petite sets. The “petite property” can also be defined using hitting times: It is not difficult to show that, for a ψ-irreducible chain, the set C is petite if for each A ∈ F^{+}, there exists n ≥ 1, and δ > 0 such that

ℙ_{x}(τ_{A} ≤ n) ≥ δ  for any x ∈ C.  (1.5) e:PetiteTime

In applications we typically find that the function s used in (^{\underline{e:accessible}} 1.4) can be taken positive everywhere, and continuous. In this case we find that every compact set is petite, in which case the Markov chain is called a T -chain. It will be convenient to restrict attention to T -chains in Section ^{\underline{s:opt}} 1.4.

1.2.2 Recurrence and stability

The fundamental stability assumption for a Markov chain is the property that the state visit ‘important’ sets with probability one from any starting point. A ψ-irreducible chain is called recurrent if there exists a set 𝕏_{0} ∈ F satisfying ψ{𝕏_{0}^{c}} = 0, and ℙ_{x}{τ_{A} < ∞} = 1 for any A ∈ F^{+}, and any x ∈ 𝕏_{0}. If in addition the chain admits an invariant probability measure μ, then the chain is called positive recurrent. Although these definitions require us to consider an infinite number of initial states and sets A, there are several equivalent characterizations of recurrence which are easier to verify.

t:recurrence Theorem 1 The following are equivalent for a ψ-irreducible Markov chain x:

(i) x is recurrent.

(ii) For some petite set C,

ℙ_{x}{τ_{C} < ∞} = 1, x ∈ C.

(iii) For any pair (s, ν) satisfying the minorization condition (^{\underline{e:accessible}} 1.4) with s ∈ F^{+},

νHs =^{Δ} ∑_{t=0}^{∞} ν(K − s ⊗ ν)^{t}s = 1.

If any of these three equivalent conditions hold, then there exists a σ-finite measure μ which is invariant for the kernel P. It is unique in the sense that any σ-finite invariant measure is a constant multiple of the measure given by

μ^{∗}(A) = νH {A} = ∑_{t=0}^{∞} ν(K − s ⊗ ν)^{t} {A},  A ∈ F.  (1.6) e:inv2

Proof. The proof can be found in ^{\underline{num84,MT}} [44, 39]. However it will be useful to explain the main ideas, and in particular show why μ^{∗} defines an invariant measure.

To show why (iii) should hold for a recurrent chain we first make some definitions. To avoid dealing with potentially infinite sums, let λ > 0 and define the kernels

H_{λ}(x) = ∑_{0}^{∞} λ^{−n−1}(K − s ⊗ ν)^{n} G_{λ}(x) = ∑_{0}^{∞} λ^{−n−1}K^{n}.

Note that H_{1} = H is the potential kernel. These kernels are uniformly bounded in x for any λ > 1. We denote α_{λ} = ν(H_{λ}s) and β_{λ} = ν(G_{λ}s).

Applying the kernel λ^{−1}K to the function H_{λ}s gives,

λ^{−1}KH_{λ}s = λ^{−1}(K − s ⊗ ν)H_{λ}s + (s ⊗ ν)H_{λ}s
= H_{λ}s − λ^{−1}(1 − α_{λ})s

Iterating this equation we find that for λ > 1,

H_{λ}s − λ^{−1}(1 − α_{λ}) ∑_{0}^{n−1} λ^{−i}K^{i}s = λ^{−n}K^{n}H_{λ}s → 0, n → ∞.

This shows that H_{λ}s = (1 − α_{λ})G_{λ}s, and hence that α_{λ} = (1 − α_{λ})β_{λ} for λ > 1. Since all of the terms are positive we also see that α_{λ} < 1. Letting λ ↓ 1 and applying the monotone convergence theorem we do see that α_{1} ≤ 1 and in fact that,

β_{1} = \frac{α_{1}}{1 − α_{1}} ,

from which we see that νHs = α_{1} = 1 if and only if β_{1} = ∞. It remains to show that an infinite value for β_{1} is equivalent to recurrence. This is not difficult to believe, but due to lack of space we ask that the reader look elsewhere, e.g.^{\underline{num84,MT}} [44, 39].

To establish invariance of μ^{∗}, first apply the kernel (K − s ⊗ ν) to μ^{∗} on the right to obtain,

μ^{∗}(K − s ⊗ ν) = ∑_{t=0}^{∞} ν(K − s ⊗ ν)^{t+1} = μ^{∗} − ν.

Now by recurrence and (iii) we have μ^{∗}(s) = 1, which shows that μ^{∗} is K-invariant: μ^{∗}K = μ^{∗}. Using the identity

P K = KP = β^{−1}K + (1 − β^{−1})I,  (1.7) e:resolvEqn

one can conclude that μ^{∗} must then be P -invariant.

The invariant measure given in (^{\underline{e:inv2}} 1.6) will be finite, so the chain x is positive recurrent, provided that the mean return time to a petite set C is bounded:

sup_{x∈C} 𝔼_{x}[τ_{C}] < ∞.  (1.8) e:mean

In terms of the variables used in the previous proof, this is equivalent to requiring that α′_{1} < ∞, where the prime denotes the left derivative of α with respect to λ.

While these definitions lead to an elegant theory, in practice one can typically take 𝕏_{0} = 𝕏 in the definition of recurrence. In this case the chain is called Harris, and it is called positive Harris if there is also an invariant probability measure. The chains we consider next exhibit a far stronger form of stability.

1.2.3 c-Regularity and Lyapunov functions

The next level of stability that we consider is more closely related to steady state performance, which moves us closer to the average cost optimality criterion. Suppose that c: 𝕏 → [1, ∞) is a measurable function on the state space. For a ψ-irreducible chain, a set S ∈ F is called c-regular if for any A ∈ F^{+},

sup_{x∈S} 𝔼_{x}[ ∑_{t=0}^{τ_{A}−1} c(x_{t})] < ∞.

From the characterization in (^{\underline{e:PetiteTime}} 1.5) we see that a c-regular set is always petite. The Markov chain is called c-regular if the state space 𝕏 admits a countable covering by c-regular sets. A c-regular chain is automatically positive Harris, and by using the previous representation of the invariant measure one can show that a c-regular chain possesses an invariant probability measure μ satisfying μ(c) =^{Δ} ∫ c(x) μ(dx) < ∞. The following result is a consequence of the f-Norm Ergodic Theorem of ^{\underline{MT}} [39, Theorem 14.0.1].

t:little-f-norm Theorem 2 Assume that c: 𝕏 → [1, ∞) and that x is c-regular. Then, for any measurable function g which satisfies

sup_{x∈𝕏}(\frac{|g(x)|}{c(x)}) < ∞,

the following ergodic theorems hold for any initial condition:

(i)  lim_{n→∞} \frac{1}{n} ∑_{t=1}^{n} g(x_{t}) = μ(g),  a.s. [ℙ_{x}].
(ii) lim_{n→∞} \frac{1}{n} ∑_{t=1}^{n} 𝔼_{x}[g(x_{t})] = μ(g).

Thus we see that one level of complexity is removed when considering c-regular chains since the steady state performance as captured by μ(g) does not depend upon the initial condition of the chain.

One approach to establishing c-regularity is through the following extension of Foster’s criterion, or Lyapunov’s second method. In general, such approaches involve the construction of a function V on the state space, taking positive values, such that V (x_{t}) is in some sense decreasing whenever the state x_{t} is ‘large’. In our context this decreasing property can be formulated as follows: Find a function V : 𝕏 → ℝ_{+} and a constant J̄ ∈ ℝ_{+} such that

P V (x) =^{Δ} 𝔼[V (x_{t+1}) | x_{t} = x] ≤ V (x) − c(x) + J̄,  x ∈ 𝕏.  (1.9) e:cFoster

However, for this to imply any form of stability, the difference c(x)−J̄ must be positive for ‘large’ x. This requires some assumptions on c. We say that c is near-monotone if the sublevel set S_{η} =^{Δ} {x ∈ 𝕏 : c(x) ≤ η} is petite for any η < ∥c∥_{∞}. The supremum norm ∥c∥_{∞} = sup_{x∈𝕏} c(x) may be infinite. The function is called norm-like if the sublevel set S_{η} is a precompact subset of the metric space 𝕏 for any η. Related assumptions on c are used in ^{\underline{araborferghomar93,bor91,MT}} [1, 5, 39].

t:LyapunovRegular Theorem 3 Assume that c: 𝕏 → [1, ∞) is near-monotone, and suppose that J̄ < ∥c∥_{∞}. Then, (i) If there exists a finite, positive-valued solution V to the inequality (^{\underline{e:cFoster}} 1.9),

then there exists d_{0} < ∞ such that for each A ∈ F^{+},

𝔼_{x}[∑_{t=0}^{τ_{A}} c(x_{t})] ≤ d_{0}V (x) + d(A),  x ∈ 𝕏,  (1.10) e:c-bdd

where d(A) < ∞ is a constant. Hence, each of the sublevel sets S_{n} = {x : V (x) ≤ n} is c-regular, and the process itself is c-regular.

(ii) If the chain is c-regular, then for any c-regular set S ∈ F^{+}, the function

V^{∗}(x) = 𝔼_{x}[∑_{t=0}^{τ_{S}} c(x_{t})],  x ∈ 𝕏,  (1.11) e:V-star

is a near-monotone solution to (^{\underline{e:cFoster}} 1.9).

Proof. The bound (^{\underline{e:cFoster}} 1.9) is equivalent to the drift condition P V_{0} ≤ V_{0} − c+b1_{S}, where S is petite: if (^{\underline{e:cFoster}} 1.9) holds, we can take V_{0} = d_{0}V and b = d_{0}J̄, with d_{0} sufficiently large. The result is then an immediate consequence of \underline{MT} [39, Theorem 14.2.3].

1.3 Performance

s:fish

1.3.1 Poisson’s equation

Poisson’s equation originated in the analysis of partial differential equations: Assuming that f is some given function on ℝ^{n}, the equation is written

Δh = −f

where h is an unknown function on ℝ^{n}, and Δ is the Laplacian. The probabilistic interpretation of this equation becomes evident when one realizes that Δ is the generator for a Brownian motion on ℝ^{n} - a similar equation can be posed for any Markov process in continuous time. When time is discrete, we then define the generator as Δ = P − I, and the Poisson equation then takes on the exact same form, where we take f = c − μ(c).

The probabilistic motivation for looking at this equation follows from our prior stability analysis. First note that the drift inequality (^{\underline{e:cFoster}} 1.9) suggests a simple approach to obtaining performance bounds. By iterating this equation one obtains,

0 ≤ 𝔼_{x}[V (x_{n})] ≤ V (x) − ∑_{t=0}^{n−1} 𝔼_{x}[c(x_{t})] + nJ̄.

Dividing by n and letting n → ∞ then gives the upper bound,

J(c) =^{Δ} lim sup_{n→∞} \frac{1}{n} 𝔼_{x}[∑_{t=0}^{n−1} c(x_{t})] ≤ J̄.

The question then is, by choosing V carefully can we get a tight upper bound? The answer is yes, provided the chain is c-regular, and in this case, the minimal upper bound J̄_{∗} is evidentally μ(c), with μ equal to the invariant probability for the chain. The following result surveys the relevant consequences of c-regulariry, and introduces the form of the Poisson equation that we will analyze in the remainder of this chapter.

t:fish Theorem 4 Assume that c: 𝕏 → [1, ∞) and that x is c-regular. Then,

(i) There exists a measurable function h: 𝕏 → ℝ satisfying the Poisson equation

P h (x) =^{Δ} 𝔼_{x}[h(x_{t+1}) | x_{t} = x] = h(x) − c(x) + J,  x ∈ 𝕏, (1.12) e:fish

where J = μ(c).

(ii) One solution to (^{\underline{e:fish}} 1.12) may be expressed,

h(x) = KHc̄(x) = ∑_{t=1}^{∞}(K − s ⊗ ν)^{t}c̄(x),  x ∈ 𝕏,  (1.13) e:fishFormula

where c̄ = c − J, and the pair (s, ν) is petite with s ∈ F^{+}.

(iii) Suppose moreover that the function c is near-monotone. Then the solution h to (^{\underline{e:fishFormula}} 1.13) is uniformly bounded from below, inf_{x∈𝕏} h(x) > −∞. It is essentially unique in the following sense: If h′ is any function on 𝕏 which is uniformly bounded from below, and solves the Poisson inequality

P h (x) ≤ h(x) − c(x) + J, x ∈ 𝕏,

with J = π(c), then there exists a constant a such that

h′(x) ≥ h(x) + a,  x ∈ 𝕏;
h′(x) = h(x) + a,  a.e. x ∈ 𝕏 [ψ].

(iv) If V is any solution to (^{\underline{e:cFoster}} 1.9) with J̄ < ∥c∥_{∞} and c near-monotone, then the solution (^{\underline{e:fishFormula}} 1.13) satisfies the uniform upper bound, for some d_{0} < ∞,

h(x) ≤ d_{0}(V (x) + 1), x ∈ 𝕏.

Proof. Again, the proof can be found elsewhere (see ^{\underline{mey97b}} [42], following ^{\underline{num91,glymey96a}} [45, 22]). However, the fact that the solution to the Poisson equation (^{\underline{e:fish}} 1.12) can be taken as in (^{\underline{e:fishFormula}} 1.13) is easy to explain, and shows some attractive symmetry with the earlier construction of an invariant measure. Consider first the function

h_{0}(x) = Hc̄ =^{Δ} ∑_{t=0}^{∞}(K − s ⊗ ν)^{t}c̄(x), x ∈ 𝕏.

Observe that by the definition of J and c̄, and the construction of μ^{∗}, we must have ν(h_{0}) = 0. Thus,

Kh_{0} (x) = (K − s ⊗ ν)h_{0} (x) = h_{0}(x) − c̄(x).

That is, h_{0} solves the Poisson equation for the kernel K. By again applying the identity (^{\underline{e:resolvEqn}} 1.7) we see that h = Kh_{0} = (K − s ⊗ ν)h_{0} solves the Poisson equation for original transition kernel P.

1.3.2 Simulation

We have now seen that the Poisson equation has a direct role in performance evaluation. Although we have not given any explicit algorithms, it is clear that an approximation of the solution h will lead to an approximation of J = μ(c). With some structure imposed on the model this idea does lead to algorithms for computing bounds on J. For example, this is the essence of the main results in ^{\underline{kummey96a,kumkum94a}} [35, 36], where performance bounds are obtained in the network scheduling problem. If the cost is linear, and if any of the linear programs constructed in these references admits a feasible solution, then the solution to Poisson’s equation is approximated by a pure quadratic function.

Perhaps the most obvious approach to estimating J is through Monte Carlo simulation via

Ĵ_{n} = \frac{1}{n} ∑_{0}^{n−1} c(x_{t}), n ∈ ℕ.

The Poisson equation again plays an important role in analysis, and in the generation of more efficient simulation approaches.

The effectiveness of the Monte Carlo method depends primarily on the magnitude of the Central Limit Theorem variance, also known as the timeaverage variance. Under suitably strong recurrence conditions on the Markov chain this can be expressed

γ_{c}^{2} = lim_{n→∞} 𝔼_{x}[(\frac{1}{\sqrt{n}} ∑_{0}^{n−1} c̄(x_{t}))^{2}]

An alternative expression for the time-average variance is computed through the formula

γ_{c}^{2} = μ(h^{2}) − μ((P h)^{2}) = 2μ(hc̄) − μ(c̄^{2}),  (1.14) e:cltVar

with h any solution to (^{\underline{e:fish}} 1.12) ^{\underline{MT}} [39].

There are many variants of the simple Monte Carlo estimate, some of which may have far smaller variance. After all, if {Δ_{t} : t ≥ 0} is any sequence of random variables satisfying \frac{1}{n} ∑_{0}^{n−1} Δ_{t} → 0, n → ∞, then the modified estimator,

Ĵ_{n}^{Δ} = \frac{1}{n} ∑_{0}^{n−1}(c(x_{t}) + Δ_{t}), n ∈ ℕ,

is another consistent estimator of J. An optimal choice for Δ_{t} is computed using the solution h to Poisson’s equation (^{\underline{e:fish}} 1.12): by setting

Δ_{t}^{∗} = P h (x_{t}) − h(x_{t}),

we obtain a time-average variance of zero. Of course, computing Δ_{t}^{∗} involves a computation of J, so this approach is nonsensical! If however an approximation g to h can be found, then the choice Δ_{t} = P g (x_{t}) − g(x_{t}) will lead to reduced variance if the approximation is sufficiently tight ^{\underline{hen97a}} [24].

This is a useful result for our purposes since we will discover such approximations when we attempt to solve some optimization problems below.

1.3.3 Examples

In this chapter we will restrict ourselves to two general examples: the linear state space model, and a family of network models. In this section we look at some special cases without control. Controlled linear systems, and controlled network models are considered as examples in the final two sections in this chapter.

The linear state space model is defined through the multidimensional recursion,

x_{t+1} = Ax_{t} + w_{t+1},  t ∈ ℕ,  (1.15) e:linear-0

where x_{t}, w_{t} ∈ ℝ^{d}, and w is i.i.d. with w_{t} ∼ N(0, Σ). Let F by any matrix of suitable dimension satisfying F F^{T} = Σ. If F is d × q for some q, then the controllability matrix is the d × (dq) matrix C =^{Δ} [A^{d−1}F |A^{d−2}F | ⋯ |AF |F], and the pair (A, F) is called controllable if the matrix C has rank d. The process is ψ-irreducible with ψ equal to Lebesgue measure if the pair (A, F) is controllable, since in this case P^{t}(x, ·) is equivalent to Lebesgue measure for any x, and any t ≥ d. By continuity of the model it is easy to check that (^{\underline{e:accessible}} 1.4) holds with s continuous, and ν equal to normalized Lebesgue measure on an open ball in ℝ^{d}. We conclude that all compact sets are petite if the controllability condition holds, so that x is a T -chain. A sample path from a particular two dimensional linear model is shown in Figure ^{\underline{f:lqg}} 1.1 where it can be seen that when the state is large, the sample path behavior appears almost deterministic.

Figure 1.1: A sample path of the linear state space model with a ‘large’ initial condition x_{0} = \binom{−100}{100}.

f:lqg

To find a stochastic Lyapunov function V with c(x) = \frac{1}{2}x^{T} Qx, first solve the Lyapunov equation

A^{T} P A = P − Q.  (1.16) e:lyapEqn

If P ≥ 0, then V (x) = \frac{1}{2}x^{T} P x is a solution to (^{\underline{e:cFoster}} 1.9). In fact, the function V is also the essentially unique solution to Poisson’s equation, with J = \frac{1}{2}trace (P Σ).

For the nonlinear state space model

x_{t+1} = F (x_{t}, w_{t+1}), t ∈ ℕ,

the ψ-irreducibility condition can still be verified under a nonlinear controllability condition called forward accessibility ^{\underline{MT}} [39, Chapter 7]. The construction of a Lyapunov function is however far more problem-specific.

Over the past five years there has been much research on algorithmic methods for constructing Lyapunov functions for network models. One is based upon linear programming methods, and is similar to the Lyapunov equation (^{\underline{e:lyapEqn}} 1.16) used for linear state space models ^{\underline{kummey96a}} [35]. We describe here a recent approach based upon a fluid model ^{\underline{mey97a,mey97b}} [41, 42]. As an example we consider here the simplest case: An uncontrolled M/M/1 queue.

When the arrival stream is renewal, and the service times are i.i.d., then the waiting time for a simple queue can be modeled as a Markov chain with state space 𝕏 = ℝ_{+}. The dynamics take the form of a one dimensional linear state space model, where the state space is constrained to the positive half line. The queue length process is itself a Markov process in the special case where the service times and interarrival times are exponentially distributed. By applying uniformization (i.e. sampling the process appropriately - see \underline{lip75} [38]), the queue length process x obeys the recursion

x_{t+1} = x_{t} + (1 − I_{t+1})π(x_{t}) + I_{t+1}, t ∈ ℕ,

where I is a Bernouli, i.i.d. random process: λ = ℙ(I(t) = 1) is the arrival rate, and μ = ℙ(I(t) = −1) is the service rate. The function π plays the role of a ‘policy’, where in this simple example we take π(x) = 1(x > 0). Time has been normalized so that λ + μ = 1.

To construct a Lyapunov function, first note that stability is a ‘large state’ property, so it may pay to consider the process starting from a large initial condition. In the left hand side of Figure ^{\underline{f:mm1}} 1.2 we see one such simulation. As was seen in the linear model, when the initial condition is large the behavior of the model is roughly deterministic.

Suppose we take the cost function c(x) = x. To construct a Lyapunov function we would ideally like to compute the expected sum given in (^{\underline{e:V-star}} 1.11), with S equal to some finite set, perhaps S = {0}. While this is computable for the M/M/1 queue, such computation can be formidable for more complex network models. However, consider the right hand side of Figure ^{\underline{f:mm1}} 1.2 which shows a sample path of the deterministic fluid, or leaky bucket model. This satisfies the differential equation ϕ̇ = (−μ + λ)π(ϕ), where π is again equal Figure 1.2: On the left is a sample path x_{t} of the M/M/1 queue with ρ = λ/μ = 0.9, and x_{0} = 400. On the right is a solution to the differential equation ϕ̇ = (−μ + λ)π(ϕ) starting from the same initial condition.

f:mm1

the indicator function of the strictly positive real axis. The behavior of the two processes look similar when viewed on this large spatial/temporal scale. It appears that a good approximation is V^{∗}(x) ≈

\begin{aligned}V (x) & =^{Δ} & ∫_{0}^{∞} ϕ(t) dt, ϕ(0) = x, \\ & = & \frac{1}{2} \frac{x^{2}}{μ − λ}. \end{aligned}  (1.17) e:fluidValue

If we apply the transition kernel P to V we find, for x ≥ 1,

P V (x) = λV (x + 1) + μV (x − 1)
= \frac{1}{2(μ − λ)}(λ(x + 1)^{2} + μ(x − 1)^{2})
= V (x) − x + \frac{1}{2(μ − λ)},

while for x = 0 we have,

P V (x) = \frac{λ}{2(μ − λ)} ≤ V (x) − x + \frac{1}{2(μ − λ)}

That is, we see that this approach works: The stochastic Lyapunov criterion (^{\underline{e:cFoster}} 1.9) does hold with this function V derived from the fluid model, where J̄ = (2(μ − λ))^{−1}, under the stability condition that ρ = λ/μ < 1. The actual steady state mean of c(x) = x is given by J = λ(μ − λ)^{−1}, which is indeed upper-bounded by J̄.

What about the more exact Poisson’s equation? Can the fluid model be used to approximate a solution?

With the cost function c(x) = x, the Poisson equation for the M/M/1 queue becomes

P h = λh(x + 1) + μh((x − 1)^{+}) = h(x) − x + J.

One solution is given by

h(x) = \frac{x^{2} + x}{2(μ − λ)},

which is similar in form to the fluid value function given in (^{\underline{e:fluidValue}} 1.17).

For a general class of network models it can be shown that the value function for the fluid model and the solution to Poisson’s equation are roughly equal for large x in the sense that

h(x) = V (x)(1 + o(1))

where the term o(1) → 0 as x → ∞. Some results of this type are described in Section ^{\underline{s:netEx}} 1.6.

The M/M/1 queue illustrates nicely the difficulties one faces in using simulation to esimate steady state performance measures since the timeaverage variance constant γ_{c}^{2} grows quickly with the ‘system load’. To obtain bounds on γ_{c}^{2} note firstly that the solution to Poisson’s equation for the M/M/1 queue is a quadratic of the form h(x) = b(1 − ρ)^{−1}(x^{2} + x), where b is a bounded function of λ. Using (^{\underline{e:cltVar}} 1.14) we conclude that γ_{c}^{2} is of order (1 − ρ)^{−4} in this example since pth moments for the M/M/1 queue are of order (1 − ρ)^{−p}.

It is shown in ^{\underline{henmey99a}} [25] that this growth rate on γ_{c}^{2} continues to hold for a general class of network models. The large variance indicates that in heavy traffic (where ρ ∼ 1) it is computationally expensive to compute the mean performance through standard Monte Carlo simulation. This provides ample motivation to find methods for approximating solution’s to Poisson’s equation to construct reduced variance estimators for network models. The use of the fluid value function is one promising approach ^{\underline{henmey99a}} [25].

1.4 Optimization

s:opt

With this background we are now ready to turn to MDP models.

We now assume that there is a control sequence taking values in the action space 𝔸 which influences the behavior of x. The state space 𝕏 and the action space 𝔸 are assumed to be locally compact, separable metric spaces, and we continue to let F denote the Borel σ-field on 𝕏. Associated with each x ∈ 𝕏 is a non-empty and closed subset 𝔸(x) ⊆ 𝔸 whose elements are the admissible actions when the state process x_{t} takes the value x. The set of admissible state-action pairs {(x, a) : x ∈ 𝕏, a ∈ 𝔸(x)} is assumed to be a measurable subset of the product space 𝕏 × 𝔸.

The transitions of x are governed by the conditional probability distributions {P_{a}(x, B)} which describe the probability that the next state is in B for any B ∈ F given that the current state is x ∈ 𝕏, and the current action chosen is a ∈ 𝔸. These are assumed to be probability measures on F for each state-action pair (x, a), and measurable functions of (x, a) for each B ∈ F. We will primarily restrict attention to stationary Markov policies. Recall that such a policy, denoted π ∈ Π^{S}, is a measurable function π : 𝕏 → 𝔸 such that π(x) ∈ 𝔸(x) for all x. When the policy π is applied to the MDP, then the action π(x) is applied whenever the MDP is in state x, independent of the past and independent of the time-period. We shall write P_{π}(x, B) = P_{π(x)}(x, B) for the transition law corresponding to a policy π.

The state process x^{π} =^{Δ} {x_{t}^{π} : t ≥ 0} of the MDP is, for each fixed policy π, a Markov chain on (𝕏, F), and we write the t-step transition probabilities for this chain as

P_{π}^{t}(x, B) = ℙ(x_{t}^{π} ∈ B | x_{0}^{π} = x), x ∈ 𝕏, B ∈ F, t ∈ ℕ.

In the controlled case we continue to use the operator-theoretic notation,

P_{π}^{t}h (x) =^{Δ} 𝔼[h(x_{t}^{π}) | x_{0}^{π} = x]

In the remander of this section we will first recall the ACOE for an MDP, and then develop two algorithms which can be used to construct solutions to these equations. The value iteration algorithm (VIA) will be considered first, which is simply the successive approximation procedure. The results describe mainly stability of the policies generated by the algorithm, and some bounds on steady state performance. These results are based on ^{\underline{chemey99a}} [10]. The policy iteration algorithm (PIA), first proposed in ^{\underline{how60}} [31], may be viewed as a version of the Newton Raphson method, and is consequently considerably more complex than value iteration. It is however far easier to analyze.

We find under the assumptions we impose that it exhibits nearly monotone convergence, and that the algorithm does generate successively better policies. These results are taken from ^{\underline{mey97b}} [42].

More background and convergence results for such algorithms may be

found in ^{\underline{araborferghomar93,bor91,cav96,cavfer95c,hermoncav91,horput87,hor77,put94,sen96a}} [1, 5, 7, 8, 27, 29, 28, 48, 54].

1.4.1 Regular policies and the ACOE

We now suppose that a one-step cost function c: 𝕏 × 𝔸 → [1, ∞) is given: we assume below that c satisfies a near-monotone condition so that the results of Section ^{\underline{s:lyapunov}} 1.2 and ^{\underline{s:fish}} 1.3 may be applied. For any policy π we denote c_{π}(y) = c(y, π(y)), and we denote the steady state average cost by

J(π, x) =^{Δ} lim sup_{n→∞} \frac{1}{n} 𝔼_{x}[∑_{t=0}^{n−1} c_{π}(x_{t}^{π})].

A policy π_{∗} will be called optimal if J(π_{∗}, x) ≤ J(π, x) for all policies π, and any initial state x.

The policy π is called regular if the controlled chain x^{π} is a c_{π}-regular Markov chain. This is a natural and highly desirable stability property for the controlled process. If the policy is regular, then necessarily an invariant probability measure μ_{π} exists such that μ_{π}(c_{π}) < ∞. Moreover, for a regular policy, the resulting average cost is independent of the starting point of the chain: J(π, x) ≡ μ_{π}(c_{π}).

When J_{∗} = J(π_{∗}, x) is independent of x, then the associated ACOE are given as follows, where the function h_{∗} is known as the relative value function.

J_{∗} + h_{∗}(x) = min_{a∈𝔸(x)}[c(x, a) + P_{a}h_{∗} (x)]  (1.18) e:OE1
π_{∗}(x) = arg min_{a∈𝔸(x)}[c(x, a) + P_{a}h_{∗} (x)],  x ∈ 𝕏. (1.19) e:OE2

If a policy π_{∗}, a measurable function h_{∗}, and a constant J_{∗} exist which solve these equations, then typically the policy π_{∗} is optimal (see for example \underline{araborferghomar93,bor91,lerlas96,put94,sen99a} [1, 5, 26, 48, 53] for a proof of this and related results).

t:OE Theorem 5 Suppose that the following conditions hold

(a) The pair (J_{∗}, h_{∗}) solve the optimality equation (^{\underline{e:OE1}} 1.18); (b) The policy π_{∗} satisfies (^{\underline{e:OE2}} 1.19), so that

c_{π_{∗}}(x) + P_{π_{∗}}h_{∗} (x) ≤ c(x, a) + P_{a}h_{∗} (x), x ∈ 𝕏, a ∈ 𝔸(x).

(c) For any x ∈ 𝕏, and any policy π satisfying J(π, x) < ∞,

\frac{1}{n}P_{π}^{n}h_{∗} (x) → 0, n → ∞.

Then π_{∗} is an optimal control, and J_{∗} is the optimal cost, in the sense that

lim_{n→∞} \frac{1}{n} ∑_{t=1}^{n} 𝔼_{x}[c_{π_{∗}}(x_{t}^{π_{∗}})] = J_{∗},

and J(π, x) ≥ J_{∗} for all policies π, and all initial states x.

The assumption (c) is unfortunate, but examples show that some additional conditions on h_{∗} are required (see e.g. page 87 of ^{\underline{bor91}} [5], Chapter 7 of \underline{sen99a} [53], or the examples in ^{\underline{araborferghomar93,dek87,put94,ros92}} [1, 16, 48, 50]). The following result gives a condition implying (c) which is often verifiable in practice, as we shall see in Section ^{\underline{s:linearEx}} 1.5 and Section ^{\underline{s:netEx}} 1.6.

If the controlled chain x^{π} is ψ_{π}-irreducible and the resulting cost J_{π} =^{Δ} μ_{π}(c_{π}) = ∫ c_{π}(x)μ_{π}(dx) is finite, let S_{π} denote any fixed c_{π}-regular set for which μ_{π}(S_{π}) > 0. We then define the function

V_{π}(x) = 𝔼_{x}[ ∑_{t=0}^{τ_{π}−1} c_{π}(x_{t}^{π})],  (1.20) e:optV

where τ_{π} = τ_{S_{π}}. Since μ_{π}(S_{π}) > 0, the function V_{π} is a.e. [μ_{π}] finite-valued \underline{MT} [39, Theorem 14.2.5]. Note that by ^{\underline{MT}} [39, Theorem 14.2.3], the particular c_{π}- regular set S_{π} chosen is not important. If S_{π}^{1} and S_{π}^{2} give rise to functions V_{π}^{1} and V_{π}^{2} of the form (^{\underline{e:optV}} 1.20), then for some constant γ ≥ 1,

γ^{−1}V_{π}^{1}(x) ≤ V_{π}^{2}(x) ≤ γV_{π}^{1}(x), x ∈ 𝕏.

t:comparable Theorem 6 (Meyn ^{\underline{mey97b}} [42]) Suppose that

(a) The optimality equations (^{\underline{e:OE1}} 1.18,^{\underline{e:OE2}} 1.19) hold for (π_{∗}, h_{∗}, J_{∗}), with h_{∗} bounded from below;

(b) For any policy π the average cost K_{π}c_{π} is finite and norm-like, and all compact sets are petite for the Markov chain x^{π};

(c) For any policy π there exists some constant d_{0} = d_{0}(w) < ∞ such that,

|h_{∗}(x)| ≤ d_{0}V_{π}(x),  x ∈ 𝕏.  (1.21) e:relativeBound

Then π_{∗} is a regular, optimal policy.

Proof. To show that π_{∗} is regular we consider the Poisson equation for the Markov chain with transition kernel K_{π_{∗}}:

K_{π_{∗}}h_{∗} = h_{∗} − K_{π_{∗}}c_{π_{∗}} + J_{∗}.

Under the assumption that h_{∗} is bounded from below, it follows from Theorem ^{\underline{t:little-f-norm}} 2 that the “K_{π_{∗}}-chain” is K_{π_{∗}}c_{π_{∗}}-regular, and regularity of π_{∗} follows.

To prove optimality suppose that π is any policy, and note first that if J(π, x) = ∞ for all x, then there is nothing to prove. If not, then since all compact sets are petite, the Markov chain x^{π} is a positive recurrent T -chain, with unique invariant probability μ_{π}, and J_{π} =^{Δ} μ_{π}(c_{π}) is finite ^{\underline{MT}} [39, Theorem 14.0.1].

Under the assumptions of the theorem, we can show inductively that

P_{π}^{n}V_{π} = V_{π} − ∑_{t=0}^{n−1} P_{π}^{t}(c_{π} − s_{π}),

where s_{π} ≥ 0 and satisfies μ_{π}(s_{π}) = μ_{π}(c_{π}). This function can be written explicitly as

s_{π}(x) = ∫_{S_{π}} P_{π}(x, dy) 𝔼_{y}[ ∑_{t=0}^{τ_{π}−1} c_{π}(x_{t})]

It follows from ^{\underline{MT}} [39, Theorem 14.0.1] that P_{π}^{n}V_{π}(x)/n → 0 as n → ∞ for a.e. [μ_{π}] x ∈ 𝕏. By (^{\underline{e:relativeBound}} 1.21), we also have P_{π}^{n}h_{∗} (x)/n → 0, and it easily follows that J(π, x) ≥ J(π_{∗}) for such x. To generalize to arbitrary x we use the norm-like assumption on c which implies that J(π, x) ≥ μ_{π}(c_{π}) for all x ∈ 𝕏 (for details see ^{\underline{mey97b}} [42]).

1.4.2 Value iteration

s:via

The ACOE can be viewed as a fixed point equation in the variables (h_{∗}, J_{∗}). By ignoring the constant term and applying successive approximation to this fixed point equation we obtain the VIA: Suppose that the positive-valued function V_{n} is given. Then the policy π_{n} is defined as

π_{n}(x) = arg min_{a∈𝔸(x)}[P_{a}V_{n} (x) + c(x, a)], x ∈ 𝕏,

and one then defines

V_{n+1}(x) = c_{π_{n}}(x) + P_{π_{n}}V_{n}(x) = min_{a∈𝔸(x)}(P_{a}V_{n} (x) + c_{a}(x)),

which then makes it possible to compute the next policy π_{n+1}.

This is in fact the standard dynamic programming approach to constructing a finite horizon optimal policy since for each n we may write,

V_{n}(x) = min 𝔼_{x}[∑_{t=0}^{n−1} c(Φ(t), a(t)) + V_{0}(Φ(n))],  (1.22) e:ValueFunction

where {a(t) : t ≥ 0} is a sequence of actions determined by some Markov policy, and the minimum in (^{\underline{e:ValueFunction}} 1.22) is with respect to all such policies.

To simplify notation we define c_{n} = c_{π_{n}}, P_{n} = P_{π_{n}}, and we define the resolvent for the nth policy by

K_{n} =^{Δ} (1 − β) ∑_{t=0}^{∞} β^{t}P_{n}^{t},  n ≥ 0,  (1.23) e:Kn-def

where β ∈]0, 1[ as before. We let 𝔼^{n} denote the expectation operator induced by the stationary policy π_{n}.

Let ν denote some fixed finite measure on F - we will impose conditions below which ensure that ν(V_{n}) is finite-valued. We then define for each n the normalized value function, and the incremental cost,

h_{n}(x) = V_{n}(x) − ν(V_{n});  γ_{n}(x) = V_{n+1}(x) − V_{n}(x),  x ∈ 𝕏, n ∈ ℕ._{(1.24) e:hg-def}

From the definitions, for each n we have the familiar looking identity P_{n}h_{n} = h_{n} − c_{n} + γ_{n}. For this to imply any form of stability for the policy π_{n} we require that γ_{n} be bounded from above. When this is the case we define J̄_{n} = sup_{x} γ_{n}(x), so that we have the following solution to (^{\underline{e:cFoster}} 1.9):

P_{n}V_{n} ≤ V_{n} − c_{n} + J̄_{n}.  (1.25) e:almostFish-a

Useful bounds on γ_{n} will in general require conditions on the initial condition V_{0}. We show here that a sufficient condition is that V_{0} act as a Lyapunov function for some policy. We assume that at least one regular policy π_{−1} exists, and that for some constant J̄ < ∞,

P_{π_{−1}}V_{0} ≤ V_{0} − c_{π_{−1}} + J̄.  (1.26) e:stable-ic

The existence of a pair (V_{0}, π_{−1}) satisfying (^{\underline{e:stable-ic}} 1.26) is a natural stabilizability assumption on the model.

The assumptions below ensure that the algorithm can be initialized to generate stabilizing policies. Condition (A2) relates the average optimality problem with the discounted optimal control problem: this assumption is satisfied if the state dependent cost V_{β} is norm-like, where V_{β}(x) denotes the optimal cost for the discounted optimal control problem when the discount factor is equal to β, and the initial condition is x. Assumption (A3) is a uniform minorization condition on the MDP model.

(A1) For each n, if the VIA yields a value function V_{n}: 𝕏 → ℝ_{+}, then the minimization

π_{n}(x) =^{Δ} arg min_{a∈𝔸(x)}[c(x, a) + P_{a}V_{n} (x)]

admits a measurable solution π_{n}.

(A2) For each fixed x, the function c(x, ·) is norm-like on 𝔸, and there exists a norm-like function \underline{c}: 𝕏 → ℝ_{+} such that for the policies π_{n} obtained through the VIA,

∞ > K_{n}c_{n} (x) ≥ \underline{c}(x), for any x ∈ 𝕏, n ∈ ℕ.

(A3) There is a fixed probability ν on F, a δ > 0, and an initial value function V_{0} with the following property: For each n ≥ 1, if the VIA yields the value function V_{n}, then for any policy π_{n} given in (A1),

K_{n}(x, A) ≥ δν(A)  for all x ∈ S, A ∈ F,  (1.27) e:nuAccessibilityVIA

where S denotes the precompact set

S = \{x : \underline{c}(x) ≤ 2J̄\}.  (1.28) e:S-def-VIA

The following result is largely taken from ^{\underline{chemey99a}} [10].

t:VIA-c-regular Theorem 7 Suppose that (A1) and (A2) hold, and suppose that for each of the policies {π_{n}} obtained through the VIA, all compact sets are petite for the Markov chain x^{π_{n}}. Assume moreover that the initialization V_{0} satisfies (^{\underline{e:stable-ic}} 1.26). Then,

(i) Each of the policies {π_{i} : i ∈ ℕ} is regular.

(ii) The upper bounds {J̄_{n}} are decreasing:

J̄_{0} ≥ J̄_{1} ≥ ⋯ ≥ J̄_{n} ≥ ⋯;

(iii) If in addition Assumption (A3) holds, then the sequence {h_{n}} is uniformly bounded from below.

Proof. The main idea is to apply the following bound on the resolvent, which follows from (^{\underline{e:almostFish-a}} 1.25):

K_{n}V_{n} ≤ V_{n} − \frac{β}{1 − β} K_{n}c_{n} + \frac{β}{1 − β} J̄_{n}.  (1.29) e:KalmostFish

This inequality is a version of (^{\underline{e:cFoster}} 1.9) since V_{n} ≥ 0, provided that J̄_{n} is finite.

The minimization in the value iteration algorithm immediately leads to the bound P_{n}γ_{n} ≥ γ_{n+1}. From this we deduce by induction that the J̄_{n} are finite and decreasing, which proves (ii). The initialization of the induction relies on the assumption that the initial condition V_{0} satisfies (^{\underline{e:stable-ic}} 1.26). Applying Theorem ^{\underline{t:LyapunovRegular}} 3, which is applicable under (A2) for the kernel K_{n}, we see that the Markov chain with transition kernel K_{n} is c-regular. This implies (i).

To prove (iii) note first of all that h_{n}(x) ≥ −ν(V_{n}) > −∞ for all x. It remains to obtain a bound independent of n. For any n we have

K_{n}h_{n} ≤ h_{n} − K_{n}c_{n} + J̄ ≤ h_{n} + J̄1_{S}

Letting s = δ1_{S} we then obtain, for some measure ν,

(K_{n} − s ⊗ ν)h_{n} ≤ h_{n} + J̄δ^{−1}s,

and by iteration, for any N,

−ν(V_{n})(K_{n} − s ⊗ ν)^{N}1 ≤ (K_{n} − s ⊗ ν)^{N}h_{n} ≤ h_{n} + J̄δ^{−1} ∑_{i=0}^{N−1}(K_{n} − s ⊗ ν)^{i}s.

By c_{n}-regularity of the nth chain one can show that for any x,

∑_{i=0}^{∞}(K_{n} − s ⊗ ν)^{i}c (x) < ∞.

Since c ≥ 1 it follows that (K_{n} − s ⊗ ν)^{N}1 (x) → 0 as N → ∞ for any x. This then gives the bound

0 ≤ h_{n} + J̄δ^{−1}H_{n}s =^{Δ} J̄δ^{−1} ∑_{i=0}^{∞}(K_{n} − s ⊗ ν)^{i}s.

The proof is completed on noting that ∑_{i=0}^{∞}(K_{n} −s⊗ν)^{i}s (x) ≤ (δ(1−δ))^{−1} for all x (this can be seen by noting that this sum is proportional to the hitting probability to an atom for the split chain as in ^{\underline{num84,MT}} [44, 39]).

Convergence of the algorithm is subtle. This is not surprising since it is rare in optimization to prove global convergence of successive approximation. The countable state space case is considered in ^{\underline{chemey99a}} [10] where it is shown that (A1), (A2), and a strengthening of (A3) do imply convergence of {h_{n}} to a solution of the ACOE. To generalize this result to general state spaces it may be necessary to impose a blanket stability condition as in ^{\underline{hor77}} [28], or the stronger stability assumption imposed in ^{\underline{dek85a,spi90a}} [15, 55].

1.4.3 Policy iteration

s:pia

The PIA, which is again a recursive algorithm for generating useful policies, follows naturally as a refinement of the VIA. Recall that the key observation in Section ^{\underline{s:via}} 1.4.2 was the drift inequality,

P_{n−1}V_{n−1} ≤ V_{n−1} − c_{n−1} + J̄_{n−1},

From this bound we discovered easily that the next policy π_{n} has cost bounded by J(π_{n}, x) ≤ J̄_{n−1}, x ∈ 𝕏. We have seen that there are an infinite number of solutions to the drift inequality (^{\underline{e:cFoster}} 1.9), and some give better bounds than others. The optimal solution is the solution to Poisson’s equation, since this gives the minimal possible value for J̄. If the function V_{n−1} is replaced by the solution to Poisson’s equation in the VIA recursion then one obtains precisely the PIA.

To give a precise description of the algorithm, suppose that at the (n − 1)th stage of the algorithm a policy π_{n−1} is given, and assume that h_{n−1} satisfies the Poisson equation

P_{n−1}h_{n−1} = h_{n−1} − c_{n−1} + J_{n−1},

where P_{n−1} = P_{π_{n−1}}, c_{n−1}(x) = c_{π_{n−1}}(x) = c(x, π_{n−1}(x)), and J_{n−1} is a constant (equal to the steady state cost with this policy).

Given h_{n−1}, one then attempts to find an improved policy π_{n} by choosing, for each x,

π_{n}(x) = arg min_{a∈𝔸(x)}[c(x, a) + P_{a}h_{n−1} (x)].  (1.30) e:fn-def

Once π_{n} is found, polices π_{n+1}, π_{n+2}, … may be computed by induction, so long as the appropriate Poisson equation may be solved, and the minimization above has a solution.

Recall that the relative value functions {h_{n}} are not uniquely defined: If h_{n} satisfies Poisson’s equation, then so does h_{n} + b for any constant b.

The main results to follow all apply to specific normalized versions of the relative value functions.

We find that the PIA is more easily analysed then the VIA since we know a great deal about the structure of solutions to Poisson’s equation. To begin, consider the pair of equations

P_{n}h_{n} = h_{n} − c̄_{n};  (1.31) e:fishn
P_{n}h_{n−1} = h_{n−1} − c̄_{n} + γ_{n},  (1.32) e:almostFish

where c̄_{n} = c_{n}−J_{n}, and γ_{n} is defined through (^{\underline{e:almostFish}} 1.32). From the minimization (^{\underline{e:fn-def}} 1.30) we have

c_{n} + P_{n}h_{n−1} ≤ c_{n−1} + P_{n−1}h_{n−1},

and from Poisson’s equation we have

c_{n−1} + P_{n−1}h_{n−1} = h_{n−1} + J_{n−1}.

Combining these two equations gives the upper bound γ_{n}(x) ≤ J_{n−1} − J_{n}, x ∈ 𝕏, which shows that the PIA automatically generates solutions to (^{\underline{e:cFoster}} 1.9), provided the functions {h_{n}} can be taken to be positive. Through these observations we will show that the sequence {J_{n}} converges monotonically, which shows that the limit superior of {γ_{n}} is bounded from above by zero. Under suitable conditions we also show that the sequence {h_{n}} is bounded from below, and this then gives a lower bound on {γ_{n}}: Since μ_{π_{n}}(c_{n}) = J_{n}, it follows from the Comparison Theorem ^{\underline{MT}} [39, p. 337] that μ_{π_{n}}(γ_{n}) ≥ 0.

Thus, for large n, the error term γ_{n} is small, and hence the function h_{n−1} almost solves the Poisson equation for P_{n}. One might then expect that h_{n} will be close to h_{n−1}. Under mild conditions, this is shown to be true in a very strong sense.

As was the case with the value iteration algorithm, much of the analysis of ^{\underline{mey97b}} [42] focuses on \{K_{n}\} rather than \{P_{n}\}, as given in (^{\underline{e:Kn-def}} 1.23). To invoke the algorithm we must again ensure that the required minimum exists. Condition (A2) is also identical to our earlier assumption on the one step cost in Section ^{\underline{s:via}} 1.4.2.

(A1) For each n, if the PIA yields a triplet (π_{n−1}, h_{n−1}, J_{n−1}) which solve Poisson’s equation

P_{n−1}h_{n−1} = h_{n−1} − c_{n−1} + J_{n−1},

with h_{n−1} bounded from below, then the minimization

π_{n}(x) =^{Δ} arg min_{a∈𝔸(x)}[c(x, a) + P_{a}h_{n−1} (x)]

admits a measurable solution π_{n}.

(A2) For each fixed x, the function c(x, ·) is norm-like on 𝔸, and there exists a norm-like function \underline{c}: 𝕏 → ℝ_{+} such that for the policies π_{n} obtained through the PIA,

∞ > K_{n}c_{n} (x) ≥ \underline{c}(x), for any x ∈ 𝕏, n ∈ ℕ.

Under Assumptions (A1) and (A2), the algorithm produces stabilizing policies recursively. The proof is identical to Theorem ^{\underline{t:VIA-c-regular}} 7.

t:PIA-c-regular Theorem 8 (Meyn ^{\underline{mey97b}} [42]) Suppose that (A1) and (A2) hold, and that for some n, the policies {π_{i} : i < n} and relative value functions {h_{i} : i < n} are defined through the policy iteration algorithm. Suppose moreover that

(a) The relative value function h_{i} is bounded from below, i ≤ n − 1.

(b) All compact sets are petite for the Markov chains {x^{π_{i}}, i ≤ n − 1}, and for x^{π_{n}}, where π_{n} is a policy given in (A1).

Then, the PIA admits a solution (π_{n}, h_{n}, J_{n}) such that

(i) The relative value function h_{n} is bounded from below;

(ii) Each of the policies {π_{i} : i ≤ n} is regular.

(iii) For all 0 ≤ i ≤ n the constant J_{i} is the cost at the ith stage: J_{i} = J(π_{i}), and the costs are decreasing:

J_{0} ≥ J_{1} ≥ ⋯ ≥ J_{n};

To obtain convergence of the algorithm, we strengthen assumption (b) of Theorem ^{\underline{t:PIA-c-regular}} 8 to the following uniform accessibility condition:

(A3) There is a fixed probability ν on F, a δ > 0, and an initial regular policy π_{0} with the following property: For each n ≥ 1, if the PIA yields a triplet (π_{n−1}, h_{n−1}, J_{n−1}) with h_{n−1} bounded from below, then for any policy π_{n} given in (A1),

K_{n}(x, A) ≥ δν(A)  for all x ∈ S, A ∈ F,  (1.33) e:nuAccessibility

where S denotes the precompact set

S = \{x : \underline{c}(x) ≤ 2J_{0}\}.  (1.34) e:S-def-PIA

If the assumptions of Theorem ^{\underline{t:PIA-c-regular}} 8 hold, so that the relative value functions are bounded from below, then the minimal solution h_{n} to Poisson’s equation (^{\underline{e:fish}} 1.12) which is bounded from below, and which satisfies ν(h_{n}) = 0, is defined through the potential kernel in equation (^{\underline{e:fishFormula}} 1.13) of Theorem ^{\underline{t:fish}} 4. These particular versions are in fact convergent:

t:Kh-converge Theorem 9 (Meyn ^{\underline{mey97b}} [42]) Suppose that (A1)-(A3) hold, and that the initial policy π_{0} is regular. Then for each n the PIA admits a solution (π_{n}, h_{n}, J_{n}) such that π_{n} is regular, and the sequence of relative value functions {h_{n}} defined in (^{\underline{e:fishFormula}} 1.13) satisfy,

(i) For some constant N < ∞,

inf_{x∈𝕏,n≥0}h_{n}(x) > −N;

(ii) There exists a sequence of functions {g_{n} : n ≥ 0} such that

g_{n}(x) ≤ g_{n−1}(x) ≤ ⋯ ≤ g_{0}(x), x ∈ 𝕏, n ≥ 0,

and for some sequence of positive numbers {α_{k}, β_{k}},

g_{n}(x) = α_{n}h_{n}(x) + β_{n}, n ≥ 0, x ∈ 𝕏,

with α_{k} ↓ 1, β_{k} ↓ 0 as k → ∞.

These properties together imply that the relative value functions are pointwise convergent to the function h(x) =^{Δ} lim_{n} g_{n}(x).

We have already remarked that the average cost optimal control problem is plagued with counterexamples. It is of some interest then to see why Theorem ^{\underline{t:Kh-converge}} 9 does not fall into any of these traps. Consider first counterexamples 1 and 2 of ^{\underline{ros92}} [50, p. 142]. In each of these examples the process, for any policy, is completely non-irreducible in the sense that ℙ(x_{t}^{π} < x_{0}^{π}) = 0 for all times t, and all policies π. It is clear then from the cost structure that the bound (^{\underline{e:nuAccessibility}} 1.33) on the resolvent cannot hold. A third example is given in the Appendix of ^{\underline{ros92}} [50]. Here (^{\underline{e:nuAccessibility}} 1.33) is directly assumed! However, the cost is not unbounded, and is in fact designed to favor large states.

The assumptions (A2) and (A3) together imply that the center of the state space, as measured by the cost criterion, possesses some minimal amount of irreducibility, at least for the policies {π_{n}}. If either the unboundedness condition or the accessibility condition is relaxed, so that the process is non-irreducible on a set where the cost is low, then we see from these counterexamples that optimal stationary policies may not exist.

Now that we know that {h_{n}} is pointwise convergent to a function h, we can show that the PIA yields a solution to the ACOE. Theorem ^{\underline{t:PIA-opt}} 10 is similar to Theorem 4.3 of ^{\underline{lerlas96}} [26] which also requires a continuity condition related to (A4). Weaker conditions are surely possible for a specific application.

(A4) The function c: 𝕏 × 𝔸 → [1, ∞) is continuous, and the functions (P_{a}h_{n} (x) : n ≥ 0) and P_{a}h (x) are continuous in a for any fixed x ∈ 𝕏.

t:PIA-opt Theorem 10 (Meyn ^{\underline{mey97b}} [42]) Suppose that (A1)-(A3) hold, and that the initial policy π_{0} is regular. Then,

(i) The PIA produces a sequence of solutions (π_{n}, h_{n}, J_{n}) such that {h_{n}} is pointwise convergent to a solution h of the optimality equation (^{\underline{e:OE1}} 1.18), and any policy π which is a pointwise limit of π_{n} satisfies (^{\underline{e:OE2}} 1.19). Moreover, the costs {J_{n}} are decreasing with n.

(ii) Any limiting policy π is c_{π}-regular, so that for any initial condition x ∈ 𝕏,

J(π, x) = lim_{n→∞} \frac{1}{n} ∑_{t=1}^{n} 𝔼_{x}[c_{π}(x_{t}^{π})] = μ_{π}(c_{π}).

Observe that Theorem ^{\underline{t:PIA-opt}} 10 still does not claim that the solution to (^{\underline{e:OE1}} 1.18) gives rise to an optimal policy. For this one must impose extra conditions. For instance, the limiting policy will be optimal provided the relative value function h satisfies the conditions of Theorem ^{\underline{t:OE}} 5 or Theorem ^{\underline{t:comparable}} 6. Because of the monotone convergence, this will be implied by a bound on the initial condition (π_{0}, h_{0}).

We now illustrate the theory with some general examples.

1.5 Linear models

s:linearEx

The controlled linear state space model is defined through the recursion,

x_{t+1} = Ax_{t} + Ba_{t} + w_{t+1},  t ∈ ℕ,  (1.35) e:linear

where w_{t}, x_{t} ∈ ℝ^{d}, and a_{t} ∈ ℝ^{p}. The cost c in the linear-quadratic control problem takes the form

c(x, w) = \frac{1}{2}x^{T} Qx + \frac{1}{2}w^{T} Rw,  (1.36) e:lqg-cost

with Q ≥ 0, and R > 0. If w is i.i.d., then this is a Markov decision process with transition function

P_{a}(x, C) = ℙ(w_{1} + Ax + Ba ∈ C).

The optimization of J(π, x) is known as the LQG (linear-quadratic-Gaussian) problem in the special case where w is Gaussian white noise. Note that the assumption c ≥ 1 fails in this example. However, c is positive, so that we can add 1 to the cost function to satisfy the desired lower bound on c and the MDP is essentially unchanged.

In the Gaussian case one may obtain a solution (π_{∗}, h_{∗}, J_{∗}) to the ACOE with h_{∗} quadratic, and π_{∗} linear in x, by solving a Riccati equation. Why then should the solution give rise to an optimal policy? This question can be answered using Theorem ^{\underline{t:comparable}} 6: Suppose that π is any (measurable) nonlinear feedback control. From the assumption that Σ > 0, the process x^{π} is ψirreducible, with ψ = Lebesgue measure. The function V_{π} given in (^{\underline{e:optV}} 1.20) then satisfies the lower bound

V_{π}(x) ≥ ∑_{t=0}^{∞} 2^{−k}P_{π}^{t}c_{π} (x) ≥ V_{\frac{1}{2}} (x) − b_{π},  (1.37) e:lqg-A4

where V_{\frac{1}{2}} is the value function for the discounted problem, and b_{π} is a finite constant. To see this, let τ_{π} = τ_{S_{π}}, and write

V_{\frac{1}{2}} (x) ≤ ∑_{t=0}^{∞} 2^{−t}P_{π}^{t}c_{π} (x) = 𝔼_{x}[ ∑_{t=0}^{τ_{π}−1} 2^{−k}c_{π}(x_{t}^{π})] + 𝔼_{x}[∑_{t=τ_{π}}^{∞} 2^{−k}c_{π}(x_{t}^{π})]

≤ V_{π}(x) + 𝔼_{x}[2^{−τ_{π}} ∑_{t=0}^{∞} 2^{−k}c_{π}(x_{τ_{π}+t}^{π})]

The lower bound (^{\underline{e:lqg-A4}} 1.37) on V_{π} then follows from the strong Markov property, and regularity of the set S_{π}. If (A, \sqrt{Q}) is observable so that V_{\frac{1}{2}} dominates a positive definite quadratic, then we see from Theorem ^{\underline{t:comparable}} 6 that the linear/quadratic solution (π_{∗}, h_{∗}) to the ACOE does yield an optimal policy over the class of all nonlinear feedback control laws.

The linear state space model gives an excellent test-case for interpreting the assumptions of Theorems ^{\underline{t:VIA-c-regular}} 7–^{\underline{t:Kh-converge}} 9. Condition (A1), which demands the existence of a minimizing policy, is satisfied because the model is continuous. The assumption of an initial regular policy is simply stabilizability of (A, B). The norm-like condition (A2) is implied by the standard observability condition on (A, \sqrt{Q}), where Q is the state weighting matrix given in (^{\underline{e:lqg-cost}} 1.36)

\underline{mey97b} [42].

To verify the uniform accessibility condition (^{\underline{e:nuAccessibilityVIA}} 1.27) in (A3), suppose the initial policy π_{0} is linear, so that each subsequent policy is of the form π_{n}(x) = −K_{n}x, and let A_{n} = A−BK_{n} denote the closed loop system matrix. Then the accessibility condition (A3) holds if Σ > 0, and the steady state costs {J_{n}} are bounded. For both the PIA and the VIA, the boundedness condition holds automatically under (A1) and (A2) through stability of the algorithm, which is guaranteed by Theorems ^{\underline{t:VIA-c-regular}} 7 and ^{\underline{t:PIA-c-regular}} 8.

From these results it follows that Theorem ^{\underline{t:Kh-converge}} 9 recovers known properties of the Newton-Raphson technique applied to the LQG problem. Consider the well known decreasing property of the solutions {Λ_{n}} to the associated Riccati equation. The proof of Theorem ^{\underline{t:Kh-converge}} 9 depends upon the bound

h_{n}(x) ≤ [1 + 2(J_{n−1} − J_{n})]h_{n−1}(x) + b_{n},  (1.38) e:h-bdds2

where b_{n} is a constant. In the case of linear controls, it can be shown that the relative value function h_{n} takes the form h_{n}(x) = h_{n}(0) + x^{T} Λ_{n}x. Letting x → ∞, it follows from the previous inequality that

Λ_{n} ≤ [1 + 2(J_{n−1} − J_{n})]Λ_{n−1},  n ≥ 1.  (1.39) e:Pdecrease

It may be shown directly that Λ_{n} ≤ Λ_{n−1} ^{\underline{duf90,whi90}} [19, 60], so the bound (^{\underline{e:h-bdds2}} 1.38) is not tight in the linear model. However, the semi-decreasing property (^{\underline{e:Pdecrease}} 1.39) is sufficient to deduce convergence of the algorithm.

There is no space here to consider the VIA in further detail. We note however that it is well known that the successive approximation procedure generates stabilizing policies for the linear state space model provided the initial policy is stabilizing and linear. Theorem ^{\underline{t:VIA-c-regular}} 7 shows that it is enough to assume only stability.

1.6 Network models

s:netEx

We now apply the general results of Section ^{\underline{s:opt}} 1.4 to the scheduling problem for multiclass queueing networks. For simplicity we discuss here only a relatively simple class of network models which can be formulated through an extension of the M/M/1 model. A treatment of general network models is given in ^{\underline{mey99a}} [43].

Consider a network composed of d single server stations, indexed by σ = 1, …, d. The network is populated by ℓ classes of customers: Class k customers require service at station s(k). An exogenous stream of customers of class 1 arrive to machine s(1), and subsequent routing of customers is deterministic. If the service times and interarrival times are assumed to be exponentially distributed, then after a suitable time scaling and sampling of the process, the dynamics of the network can be described by the random linear system,

x_{t+1} = x_{t} + ∑_{k=0}^{ℓ} I_{t+1}(k)[e^{k+1} − e^{k}]a_{t}(k),  (1.40) e:Net

where the state process x evolves on 𝕏 = ℕ^{ℓ}, and x_{t}(k) denotes the number of class k customers in the system at time t. An example of a two station network is illustrated in Figure ^{\underline{f:net}} 1.3.

The random variables {I_{n} : n ≥ 0} are i.i.d. on {0, 1}^{ℓ+1}, with

ℙ{∑_{i} I_{n}(k) = 1} = 1, and 𝔼[I_{n}(k)] = μ_{k}.

For 1 ≤ k ≤ ℓ, μ_{k} denotes the service rate for class k customers. For k = 0, we let μ_{0} =^{Δ} λ denote the arrival rate of customers of class 1. For 1 ≤ k ≤ ℓ we let e^{k} denote the kth basis vector in ℝ^{ℓ}, and we set e^{0} = e^{ℓ+1} =^{Δ} 0.

The sequence {a_{t} : t ≥ 0} is the control, which takes values in {0, 1}^{ℓ+1}. We define a_{t}(0) ≡ 1. The set of admissible control actions 𝔸(x) is defined in an obvious manner: for a ∈ 𝔸(x),

(i) For any 1 ≤ k ≤ ℓ, a_{k} = 0 or 1;

(ii) For any 1 ≤ k ≤ ℓ, x_{k} = 0 ⇒ a_{k} = 0;

(iii) For any station σ, 0 ≤ ∑_{k:s(k)=σ} a_{k} ≤ 1;

(iv) For any station σ, ∑_{k:s(k)=σ} a_{k} = 1 whenever ∑_{k:s(k)=σ} x_{k} > 0.

If a_{k} = 1, then buffer k is chosen for service. Condition (ii) then imposes the physical constraint that a customer cannot be serviced at a buffer if that buffer is empty. Condition (iii) means that only one customer may be served at a given instant at a single machine σ. The non-idling condition (iv) is satisfied by any optimal policy with this cost criterion: An inductive proof is given in ^{\underline{mey97a}} [41] based upon value iteration.

Since the control is bounded, a reasonable cost function is c(x, a) = c^{T} x, where c ∈ ℝ^{ℓ} is a vector with strictly positive entries. For concreteness, we take c(x, a) = |x| =^{Δ} ∑_{k} x(k). Since 𝔸(x) is a finite set for any x, it follows that (A1) holds with this cost function.

Figure 1.3: A multiclass network with d = 2 and ℓ = 3.  f:net

The transition function has the simple form,

ℙ_{a}(x, x + e^{k+1} − e^{k}) = μ_{k}a_{k},  0 ≤ k ≤ ℓ.
ℙ_{a}(x, x) = 1 − ∑_{0}^{ℓ} μ_{k}a_{k}

The accessibility condition (^{\underline{e:accessible}} 1.4) holds with s everywhere positive, and ν = δ_{θ}, with θ equal to the empty state θ = (0, …, 0)^{T} ∈ 𝕏. This follows from the non-idling assumption.

Associated with this network is a fluid model. For each initial condition x_{0} ≠ 0, we construct a continuous time process ϕ^{x}(t) as follows. If m = |x_{0}|, and if tm is an integer, we set

ϕ^{x}(t) = \frac{1}{m}x_{mt}.

For all other t ≥ 0, we define ϕ^{x}(t) by linear interpolation, so that it is continuous and piecewise linear in t. Note that |ϕ^{x}(0)| = 1, and that ϕ^{x} is Lipschitz continuous. The collection of all “fluid limits” is defined by

L =^{Δ} ⋂_{n=1}^{∞} \overline{\{ϕ^{x} : |x| > n\}}

where the overbar denotes weak closure. This set of stochastic process of course depends on the particular policy π which has been applied. The process ϕ evolves on the state space ℝ_{+}^{ℓ} and, for a wide class of scheduling policies, satisfies a differential equation of the form

\frac{d}{dt}ϕ(t) = ∑_{k=0}^{ℓ} μ_{k}[e^{k+1} − e^{k}]u_{t}(k),  (1.41) e:fluidNet

where the function u_{t} is analogous to the discrete control, and satisfies similar constraints (see the M/M/1 queue model described earlier, or ^{\underline{dai95a,daiwei96}} [13, 12] for more general examples).

Stability of (^{\underline{e:Net}} 1.40) in terms of c-regularity is closely connected with the stability of the fluid model ^{\underline{dai95a,kummey96a,daimey95a}} [13, 35, 14]. The fluid model L is called L_{p}-stable if

lim_{t→∞} sup_{ϕ∈L} 𝔼[|ϕ(t)|^{p}] = 0.

It is shown in ^{\underline{kummey96a}} [35] that L_{2}-stability of the fluid model is equivalent to a form of c-regularity for the network.

t:equi Theorem 11 (Kumar and Meyn ^{\underline{kummey96a}} [35]) The following stability criteria are equivalent for the network under any nonidling policy.

(i) The drift condition (^{\underline{e:cFoster}} 1.9) holds for some function V. The function V is equivalent to a quadratic in the sense that, for some γ > 0,

1 + γ|x|^{2} ≤ V (x) ≤ 1 + γ^{−1}|x|^{2},  x ∈ 𝕏.  (1.42) e:QuadEquiv

(ii) For some quadratic function V,

𝔼_{x}[∑_{n=0}^{σ_{θ}} |x_{n}|] ≤ V (x),  x ∈ 𝕏,

where σ_{θ} is the first entrance time to θ = 0.

(iii) For some quadratic function V and some J̄ < ∞,

∑_{n=1}^{N} 𝔼_{x}[|x_{n}|] ≤ V (x) + NJ̄, for all x and N ≥ 1.

(iv) The fluid model L is L_{2}-stable.

The previous result can be strengthened: if the fluid model is L_{2}-stable, then in fact ϕ(t) = 0 for all t sufficiently large ^{\underline{mey97a}} [41].

Using this result it is shown in ^{\underline{mey97b}} [42] that when applying policy iteration to a network model, on performing the fluid scaling one obtains a sequence of fluid models which are the solutions of a policy iteration scheme for the fluid model. Moreover, the algorithm convergence to yield a policy which is optimal for both the network and its fluid model.

Let us turn to the VIA: In view of Theorem ^{\underline{t:equi}} 11, how should we initialize the algorithm? Two possibilities are suggested:

Figure 1.4: Convergence of the VIA with V_{0} taken as the value function for the associated fluid control problem, or a pure quadratic function obtained

through a linear program.  f:Fluid

(i) Given the previous analysis of the M/M/1 queue it appears natural to set V_{0} equal to the value function for a fluid model,

V_{0}(x) = min ∫_{0}^{∞} |ϕ(t)| dt ϕ(0) = x, x ∈ 𝕏,

where the minimum is with respect to all policies for the fluid model. One can show that for large x, V_{0} does approximate the relative value

function ^{\underline{mey97a,mey97b,henmey99a}} [41, 42, 25].

(ii) The conclusion that the relative value function is ‘nearly quadratic’ suggests that we search for a pure quadratic form satisfying (^{\underline{e:cFoster}} 1.9),

V_{0}(x) = x^{T} Qx, x ∈ 𝕏.

In ^{\underline{kummey96a}} [35] a linear program is constructed to compute a quadratic solution to (^{\underline{e:cFoster}} 1.9) for network models, based on prior results of ^{\underline{kumkum94a,berpastsi92a}} [36, 46].

We conclude with a numerical experiment to show how a careful initialization can dramatically speed convergence of the VIA. We consider the three buffer model illustrated in Figure ^{\underline{f:net}} 1.3 with the following parameters: λ/μ_{2} = 9/10; λ/μ_{1} + λ/μ_{3} = 9/11; and μ_{1} = μ_{3}. The optimal value function can be computed explicitly in this case, and a pure quadratic Lyapunov function can also be computed easily.

Two experiments were performed to compare the performance of the VIA initialized with these two value functions. To apply value iteration the buffer levels were truncated so that x_{i} < 45 for all i. This gives rise to a finite state space MDP with 45^{3} = 91, 125 states. The results from two experiments are shown in Figure ^{\underline{f:Fluid}} 1.4. For comparison, data from the standard VIA with V_{0} ≡ 0 is also given. We have taken 300 steps of value iteration, saving data for n = 10, …, 300. The convergence is exceptionally fast in both experiments. Note that the convergence of J_{n} is not monotone in the experiment shown using the fluid value function initialization. However, this initialization leads to fast convergence to the optimal cost J_{∗} ≈ 10.9.

1.7 Extensions and Open Problems

It is hoped that the development in this chapter has suggested to the reader some interesting topics for further research. We list here some areas which have been of interest to the author.

Continuous time There is very little of interest to say in the continuous time framework unless one specializes to an interesting class of models. In this chapter the analysis has been restricted to a resolvent kernel, and the same approach can be followed in continuous time where the resolvent becomes

R_{β} = ∫_{0}^{∞} βe^{−βt}P^{t} dt

with β > 0. Again one can show that any variable of interest (the invariant measures, solutions to Poisson’s equation, or solutions to (^{\underline{e:cFoster}} 1.9)) can be mapped between the resolent and the continuous time process. Further discussion may be found in ^{\underline{mey97b,sen99a}} [42, 53].

Geometric ergodicity and risk sensitive control The risk sensitive control criterion is given via

J_{γ}(π, x) =^{Δ} lim sup_{n→∞} \frac{1}{n} log(𝔼_{x}[exp(γ ∑_{t=0}^{n−1} c_{π}(x_{t}^{π}))]).

where the ‘risk factor’ γ is typically assumed to be a small, positive number. Models of this sort were first considered in ^{\underline{bel57}} [3, p. 329] for finite state space models, and complete solutions to the optimization problem were given in \underline{howmat72a,rot84a} [32, 51], again in the finite state space case. This control problem has attracted more recent attention because of the interesting connections between risk sensitive control and game theory ^{\underline{jac73,flemce92,whi90}} [33, 21, 60].

Under a norm-like condition on the model it can be shown that when this cost is finite valued, the Markov chain exhibits a strong form of stability known as geometric ergodicity ^{\underline{balmey98a,bormey98a}} [2, 6]. Conversely, such stability assumptions imply that the cost is finite, and ensure that an optimal policy does exist

\underline{hermar96,cavfer98a} [23, 9].

Our present understanding of the optimization problem for Markov chains on an infinite state space is currently very weak, but this appears to be an area worthy of further study.

Simulation The use of simulation will become increasingly important in both evaluating and synthesizing policies. Much of the burden of finding an optimal policy surrounds the solution of Poisson’s equation, for which now there are several simulation based algorithms such as temporal difference learning, and there are also simulation based versions of both value and policy iteration (see ^{\underline{bertsi96a}} [4]).

Complexity This has always been one of the most challenging issues in optimal control. Markovian models are frequently too ‘fine grained’ to be useful in optimization. One solution then is to seek some form of aggregation. For general MDP models one can directly discretize the state space to obtain a finite state space model. In the case of network models, either fluid models or Brownian motion models provide approaches to aggregation which deserve further study.

Bibliography

araborferghomar93 [1] A. Arapostathis, V. S. Borkar, E. Fernandez-Gaucherand, M. K. Ghosh, and S. I. Marcus. Discrete-time controlled Markov processes with average cost criterion: a survey. SIAM J. Control Optim., 31:282– 344, 1993.

balmey98a [2] S. Balaji and S.P. Meyn. Multiplicative ergodic theorems for an irreducible Markov chain. Stoch. Processes and their Appl., 1999. Submitted for publication.

bel57 [3] R. Bellman. Dynamic Programming. Princeton University Press, Princeton, NJ, 1957.

bertsi96a [4] Bertsekas, D., Tsitsiklis, J. Neuro-Dynamic Programming, Athena Scientific, Belmont, MA, 1996.

bor91 [5] V. S. Borkar. Topics in controlled Markov chains. Pitman Research Notes in Mathematics Series # 240, Longman Scientific & Technical, UK, 1991.

bormey98a [6] V. Borkar and S.P. Meyn. Risk Sensitive Optimal Control: Existence and Synthesis for Models with Unbounded Cost. Mathematics of O.R., 1999. Submitted for publication.

cav96 [7] R. Cavazos-Cadena. Value iteration in a class of communicating Markov decision chains with the average cost criterion. Technical report, Universidad Autónoma Agraria Anonio Narro, 1996.

cavfer95c [8] R. Cavazos-Cadena and E. Fernandez-Gaucherand. Value iteration in a class of average controlled Markov chains with unbounded costs: Necessary and sufficient conditions for pointwise convergence. J. Applied Probability, 33:986–1002, 1996. cavfer98a [9] R. Cavazos-Cadena and E. Fernandez-Gaucherand. Controlled Markov

chains with risk-sensitive criteria: Average cost, optimality equations, and optimal solutions. Mathematical Methods of Operations Research (1999) 49:299-324.

chemey99a [10] R-R. Chen and S.P. Meyn. Value iteration and optimization of multiclass queueing networks. to appear, QUESTA, 1999.

enghummey96a [11] J. Humphrey D. Eng and S.P. Meyn. Fluid network models: Linear programs for control and performance bounds. In J. Cruz J. Gertler and M. Peshkin, editors, Proceedings of the 13th IFAC World Congress, volume B, pages 19–24, San Francisco, California, 1996.

daiwei96 [12] J. Dai and G. Weiss. Stability and instability of fluid models for certain re-entrant lines. Mathematics of Operations Research, 21(1):115–134, February 1996.

dai95a [13] J. G. Dai. On the positive Harris recurrence for multiclass queueing networks: A unified approach via fluid limit models. Ann. Appl. Probab., 5:49–77, 1995.

daimey95a [14] J. G. Dai and S.P. Meyn. Stability and convergence of moments for multiclass queueing networks via fluid limit models. IEEE Trans. Automat. Control, 40:1889–1904, November 1995.

dek85a [15] R. Dekker. Denumerable Markov Decision Chains: Optimal Policies for Small Interest Rates, PhD thesis, University of Leiden, Leiden, the Netherlands, 1985.

dek87 [16] R. Dekker. Counterexamples for compact action Markov decision chains with average reward criteria. Comm. Statist.-Stoch. Models, 3:357–368, 1987.

der66 [17] C. Derman. Dunumerable state MDPs. Ann. Amth. Statist., 37:1545– 1554, 1966.

dowmeytwe96a [18] D. Down, S. P. Meyn, and R. L. Tweedie. Geometric and uniform ergodicity of Markov processes. Ann. Probab., 23(4):1671–1691, 1996.

duf90 [19] M. Duflo. Méthodes Récursives Aléatoires. Masson, 1990. dynyus79 [20] E. B. Dynkin and A. A. Yushkevich. Controlled Markov Processes,

volume Grundlehren der mathematischen Wissenschaften 235 of A Series of Comprehensive Studies in Mathematics. Springer-Verlag, New York, NY, 1979.

flemce92 [21] W.H. Fleming and W.M. McEneaney. Risk-sensitive control and differential games. volume 84 of Lecture Notes in Control and Info. Sciences, pages 185–197. Springer-Verlag, Berlin; New York, 1992.

glymey96a [22] P. W. Glynn and S. P. Meyn. A Lyapunov bound for solutions of Poisson’s equation. Ann. Probab., 24, April 1996.

hermar96 [23] D. Hernández-Hernández and S.I. Marcus. Risk sensitive control of Markov processes in countable state space. Systems Control Lett., 29:147–155, July 1996. correction in Systems and Control Letters, 34(1-2), 1998, pp. 105-106.

hen97a [24] Henderson, S. G. Variance Reduction Via an Approximating Markov Process. Ph.D. thesis. Department of Operations Research, Stanford University. Stanford, California, USA, 1997.

henmey99a [25] S.G. Henderson and S.P. Meyn. Variance reduction for simulation in multiclass queueing networks. submitted to the IIE Transactions on Operations Engineering: special issue honoring Alan Pritsker on simulation in industrial engineering, 1999.

lerlas96 [26] O. Hernández-Lerma and J. B. Lasserre. Discrete time Markov control processes I. Springer-Verlag, New York, 1996.

hermoncav91 [27] O. Hernández-Lerma, R. Montes-de-Oca, and R. Cavazos-Cadena. Recurrence conditions for Markov decision processes with Borel state space: A survey. Ann. Operations Res., 28:29–46, 1991.

hor77 [28] A. Hordijk. Dynamic Programming and Markov Potential Theory. 1977.

horput87 [29] A. Hordijk and M. L. Puterman. On the convergence of policy iteration. Math. Op. Res., 12:163–176, 1987.

horspitwe95 [30] A. Hordijk, F. M. Spieksma, and R. L. Tweedie. Uniform stability conditions for general space Markov decision processes. Technical report, Leiden University and Colorado State University, 1995. Technical Report. how60 [31] R. A. Howard. Dynamic Programming and Markov Processes. John

Wiley and Sons/MIT Press, New York, NY, 1960.

howmat72a [32] R.A. Howard and J.E. Matheson. Risk-sensitive Markov decision processes. Management Sci., 8:356–369, 1972.

jac73 [33] D. H. Jacobson. Optimal stochastic linear systems with exponential performance criteria and their relation to deterministic differential games. IEEE Trans. Automat. Control, AC-18:124–131, 1973.

konbor96 [34] V. R. Konda and V. S. Borkar. Learning algorithms for Markov decision processes. Technical report, Indian Institute of Science, Bangalore, 1996.

kummey96a [35] P.R. Kumar and S.P. Meyn. Duality and linear programs for stability and performance analysis queueing networks and scheduling policies. IEEE Transactions on Automatic Control, 41(1):4–17, January 1996.

kumkum94a [36] S. Kumar and P. R. Kumar. Performance bounds for queueing networks and scheduling policies. IEEE Trans. Automat. Control, AC-39:1600–1611, August 1994.

kwasiv72 [37] H. Kwakernaak and R. Sivan. Linear Optimal Control Systems. Wiley-Interscience, New York, NY, 1972.

lip75 [38] S. Lippman. Applying a new device in the optimization of exponential queueing systems. Operations Research, 23:687–710, 1975.

MT [39] S. P. Meyn and R. L. Tweedie. Markov Chains and Stochastic Stability. Springer-Verlag, London, 1993.

meytwe93b [40] S. P. Meyn and R. L. Tweedie. Stability of Markovian processes III: Foster-Lyapunov criteria for continuous time processes. Adv. Appl. Probab., 25:518–548, 1993.

mey97a [41] S.P. Meyn. Stability and optimization of multiclass queueing networks and their fluid models. In proceedings of the summer seminar on “The Mathematics of Stochastic Manufacturing Systems”. American Mathematical Society, 1997.

mey97b [42] S.P. Meyn. The policy improvement algorithm for Markov decision processes with general state space. IEEE Trans. Automat. Control, AC-42:191–196, 1997. mey99a [43] S.P. Meyn. Feedback Regulation for Sequencing and Routing in Mul-

ticlass Queueing Networks. submitted for publication, 1999.

num84 [44] E. Nummelin. General Irreducible Markov Chains and Non-Negative Operators. Cambridge University Press, Cambridge, 1984.

num91 [45] E. Nummelin. On the Poisson equation in the potential theory of a single kernel. Math. Scand., 68:59–82, 1991.

berpastsi92a [46] I. Ch. Paschalidis D. Bertsimas and J. N. Tsitsiklis. Scheduling of multiclass queueing networks: Bounds on achievable performance. In Workshop on Hierarchical Control for Real–Time Scheduling of Manufacturing Systems, Lincoln, New Hampshire, October 16–18, 1992.

per93a [47] J. Perkins. Control of Push and Pull Manufacturing Systems. PhD thesis, University of Illinois, Urbana, IL, September 1993. Technical report no. UILU-ENG-93-2237 (DC-155).

put94 [48] M. L. Puterman. Markov Decision Processes. Wiley, New York, 1994.

ritsen93 [49] R. K. Ritt and L. I. Sennott. Optimal stationary policies in general state space Markov decision chains with finite action set. Mathematics of Operations Research, 17(4):901–909, November 1993.

ros92 [50] S. M. Ross. Applied probability models with optimization applications. Dover books on advanced Mathematics, 1992. Republication of the work first published by Holden-Day, 1970.

rot84a [51] U.G. Rothblum. Multiplicative Markov decision chains. Math. Operations Res., 9:6–24, 1984.

sen86 [52] L. I. Sennott. A new condition for the existence of optimal stationary policies in average cost Markov decision processes. Operations Research Letters, 5:17–23, 1986.

sen99a [53] L.I. Sennott. Stochastic Dynamic Programming and the Control of Queueing Systems. Wiley, 1999.

sen96a [54] L.I. Sennott. The convergence of value iteration in average cost Markov decision chains. Operations Research Letters, 19:11–16, 1996.

spi90a [55] F.M. Spieksma. Geometrically Ergodic Markov Chains and the Optimal Control of Queues, PhD thesis, University of Leiden, Leiden, the Netherlands, 1990. tsivan96 [56] J. N. Tsitsiklis and B. Van Roy. An analysis of temporal-difference

learning with function approximation. technical report LIDS-P-2322, Massachusetts Institute of Technology, Cambridge, MA, March 1996. To appear in the IEEE Transactions on Automatic Control.

tuotwe94a [57] P. Tuominen and R.L. Tweedie. Subgeometric rates of convergence of f-ergodic Markov chains. Adv. Appl. Probab., 26:775–798, 1994.

websti87 [58] R. Weber and S. Stidham. Optimal control of service rates in networks of queues. Adv. Appl. Probab., 19:202–218, 1987.

wei94a [59] G. Weiss. On the optimal draining of re-entrant fluid lines. Technical report, Georgia Georgia Institute of Technology and Technion, 1994.

whi90 [60] P. Whittle. Risk-Sensitive Optimal Control. John Wiley and Sons, Chichester, NY, 1990.