Continuous-Time Markov Chains · Study tools

SimulationYin–Zhang, Ch. 1: Introduction and Overview · Ch. 2: Mathematical Preliminaries

Generators and sample paths

A continuous-time Markov chain on finitely many states is determined by one matrix, its generator $Q$ . The diagonal says how long the chain waits in each state, and the off-diagonal entries say where it goes next. Set the rates below and run the chain.

The idea

Generator. An

m\times m

matrix

Q=(q_{ij})

with

q_{ij}\ge 0

for

i\neq j

and every row summing to zero, so

q_{ii}=-\sum_{j\neq i}q_{ij}

. Write

q_i=-q_{ii}

for the total rate of leaving

i

. A time-dependent

Q(t)

must have this form at each

t

(the book also asks for measurability and bounded entries).

Construction from rates. In state

i

, wait a time that is exponential with rate

q_i

(mean

1/q_i

), then jump to

j\neq i

with probability

q_{ij}/q_i

, independently of the past. The sequence of states visited is the jump chain, a discrete-time chain with those transition probabilities. If

q_i=0

, state

i

is absorbing.

The transition matrix $P(t)=\big(P(\alpha(t)=j\mid\alpha(0)=i)\big)$ then satisfies $P(t+s)=P(t)P(s)$ and

P(\alpha(t+h)=j\mid \alpha(t)=i) = q_{ij}\,h + o(h)\quad (j\neq i), \qquad P(t)=e^{Qt}.

Martingale characterisation.

Q

generates

\alpha(\cdot)

exactly when, for every function

f

on the state space,

f(\alpha(t)) - \int_0^t (Qf)(\alpha(s))\,ds, \qquad (Qf)(i)=\sum_{j\neq i} q_{ij}\big(f(j)-f(i)\big),

is a martingale. This is the definition the book works with, and it is what makes the expected-value equations of later modules come out in one line.

Long-run time fractions. If

Q

is irreducible (every state can reach every other), then with probability one

\frac1t\int_0^t I_{\{\alpha(s)=i\}}\,ds \to \nu_i

, where

\nu

is the unique probability vector with

\nu Q = 0

.

Why two time scales. In the systems the book has in mind (machines that fail and get repaired, queues, networks, economic regimes) some rates are much larger than others. Writing the generator as $Q^\varepsilon=\tilde Q/\varepsilon+\hat Q$ , with a small $\varepsilon>0$ , separates the fast moves from the slow ones. Modules 4 to 13 are about that limit. Note that $Q/\varepsilon$ generates the same chain as $Q$ , only run $1/\varepsilon$ times faster.

Explore

Speed

Arrow width is the rate

q_{ij}

. The filled node is the current state, and its outgoing arrows are highlighted.

fraction of time so far stationary

\nu_i

. Below: the last 12 time units of the path.

Worked example

Take the generator

Q=\begin{pmatrix} -3 & 2 & 1\\ 1 & -1 & 0\\ 2 & 2 & -4\end{pmatrix}.

Describe how the chain moves, find its stationary distribution, and compare the fraction of time in each state with the fraction of jumps into it.

What are the mean holding times in states 1, 2, 3?

The leaving rates are $q_1=3,\ q_2=1,\ q_3=4$ , so the mean holding times are $1/3,\ 1,\ 1/4$ .
Write down the transition matrix of the jump chain.

Divide each off-diagonal entry by its row's leaving rate and put zeros on the diagonal:

$\Pi=\begin{pmatrix} 0 & 2/3 & 1/3\\ 1 & 0 & 0\\ 1/2 & 1/2 & 0\end{pmatrix}.$

From state 2 the chain always returns to 1.
Solve $\nu Q=0$ , $\nu_1+\nu_2+\nu_3=1$ .

The third column gives $\nu_1-4\nu_3=0$ , so $\nu_1=4\nu_3$ . The first column gives $-3\nu_1+\nu_2+2\nu_3=0$ , so $\nu_2=10\nu_3$ . The second column, $2\nu_1-\nu_2+2\nu_3=8\nu_3-10\nu_3+2\nu_3=0$ , checks out. Normalising,

$\nu=\tfrac1{15}(4,\ 10,\ 1).$
The jump chain has its own stationary distribution $\pi$ (with $\pi\Pi=\pi$ ). How is it related to $\nu$ ?

The chain enters state $i$ at long-run rate $\nu_i q_i$ (time fraction times leaving rate), so $\pi_i\propto\nu_iq_i$ :

$\pi\propto(4\cdot3,\ 10\cdot1,\ 1\cdot4)=(12,10,4),\qquad \pi=\tfrac1{13}(6,5,2).$

Check: $\pi_1=\pi_2+\tfrac12\pi_3=\tfrac{5}{13}+\tfrac1{13}=\tfrac6{13}$ . State 1 gets the most jumps but state 2 holds the most time, because the chain lingers there.

Check your understanding

Row $i$ of a generator is $(\,1,\ -4,\ 3\,)$ , with $q_{ii}=-4$ in the middle. What is the chance the next jump goes to the third state, and how long does the chain wait on average?

Probability $3/4$ ; mean wait $1/4$ .
Probability $3/4$ ; mean wait $4$ .The holding time is exponential with rate 4. Its mean is the reciprocal of the rate.
Probability $3\,h$ in a short time $h$ ; there is no single "next jump" probability. $3h$ is the chance of a jump to state 3 within time $h$ . Conditioned on a jump happening, the destination has a fixed distribution.

The two exits race as independent exponential clocks with rates 1 and 3. The first ring comes after an

\mathrm{Exp}(4)

time, and the rate-3 clock wins with probability

3/(1+3)

.

Why must every row of a generator sum to zero?

Because the stationary distribution satisfies $\nu Q=0$ . $\nu Q=0$ is a statement about the columns, weighted by $\nu$ . It is a property of one special vector, not of the rows.
Because each row of $P(t)=e^{Qt}$ sums to 1 for all $t$ ; differentiate at $t=0$ .
Because probabilities are nonnegative.Nonnegativity is why the off-diagonal entries are $\ge0$ . It says nothing about the sum.

P(t)\mathbf 1=\mathbf 1

for all

t

, and

\frac{d}{dt}P(t)\big|_{t=0}=Q

, so

Q\mathbf 1=0

. In words: the rate of leaving

i

equals the sum of the rates of going to each other state.

Chain A has generator $Q$ and chain B has generator $10\,Q$ . Which is true?

B has a different stationary distribution, ten times larger. $\nu(10Q)=10\,\nu Q=0$ : the same $\nu$ solves both systems, and it must still sum to 1.
B has a different jump chain.The jump probabilities are ratios $q_{ij}/q_i$ ; multiplying every rate by 10 cancels.
Same jump chain and same $\nu$ , but B runs ten times faster: $\alpha_B(t)$ has the law of $\alpha_A(10t)$ .

Scaling a generator rescales time and nothing else. That is the starting point of the two-time-scale theory:

Q/\varepsilon

is the chain of

Q

run at speed

1/\varepsilon

, which mixes in time of order

\varepsilon

.

Practice

The two-state chainAn explicit transition matrix

A two-state chain jumps $1\to2$ at rate $\lambda>0$ and $2\to1$ at rate $\mu>0$ . Find $P(t)=e^{Qt}$ explicitly and its limit as $t\to\infty$ .

Hint

$p_{11}(t)$ solves $p_{11}'=-\lambda p_{11}+\mu(1-p_{11})$ . Alternatively, $Q$ has eigenvalues $0$ and $-(\lambda+\mu)$ .

Solution

The equation $p_{11}'=\mu-(\lambda+\mu)p_{11}$ with $p_{11}(0)=1$ gives

P(t)=\frac1{\lambda+\mu}\begin{pmatrix}\mu+\lambda e^{-(\lambda+\mu)t} & \lambda-\lambda e^{-(\lambda+\mu)t}\\ \mu-\mu e^{-(\lambda+\mu)t} & \lambda+\mu e^{-(\lambda+\mu)t}\end{pmatrix}.

Both rows tend to $\nu=(\mu,\lambda)/(\lambda+\mu)$ , at the exponential rate $\lambda+\mu$ .

From the martingale to the forward equationTaking expectations

Let $p_j(t)=P(\alpha(t)=j)$ . Using the martingale characterisation with $f=I_{\{\cdot=j\}}$ , show that $p_j'(t)=\sum_i p_i(t)\,q_{ij}$ .

Solution

For $f(k)=I_{\{k=j\}}$ we get $(Qf)(i)=\sum_k q_{ik}f(k)=q_{ij}$ (using $\sum_k q_{ik}=0$ ). A martingale has constant expectation, so

p_j(t)-p_j(0)=\mathbb E\!\int_0^t q_{\alpha(s)j}\,ds=\int_0^t\sum_i p_i(s)q_{ij}\,ds.

Differentiating gives $p'(t)=p(t)Q$ , the forward equation of the next module.

InteractiveYin–Zhang, Ch. 2: Mathematical Preliminaries

The forward equation and irreducibility

The distribution of a Markov chain solves a linear differential equation, $p'(t)=p(t)Q$ . For three states every solution is a curve in a triangle, and you can drag its starting point. When the generator is weakly irreducible, every curve ends at the same point.

The idea

Kolmogorov equations. Write

P(t,s)

for the matrix of probabilities

P(\alpha(t)=j\mid\alpha(s)=i)

,

t\ge s

. Then

\frac{\partial}{\partial t}P(t,s)=P(t,s)\,Q(t)\ \ \text{(forward)},\qquad \frac{\partial}{\partial s}P(t,s)=-Q(s)\,P(t,s)\ \ \text{(backward)},

with

P(s,s)=I

. The row vector

p(t)=(P(\alpha(t)=1),\dots,P(\alpha(t)=m))

solves

p'(t)=p(t)Q(t)

. For constant

Q

,

P(t,s)=e^{Q(t-s)}

.

Weak irreducibility.

Q

is weakly irreducible if the system

\nu Q=0,\ \sum_i\nu_i=1

has exactly one solution. That solution is then automatically

\ge 0

.

Q

is irreducible if in addition every

\nu_i>0

. Equivalently,

Q

is weakly irreducible when

fQ=0

and

f\mathbf 1=0

force

f=0

; its rank is then

m-1

.

Quasi-stationary distribution. For a time-dependent

Q(t)

, the solution

\nu(t)

of

\nu(t)Q(t)=0,\ \nu(t)\mathbf 1=1

. It is the distribution the chain would settle into if the rates were frozen at their time-

t

values. It is not, in general, the distribution of

\alpha(t)

.

Exponential convergence. If a constant

Q

is weakly irreducible, there are

K,\kappa>0

with

\big|e^{Qt}-\mathbf 1\nu\big|\le Ke^{-\kappa t}

. Every row of

e^{Qt}

tends to

\nu

, and

\kappa

can be taken just below the smallest

|\mathrm{Re}\,\lambda|

over the nonzero eigenvalues

\lambda

of

Q

.

Fredholm alternative. If

Q

is weakly irreducible,

fQ=g

has a solution exactly when

g\mathbf 1=0

, and fixing

f\mathbf 1

makes it unique. Every expansion in Modules 4 to 6 is built by solving equations of this form.

Explore

Generator

Time t 1.00

Each point of the triangle is a distribution

(p_1,p_2,p_3)

; a corner is a sure state. Drag p(0). Dots mark

t=1,2,\dots

p_1

p_2

p_3

, with

\nu_i

dashed.

Worked example

A two-state chain has time-varying rates: $1\to2$ at rate $\lambda(t)$ and $2\to1$ at rate $\mu(t)$ . Compare the quasi-stationary distribution with the actual distribution $p(t)$ .

Find $\nu(t)$ . When is $Q(t)$ weakly irreducible, and when irreducible?

$\nu(t)Q(t)=0$ reads $-\lambda\nu_1+\mu\nu_2=0$ , so $\nu(t)=\big(\mu(t),\lambda(t)\big)/\big(\lambda(t)+\mu(t)\big)$ . This needs $\lambda+\mu>0$ (weak irreducibility); both entries are positive when both rates are positive (irreducibility).
Write the forward equation for $p_1(t)$ alone.

Using $p_2=1-p_1$ : $\ p_1'=-\lambda p_1+\mu(1-p_1)=-(\lambda+\mu)\,p_1+\mu.$
Solve it when the rates are constant.

$p_1(t)=\nu_1+\big(p_1(0)-\nu_1\big)e^{-(\lambda+\mu)t}.$

The gap to $\nu$ closes at rate $\lambda+\mu$ , the absolute value of the nonzero eigenvalue of $Q$ .
Now let the rates vary. What equation does the gap $d(t)=p_1(t)-\nu_1(t)$ satisfy, and what does it say?

Since $\mu=(\lambda+\mu)\nu_1$ , the equation becomes $p_1'=-(\lambda+\mu)(p_1-\nu_1)$ , so

$d'(t)=-\big(\lambda(t)+\mu(t)\big)\,d(t)-\nu_1'(t).$

The chain is pulled toward $\nu(t)$ at rate $\lambda+\mu$ while $\nu(t)$ moves away at speed $\nu_1'$ . If the rates are $\lambda_0(t)/\varepsilon$ and $\mu_0(t)/\varepsilon$ , then $\nu$ does not depend on $\varepsilon$ , the pull wins, and $d$ settles near $-\varepsilon\,\nu_1'/(\lambda_0+\mu_0)$ , which is $O(\varepsilon)$ . That small lag is exactly what Module 4 computes.

Check your understanding

$Q=\begin{pmatrix}-1&1&0\\0&0&0\\0&0&0\end{pmatrix}$ . Is $Q$ weakly irreducible?

Yes: state 1 is transient, so its stationary probability is 0 and $\nu=(0,1,0)$ . $(0,1,0)$ is a solution, but is it the only one? Look at state 3.
No: $(0,s,1-s)$ solves $\nu Q=0,\ \nu\mathbf 1=1$ for every $s\in[0,1]$ .
Yes, because $Q$ has rank $m-1$ .Two rows of $Q$ are zero and the first is nonzero, so the rank is 1, not 2.

States 2 and 3 are both absorbing, two closed classes, so the long-run distribution depends on where the chain starts. Compare the "Two closed classes" generator in the explorer.

In the "Cyclic" example the path of $p(t)$ spirals into $\nu$ . What causes the spiral?

The chain is periodic.A continuous-time chain is never periodic: $p_{ii}(t)\ge e^{-q_it}>0$ for every $t$ .
$Q$ has a pair of complex eigenvalues $-a\pm ib$ , so $e^{Qt}$ contains terms $e^{-at}\cos bt$ and $e^{-at}\sin bt$ .
Rounding error in the numerical solution.The spiral is real. Try smaller and smaller steps and it stays.

The probabilities still converge, at rate

a

, but they overshoot and come back as they go. Generators of reversible chains (detailed balance) have only real eigenvalues, so their paths cannot spiral.

$Q(t)$ is weakly irreducible for every $t$ and varies slowly. Is $p(t)=\nu(t)$ for all $t$ if $p(0)=\nu(0)$ ?

Yes, because $\nu(t)$ is stationary. $\nu(t)$ would be stationary for a chain whose rates stayed frozen at time $t$ . Here they keep moving.
No: $p(t)$ lags behind $\nu(t)$ , unless $\nu(t)$ happens to be constant.

Differentiate

\nu(t)Q(t)=0

: if

p\equiv\nu

then

\nu'=p'=pQ=\nu Q=0

, so

\nu

would have to be constant. The worked example measures the lag.

Practice

The solvability conditionOne direction of the Fredholm alternative

Let $Q$ be a generator. Show that if $fQ=g$ has a solution then $g\mathbf 1=0$ . If $Q$ is weakly irreducible, show the solution is unique once $f\mathbf 1$ is fixed.

Solution

$g\mathbf 1=fQ\mathbf 1=f\cdot0=0$ . If $f_1Q=f_2Q=g$ with $f_1\mathbf 1=f_2\mathbf 1$ , then $h=f_1-f_2$ has $hQ=0$ and $h\mathbf 1=0$ , so $h=0$ by the second form of weak irreducibility.

Why

e^{Qt}

is a stochastic matrixPositivity without probability

Show directly from the matrix exponential that $e^{Qt}$ has nonnegative entries and rows summing to 1 when $Q$ is a generator.

Hint

Choose $c\ge\max_i q_i$ . Then $Q+cI$ has nonnegative entries, and $Q$ and $cI$ commute.

Solution

$e^{Qt}=e^{-ct}e^{(Q+cI)t}=e^{-ct}\sum_n\frac{t^n}{n!}(Q+cI)^n$ , a sum of nonnegative matrices. And $e^{Qt}\mathbf 1=\sum_n\frac{t^n}{n!}Q^n\mathbf 1=\mathbf 1$ since $Q\mathbf 1=0$ . This also gives a way to simulate: run a Poisson clock of rate $c$ and move by $I+Q/c$ at each tick (uniformisation).

SimulationYin–Zhang, Ch. 3: Markovian Models

Birth–death chains and a finite queue

Most Markov models are built by listing what can happen next and at what rate. In a birth–death chain the state only moves up or down by one, and that is enough to solve for the stationary distribution in closed form. Here it is for a queue with a waiting room of size $K$ .

The idea

Birth–death chain. On

\{0,1,\dots,K\}

, the only nonzero off-diagonal rates are

q_{n,n+1}=\lambda_n

(a birth) and

q_{n,n-1}=\mu_n

(a death).

Detailed balance. If all

\lambda_n,\mu_n>0

in range, the stationary distribution satisfies

\nu_n\lambda_n=\nu_{n+1}\mu_{n+1}

, so

\nu_n=\nu_0\prod_{k=0}^{n-1}\frac{\lambda_k}{\mu_{k+1}},\qquad \nu_0 \text{ fixed by } \textstyle\sum_n\nu_n=1.

The reason: the cut between

\{0,\dots,n\}

and

\{n+1,\dots,K\}

is crossed by a single pair of edges, and in equilibrium the flow across a cut is zero.

A finite queue. Customers arrive at rate

\lambda

, one server works at rate

\mu

, and at most

K

customers fit in the system; arrivals to a full system are turned away. This is the birth–death chain with

\lambda_n=\lambda

(

n<K

),

\mu_n=\mu

(

n\ge1

). With

\rho=\lambda/\mu\neq1

,

\nu_n=\frac{(1-\rho)\,\rho^n}{1-\rho^{K+1}},\qquad n=0,\dots,K,

and

\nu_n=1/(K+1)

when

\rho=1

. No condition

\rho<1

is needed: the state space is finite.

Arrivals see time averages. Poisson arrivals are independent of the queue's past, so the fraction of arrivals that find the system full equals the fraction of time it is full,

\nu_K

.

The book's catalogue of models goes on: system reliability (components failing and being repaired), competing risks, Cox processes, random evolutions, seasonal models, and chains inside optimisation algorithms and controlled systems. The common thread is time-scale separation: in each, some rates are naturally much larger than others, and $Q^\varepsilon=\tilde Q/\varepsilon+\hat Q$ makes that explicit.

Explore

Arrival rate λ Service rate μ Capacity K

Speed

Queue length over the last 40 time units.

fraction of time at each length so far

\nu_n

from detailed balance

Press Run. The blue bars converge to the gold ones, and the share of arrivals turned away converges to

\nu_K

, the share of time the queue is full. Try

\lambda>\mu

: the queue does not blow up, it sits near

K

and turns many customers away. Changing a rate restarts the statistics.

Worked example

Two identical machines each fail at rate $\lambda$ while running. One repairer fixes failed machines one at a time at rate $\mu$ . Let $\alpha(t)\in\{0,1,2\}$ be the number of failed machines. Find the long-run availability (at least one machine working) for $\lambda=1,\ \mu=4$ .

Write the generator.

With 0 failed, two machines can fail: rate $2\lambda$ . With 1 failed, one can fail (rate $\lambda$ ) or the repair finishes (rate $\mu$ ). With 2 failed, only repair (rate $\mu$ ):

$Q=\begin{pmatrix}-2\lambda & 2\lambda & 0\\ \mu & -(\lambda+\mu) & \lambda\\ 0 & \mu & -\mu\end{pmatrix}.$
Use detailed balance to express $\nu_1,\nu_2$ in terms of $\nu_0$ .

$\nu_1=\frac{2\lambda}{\mu}\nu_0$ and $\nu_2=\frac{\lambda}{\mu}\nu_1=\frac{2\lambda^2}{\mu^2}\nu_0$ .
Normalise with $\lambda=1,\ \mu=4$ , and give the availability.

$\nu\propto(1,\ \tfrac12,\ \tfrac18)$ , and the entries sum to $\tfrac{13}8$ , so $\nu=\tfrac1{13}(8,4,1)$ . Availability $=1-\nu_2=\tfrac{12}{13}\approx0.923$ .
Suppose repairs are much faster than failures: $\mu=\mu_0/\varepsilon$ . How does the downtime scale?

$\nu_2=\dfrac{2\lambda^2/\mu^2}{1+2\lambda/\mu+2\lambda^2/\mu^2}=2\lambda^2\varepsilon^2/\mu_0^2+O(\varepsilon^3)$ . A machine is under repair a fraction $O(\varepsilon)$ of the time, and the second one must fail during that short repair, which has probability $O(\varepsilon)$ . Fast repairs make the chain spend almost all its time in state 0: a first glimpse of a two-time-scale model.

Check your understanding

In the finite queue, take $\lambda=\mu$ . What is the stationary distribution?

There is none: the queue is critically loaded.That is true for an infinite waiting room. With room for only $K$ customers the chain is finite and irreducible.
Uniform: $\nu_n=1/(K+1)$ .
Concentrated at $n=K$ .Detailed balance with $\lambda=\mu$ gives $\nu_{n+1}=\nu_n$ : no state is preferred.

\nu_{n+1}/\nu_n=\lambda/\mu=1

. Set

\lambda=\mu

in the explorer: all the gold bars have the same height.

Why does the fraction of arrivals that are turned away equal $\nu_K$ , the fraction of time the system is full?

Because the server is idle exactly $\nu_0$ of the time.True, but that is about state 0. What is special about the moments when customers arrive?
Arrivals come from a Poisson process independent of the past of the queue, so the state an arrival sees has the time-stationary distribution.
It doesn't: arrivals are more likely when the queue is short.Here the arrival rate is $\lambda$ in every state, so arrivals do not prefer any state.

This is "Poisson arrivals see time averages". It fails when the arrival rate depends on the state: then arrivals sample states in proportion to

\nu_n\lambda_n

, not

\nu_n

.

A birth–death chain on $\{0,\dots,K\}$ has every $\lambda_n$ and $\mu_n$ positive. Is it reversible (does detailed balance hold)?

Always.
Only when $\lambda_n=\mu_{n+1}$ for all $n$ .That is the special case of a uniform $\nu$ . Detailed balance only needs $\nu_n\lambda_n=\nu_{n+1}\mu_{n+1}$ , which can always be arranged.
Only for the finite queue.The cut argument uses nothing about the particular rates.

Since the chain lives on a path graph, each cut is crossed only by

n\to n+1

and

n+1\to n

, so zero net flow across every cut is detailed balance. Chains with cycles, like the cyclic generator of Module 2, need not be reversible.

Practice

A tiny queue by handLittle's law

Take $K=2$ , $\lambda=1$ , $\mu=2$ . Find $\nu$ , the mean number in the system $L$ , the throughput, and the mean time a customer who gets in spends in the system.

Solution

$\rho=\tfrac12$ , so $\nu\propto(1,\tfrac12,\tfrac14)$ , giving $\nu=\tfrac17(4,2,1)$ . $L=\tfrac17(0+2+2)=\tfrac47$ . Accepted customers arrive at rate $\lambda(1-\nu_2)=\tfrac67$ , which is also the throughput. Little's law: $W=L/\tfrac67=\tfrac23$ .

Zero flow across a cutDeriving detailed balance

For a birth–death chain with stationary $\nu$ , add up the equations $(\nu Q)_k=0$ for $k=0,\dots,n$ and conclude $\nu_n\lambda_n=\nu_{n+1}\mu_{n+1}$ .

Solution

$\sum_{k\le n}(\nu Q)_k=\sum_i\nu_i\sum_{k\le n}q_{ik}$ . For $i\le n$ , the inner sum is minus the rate of leaving $\{0,\dots,n\}$ from $i$ , which is $-\lambda_n$ for $i=n$ and 0 otherwise. For $i>n$ it is the rate of entering $\{0,\dots,n\}$ , which is $\mu_{n+1}$ for $i=n+1$ and 0 otherwise. So the sum is $-\nu_n\lambda_n+\nu_{n+1}\mu_{n+1}=0$ .

InteractiveYin–Zhang, Ch. 4: Asymptotic Expansions of Solutions for Forward Equations

A fast chain: outer expansion and initial layer

Speed a chain up by a factor $1/\varepsilon$ while its rates drift in ordinary time. Its distribution snaps to the quasi-stationary distribution within a time of order $\varepsilon$ , then follows it. A matched asymptotic expansion describes both phases, with an error you can watch shrink in proportion to $\varepsilon$ .

The idea

Let $\alpha^\varepsilon(t)$ be generated by $Q(t)/\varepsilon$ , with $Q(t)$ weakly irreducible and smooth on $[0,T]$ . Its distribution solves

\frac{dp^\varepsilon(t)}{dt}=\frac1\varepsilon\,p^\varepsilon(t)\,Q(t),\qquad p^\varepsilon(0)=p^0.

Two-part expansion. Look for

p^\varepsilon(t)=\sum_{i=0}^n\varepsilon^i\varphi_i(t)+\sum_{i=0}^n\varepsilon^i\psi_i(t/\varepsilon)+e^\varepsilon_n(t)

. The outer terms

\varphi_i

describe

t

away from 0; the initial-layer terms

\psi_i(\tau)

, in the stretched time

\tau=t/\varepsilon

, repair the start and decay exponentially in

\tau

.

Outer terms. Substitute and match powers of $\varepsilon$ :

\varphi_0(t)Q(t)=0,\qquad \varphi_k(t)Q(t)=\frac{d\varphi_{k-1}(t)}{dt}\quad(k\ge1),

with $\varphi_0\mathbf 1=1$ and $\varphi_k\mathbf 1=0$ for $k\ge1$ . So $\varphi_0(t)=\nu(t)$ , the quasi-stationary distribution, and each later equation is solvable by the Fredholm alternative because $\varphi_{k-1}'\mathbf 1=0$ .

Initial layer. In the stretched time, $\frac{d\psi_0}{d\tau}=\psi_0(\tau)Q(0)$ and $\psi_0(0)=p^0-\nu(0)$ , so

\psi_0(\tau)=\big(p^0-\nu(0)\big)\,e^{Q(0)\tau},

which decays like $e^{-\kappa\tau}$ because $(p^0-\nu(0))\mathbf 1=0$ .

Error bound. Under these assumptions (with

Q

smooth enough for the order

n

),

\ \sup_{0\le t\le T}\Big|p^\varepsilon(t)-\sum_{i\le n}\varepsilon^i\varphi_i(t)-\sum_{i\le n}\varepsilon^i\psi_i(t/\varepsilon)\Big|\le K\varepsilon^{n+1}.

For

n=0

:

p^\varepsilon(t)=\nu(t)+\psi_0(t/\varepsilon)+O(\varepsilon)

, and

p^\varepsilon(t)\to\nu(t)

for every

t>0

.

Watch out. The convergence is not uniform near

t=0

unless

p^0=\nu(0)

, and even then the higher initial-layer terms

\psi_1,\psi_2,\dots

need not vanish.

Explore

Time-scale parameter ε Starting probability of state 1, p₁(0)

Two states. Rates $\lambda(t)/\varepsilon$ for $1\to2$ and $\mu/\varepsilon$ for $2\to1$ , with $\lambda(t)=1.5+\sin\pi t$ and $\mu=1$ , so $\nu_1(t)=1/(\lambda(t)+1)$ .

p^\varepsilon_1(t)

, exact

\nu_1(t)

\nu_1(t)+\psi_0(t/\varepsilon)

, on

0\le t\le2

The initial layer in stretched time

\tau=t/\varepsilon

. Here

\psi_0(\tau)=(p_1(0)-0.4)\,e^{-2.5\tau}

for state 1.

Worked example

For the two-state chain of the explorer, with rates $\lambda(t)/\varepsilon$ and $\mu(t)/\varepsilon$ , find $\varphi_0$ , $\psi_0$ and the first outer correction $\varphi_1$ , and evaluate them at $t=1$ .

What is $\varphi_0(t)$ ?

The quasi-stationary distribution: $\varphi_0(t)=\nu(t)=\big(\mu,\lambda\big)/(\lambda+\mu)$ . With $\mu=1$ , $\nu_1(t)=1/(\lambda(t)+1)$ .
Find $\psi_0(\tau)$ .

A zero-sum row vector $(x,-x)$ times $Q(0)$ is $-(\lambda(0)+\mu(0))(x,-x)$ , so

$\psi_0(\tau)=\big(p^0_1-\nu_1(0)\big)\,(1,-1)\,e^{-(\lambda(0)+\mu(0))\tau}.$

In the explorer $\lambda(0)+\mu=2.5$ and $\nu_1(0)=0.4$ .
Solve $\varphi_1Q=\varphi_0'$ with $\varphi_1\mathbf 1=0$ .

Write $\varphi_1=(c,-c)$ . Then $\varphi_1Q=-c(\lambda+\mu)(1,-1)$ and $\varphi_0'=\nu_1'(1,-1)$ , so

$\varphi_1(t)=-\frac{\nu_1'(t)}{\lambda(t)+\mu(t)}\,(1,-1).$

When $\nu_1$ rises, $p^\varepsilon_1$ sits slightly below it: the chain lags behind its moving target by about $\varepsilon\nu_1'/(\lambda+\mu)$ , as in Module 2.
Evaluate at $t=1$ , where $\lambda=1.5$ and $\lambda'=\pi\cos\pi=-\pi$ .

$\nu_1(1)=1/2.5=0.4$ and $\nu_1'=-\lambda'/(\lambda+1)^2=\pi/6.25\approx0.5027$ , so $\varphi_{1,1}(1)\approx-0.5027/2.5\approx-0.2011$ . For $\varepsilon=0.05$ : $p^\varepsilon_1(1)\approx0.4-0.0101=0.3899$ . The explorer's last two boxes compare this with the exact value.

Check your understanding

Why can't the outer expansion $\sum\varepsilon^i\varphi_i(t)$ alone match the initial condition $p^\varepsilon(0)=p^0$ ?

$\varphi_0(0)=\nu(0)$ is forced by $\varphi_0Q=0$ ; it is not free to equal $p^0$ .
It can, if enough terms are kept.Every $\varphi_k$ with $k\ge1$ carries a factor $\varepsilon^k$ . No combination of them can close an $O(1)$ gap $p^0-\nu(0)$ .
Because the outer terms are not smooth at $t=0$ .They are as smooth as $Q$ . The trouble is their value at 0, not their smoothness.

The outer equations are algebraic at leading order, so they lose the initial condition. The initial-layer term, which lives on the fast time

\tau=t/\varepsilon

, carries the mismatch and dies out within a few multiples of

\varepsilon

.

$Q(0)$ is a 3-state generator with nonzero eigenvalues $-2$ and $-5$ . For small $\varepsilon$ , how does the initial layer $\psi_0(t/\varepsilon)$ decay?

Like $e^{-2t}$ .The layer lives on the fast time $t/\varepsilon$ .
Like $e^{-5t/\varepsilon}$ .The fast mode dies first. What is left after that?
Like $e^{-2t/\varepsilon}$ (for a typical $p^0$ ).

\psi_0(\tau)=(p^0-\nu(0))e^{Q(0)\tau}

is a combination of

e^{-2\tau}

and

e^{-5\tau}

, and the slower mode dominates. With

\varepsilon=0.01

the layer is gone, to a few percent, by

t\approx0.02

.

Suppose $p^0=\nu(0)$ . Which is true?

$p^\varepsilon(t)=\nu(t)$ exactly for all $t$ .Only if $\nu(t)$ is constant. Otherwise the chain lags behind its moving target.
$\psi_0\equiv0$ , but $p^\varepsilon(t)-\nu(t)$ is still of order $\varepsilon$ .
The error becomes $O(\varepsilon^2)$ .The $\varepsilon\varphi_1$ term is still there; only the zeroth-order layer disappeared.

With

p^0=\nu(0)

the zeroth-order layer vanishes, but

\varepsilon\varphi_1(t)

does not, and

\varphi_1(0)\neq0

when

\nu'(0)\neq0

, so

\psi_1(0)=-\varphi_1(0)

is nonzero too. In the explorer, set

p_1(0)=0.40

: the gold curve coincides with the green one, but the blue one still lags by

O(\varepsilon)

.

Practice

Constant generatorsWhen the expansion terminates

If $Q$ does not depend on $t$ , show that $p^\varepsilon(t)=\nu+(p^0-\nu)e^{Qt/\varepsilon}$ exactly. Which terms of the expansion are nonzero?

Solution

$p^\varepsilon(t)=p^0e^{Qt/\varepsilon}=\nu e^{Qt/\varepsilon}+(p^0-\nu)e^{Qt/\varepsilon}$ , and $\nu e^{Qs}=\nu$ since $\nu Q=0$ . So $\varphi_0=\nu$ , $\psi_0(\tau)=(p^0-\nu)e^{Q\tau}$ , and all other terms vanish: $\varphi_1Q=\nu'=0$ with $\varphi_1\mathbf 1=0$ forces $\varphi_1=0$ , and so on.

Uniform away from zeroReading the error bound

Use the $n=0$ bound to show that $\sup_{\delta\le t\le T}|p^\varepsilon(t)-\nu(t)|=O(\varepsilon)$ for each fixed $\delta>0$ , while $\sup_{0\le t\le T}|p^\varepsilon(t)-\nu(t)|$ does not tend to 0 unless $p^0=\nu(0)$ .

Solution

$|p^\varepsilon-\nu|\le|\psi_0(t/\varepsilon)|+K\varepsilon\le K'e^{-\kappa\delta/\varepsilon}+K\varepsilon$ on $[\delta,T]$ , and $e^{-\kappa\delta/\varepsilon}=o(\varepsilon)$ . At $t=0$ the difference is $|p^0-\nu(0)|$ for every $\varepsilon$ .

InteractiveYin–Zhang, Ch. 4: Asymptotic Expansions of Solutions for Forward Equations

Weak and strong interactions

Now the fast part of the generator splits the states into groups that do not communicate, and a slow part moves the chain between groups. Inside each group the distribution equilibrates almost at once. The weights of the groups follow a smaller Markov chain whose generator is an average of the slow rates.

The idea

Setup.

Q^\varepsilon(t)=\tilde Q(t)/\varepsilon+\hat Q(t)

with

\tilde Q=\mathrm{diag}(\tilde Q^1,\dots,\tilde Q^l)

. Block

\tilde Q^k

acts on the group

\mathcal M_k=\{s_{k1},\dots,s_{km_k}\}

, is weakly irreducible, and has quasi-stationary distribution

\nu^k(t)

.

\hat Q(t)

is any generator on all

m=m_1+\dots+m_l

states.

Aggregated generator. With

\tilde{\mathbf 1}=\mathrm{diag}(\mathbf 1_{m_1},\dots,\mathbf 1_{m_l})

, an

m\times l

matrix,

\bar Q(t)=\mathrm{diag}\big(\nu^1(t),\dots,\nu^l(t)\big)\,\hat Q(t)\,\tilde{\mathbf 1}.

Entry

\bar q_{kj}

is the slow rate from group

k

into group

j

, averaged over the states of group

k

with weights

\nu^k

.

Leading-order behaviour. Uniformly on

[0,T]

,

p^\varepsilon(t)=\big(\vartheta_1(t)\nu^1(t),\dots,\vartheta_l(t)\nu^l(t)\big)+\psi_0(t/\varepsilon)+O(\varepsilon),

where

\vartheta=(\vartheta_1,\dots,\vartheta_l)

solves

\vartheta'(t)=\vartheta(t)\bar Q(t)

with

\vartheta(0)=p^0\tilde{\mathbf 1}

, the vector of initial group weights. The initial layer

\psi_0(\tau)=\big(p^0-\varphi_0(0)\big)e^{\tilde Q(0)\tau}

, where

\varphi_0(0)=(\vartheta_1(0)\nu^1(0),\dots,\vartheta_l(0)\nu^l(0))

, decays exponentially.

So $P(\alpha^\varepsilon(t)\in\mathcal M_k)\to\vartheta_k(t)$ , and given the group, the chain is spread according to $\nu^k(t)$ . The group weights have no zeroth-order layer: $e^{\tilde Q\tau}$ only moves probability inside groups, so $\psi_0(\tau)\tilde{\mathbf 1}=0$ .

What survives the limit.

\hat Q

enters only through

\bar Q

. Slow rates between two states of the same group change nothing at leading order, and two slow generators with the same

\bar Q

give the same

\vartheta

.

Explore

Time-scale parameter ε

Start

Dashed curves include ψ₀

Fast rates: $s_{11}\to s_{12}$ 2, $s_{12}\to s_{11}$ 1, $s_{21}\to s_{22}$ 1, $s_{22}\to s_{21}$ 3 (all divided by ε). Slow rates: $s_{11}\to s_{21}$ 0.9, $s_{12}\to s_{22}$ 0.3, $s_{21}\to s_{12}$ 0.4, $s_{22}\to s_{11}$ 1.2.

Solid:

P(\alpha^\varepsilon(t)=s_{kj})

. Dashed:

\vartheta_k(t)\nu^k_j

, with or without the initial layer.

P(\alpha^\varepsilon(t)\in\mathcal M_1)

P(\alpha^\varepsilon(t)\in\mathcal M_2)

; dashed:

\vartheta(t)

from

\vartheta'=\vartheta\bar Q

Worked example

Compute everything the explorer uses: $\nu^1,\nu^2$ , $\bar Q$ , and $\vartheta(t)$ when the chain starts in $s_{11}$ . Then find the long-run distribution of the full chain to leading order.

Find $\nu^1$ and $\nu^2$ .

For a two-state block with rates $a$ (first to second) and $b$ (back), $\nu=(b,a)/(a+b)$ . So $\nu^1=(1,2)/3$ and $\nu^2=(3,1)/4$ .
Compute $\bar q_{12}$ and $\bar q_{21}$ .

$\bar q_{12}=\nu^1_1\cdot0.9+\nu^1_2\cdot0.3=0.3+0.2=0.5$ and $\bar q_{21}=\nu^2_1\cdot0.4+\nu^2_2\cdot1.2=0.3+0.3=0.6$ . So

$\bar Q=\begin{pmatrix}-0.5 & 0.5\\ 0.6 & -0.6\end{pmatrix}.$
Solve for $\vartheta(t)$ with $\vartheta(0)=(1,0)$ .

By the two-state formula from Module 1, $\vartheta_1(t)=\tfrac6{11}+\tfrac5{11}e^{-1.1t}$ and $\vartheta_2=1-\vartheta_1$ .
What is the stationary distribution of the full chain, to leading order in $\varepsilon$ ?

Let $t\to\infty$ in the leading term: $\vartheta\to(\tfrac6{11},\tfrac5{11})$ , so

$\nu^\varepsilon\approx\big(\tfrac6{11}\cdot\tfrac13,\ \tfrac6{11}\cdot\tfrac23,\ \tfrac5{11}\cdot\tfrac34,\ \tfrac5{11}\cdot\tfrac14\big)=\big(\tfrac2{11},\tfrac4{11},\tfrac{15}{44},\tfrac5{44}\big).$

A direct computation with $\varepsilon=0.01$ gives $(0.1817,0.3637,0.3410,0.1136)$ , against $(0.1818,0.3636,0.3409,0.1136)$ here. The 4-state problem has been reduced to three 2-state problems: one for each group, and one for the aggregated chain.

Check your understanding

Group 1 has $\nu^1=(\tfrac12,\tfrac12)$ . From $s_{11}$ the slow rate into group 2 is 2; from $s_{12}$ it is 0. What is $\bar q_{12}$ ?

2That is the rate from $s_{11}$ only. The chain spends half of its group-1 time in $s_{12}$ , where it cannot leave.
1
0, the smaller rateThe aggregated rate is an average, not a minimum.

\bar q_{12}=\tfrac12\cdot2+\tfrac12\cdot0=1

. Because the fast chain visits both states many times between slow jumps, the slow clock sees the time-averaged exit rate.

In the right-hand plot the group probabilities show no fast transient at the start, even when the chain starts in a single state. Why?

$\tilde Q$ is block diagonal, so the fast dynamics only move probability inside a group; each group's total changes at rate $O(1)$ .
Because $\nu^k$ sums to 1.That is true but says nothing about time evolution. What does $e^{\tilde Q\tau}$ do to the total mass of a group?
They do have one, it is just too small to see.At leading order it is exactly zero: $\psi_0\tilde{\mathbf 1}=(p^0-\varphi_0(0))e^{\tilde Q\tau}\tilde{\mathbf 1}=(p^0-\varphi_0(0))\tilde{\mathbf 1}=0$ .

The initial layer reshuffles mass within each group toward

\nu^k

; the group totals

\vartheta_k

start at their correct values

p^0\tilde{\mathbf 1}

and evolve slowly. (An

O(\varepsilon)

layer can remain.)

What happens if $\hat Q=0$ ?

$p^\varepsilon(t)$ converges to a unique stationary distribution.With $\hat Q=0$ the groups never communicate, so the chain is not irreducible.
$\vartheta$ stays at its initial value: the group weights are frozen, and inside each group the distribution tends to $\nu^k$ .
The theorem fails.It still applies, with $\bar Q=0$ .

\bar Q=0

gives

\vartheta(t)=\vartheta(0)

. The weak interaction

\hat Q

is what lets probability move between groups, on the slow time scale.

Practice

\bar Q

is a generatorTwo lines

Show that $\bar Q=\mathrm{diag}(\nu^1,\dots,\nu^l)\hat Q\tilde{\mathbf 1}$ has nonnegative off-diagonal entries and zero row sums.

Solution

For $k\neq j$ , $\bar q_{kj}=\sum_i\nu^k_i\sum_{s\in\mathcal M_j}\hat q_{s_{ki},s}\ge0$ , a nonnegative combination of off-diagonal entries of $\hat Q$ . Row sums: $\bar Q\mathbf 1_l=\mathrm{diag}(\nu^k)\hat Q\tilde{\mathbf 1}\mathbf 1_l=\mathrm{diag}(\nu^k)\hat Q\mathbf 1_m=0$ .

Where

\bar Q

comes fromAveraging the forward equation

Multiply the forward equation $p^{\varepsilon\prime}=p^\varepsilon(\tilde Q/\varepsilon+\hat Q)$ on the right by $\tilde{\mathbf 1}$ . Assuming $p^\varepsilon\approx(\vartheta_1\nu^1,\dots,\vartheta_l\nu^l)$ , derive $\vartheta'=\vartheta\bar Q$ .

Solution

Since $\tilde Q\tilde{\mathbf 1}=0$ , the $1/\varepsilon$ term drops out: $(p^\varepsilon\tilde{\mathbf 1})'=p^\varepsilon\hat Q\tilde{\mathbf 1}$ . The left side is the vector of group weights, approximately $\vartheta$ . On the right, $(\vartheta_1\nu^1,\dots,\vartheta_l\nu^l)=\vartheta\,\mathrm{diag}(\nu^1,\dots,\nu^l)$ , so the right side is $\vartheta\,\mathrm{diag}(\nu^k)\hat Q\tilde{\mathbf 1}=\vartheta\bar Q$ .

InteractiveYin–Zhang, Ch. 4: Asymptotic Expansions of Solutions for Forward Equations

Transient and absorbing states

Add states that the fast dynamics leaves and never returns to. Their probability dies out within a time of order $\varepsilon$ , and the mass they held is handed to the recurrent groups in proportions that one linear solve gives. Absorbing states are the simplest recurrent groups of all.

The idea

Setup. After relabelling, the fast generator has the form

\tilde Q=\begin{pmatrix}\tilde Q^1 & & & \\ & \ddots & & \\ & & \tilde Q^l & \\ \tilde Q^1_* & \cdots & \tilde Q^l_* & \tilde Q_*\end{pmatrix},

where each

\tilde Q^k

is a weakly irreducible generator on a recurrent group

\mathcal M_k

, the last block row belongs to the transient states

\mathcal M_*

, and

\tilde Q_*

is Hurwitz: all its eigenvalues have negative real part. A group with one state and

\tilde Q^k=(0)

is an absorbing state, with

\nu^k=1

.

Exit probabilities.

a_k=-\tilde Q_*^{-1}\tilde Q^k_*\mathbf 1_{m_k}\in\R^{m_*}

. Its entry for a transient state

s

is the probability that the fast chain started at

s

first enters

\mathcal M_k

. Since

\tilde Q\mathbf 1=0

,

\sum_ka_k=\mathbf 1_{m_*}

.

Leading order.

p^\varepsilon(t)=\big(\vartheta_1(t)\nu^1,\dots,\vartheta_l(t)\nu^l,\ 0_{m_*}\big)+\psi_0(t/\varepsilon)+O(\varepsilon)

, where

\vartheta_k(0)=p^{0,k}\mathbf 1_{m_k}+p^{0,*}a_k,\qquad \vartheta'=\vartheta\bar Q_*,\qquad \bar Q_*=\mathrm{diag}(\nu^1,\dots,\nu^l)\big(\hat Q^{11}\tilde{\mathbf 1}+\hat Q^{12}(a_1,\dots,a_l)\big).

Here

\hat Q^{11}

holds the slow rates among recurrent states and

\hat Q^{12}

the slow rates from recurrent into transient states.

Read $\bar Q_*$ as a routing rule. A slow jump straight into another group counts as it is. A slow jump into a transient state is followed, within time $O(\varepsilon)$ , by a fast exit into group $k$ with probability given by $a_k$ ; only the part that lands in a different group is a real change of group.

Explore

Time-scale parameter ε

Start

Recurrent group $\mathcal M_1=\{s_{11},s_{12}\}$ (fast rates 1 and 2 between them), absorbing state $\mathcal M_2=\{s_{21}\}$ , transient states $s^*_1,s^*_2$ . Fast: $s^*_1\to s_{11}$ 2, $s^*_1\to s^*_2$ 1, $s^*_2\to s^*_1$ 1, $s^*_2\to s_{12}$ 0.5, $s^*_2\to s_{21}$ 1. Slow: $s_{11}\to s_{21}$ 0.3, $s_{12}\to s^*_1$ 0.4, $s_{21}\to s^*_2$ 0.5.

P(\mathcal M_1)

P(\mathcal M_2)

P(\mathcal M_*)

; dashed:

\vartheta(t)

The start in stretched time

\tau=t/\varepsilon

. Dashed: the predicted split

\vartheta(0)

.

Worked example

For the chain in the explorer, compute the exit probabilities $a_1,a_2$ and the aggregated generator $\bar Q_*$ .

Write $\tilde Q_*$ and check that it is Hurwitz.

On $(s^*_1,s^*_2)$ : $\tilde Q_*=\begin{pmatrix}-3&1\\1&-2.5\end{pmatrix}$ . Its row sums are $-2$ and $-1.5$ , so each diagonal entry outweighs the rest of its row and, by Gershgorin's theorem, every eigenvalue has negative real part. (They are $\approx-1.72$ and $-3.78$ .)
Find $h(s)=P(\text{enter }\mathcal M_1\text{ first}\mid s)$ by a first-step argument.

From $s^*_1$ : with probability $\tfrac23$ go straight to $s_{11}$ , with $\tfrac13$ to $s^*_2$ . From $s^*_2$ : $\tfrac{0.5}{2.5}=\tfrac15$ to $s_{12}$ , $\tfrac25$ to $s_{21}$ , $\tfrac25$ to $s^*_1$ . So $h_1=\tfrac23+\tfrac13h_2$ and $h_2=\tfrac15+\tfrac25h_1$ , giving $h_1=\tfrac{11}{13}$ , $h_2=\tfrac7{13}$ . This is $a_1=-\tilde Q_*^{-1}\tilde Q^1_*\mathbf 1$ written out; and $a_2=\mathbf 1-a_1=(\tfrac2{13},\tfrac6{13})$ .
Starting in $s^*_1$ , what is $\vartheta(0)$ ?

$\vartheta(0)=(\tfrac{11}{13},\tfrac2{13})\approx(0.846,0.154)$ . The explorer's right panel shows the mass settling to these levels within a few multiples of $\varepsilon$ .
Compute $\bar q_{12}$ and $\bar q_{21}$ .

$\nu^1=(\tfrac23,\tfrac13)$ . From $\mathcal M_1$ : $s_{11}\to s_{21}$ is a direct change of group, weight $\tfrac23\cdot0.3=0.2$ ; $s_{12}\to s^*_1$ changes group only if the chain then exits to $\mathcal M_2$ , weight $\tfrac13\cdot0.4\cdot\tfrac2{13}=\tfrac{0.8}{39}$ . So $\bar q_{12}=\tfrac{43}{195}\approx0.2205$ . From $\mathcal M_2$ : $s_{21}\to s^*_2$ , then on to $\mathcal M_1$ with probability $\tfrac7{13}$ : $\bar q_{21}=0.5\cdot\tfrac7{13}=\tfrac7{26}\approx0.2692$ .

Check your understanding

A transient state $s^*$ has fast rate 3 into group 1 and fast rate 1 into group 2, and no other moves. The chain starts at $s^*$ . What is $\vartheta(0)$ ?

$(\tfrac12,\tfrac12)$ The two exits are not equally likely: they race at rates 3 and 1.
$(\tfrac34,\tfrac14)$
$(3,1)$ These are rates; $\vartheta(0)$ is a probability vector.

\tilde Q_*=(-4)

, so

a_1=-(-4)^{-1}\cdot3=\tfrac34

and

a_2=\tfrac14

: the first exit clock to ring decides.

Transient states have probability $O(\varepsilon)$ after the initial layer. Why not simply delete them and use $\mathrm{diag}(\nu^k)\hat Q^{11}\tilde{\mathbf 1}$ as the aggregated generator?

That would be fine: states with vanishing probability cannot matter.Small probability is not the same as small flow. Slow jumps into transient states happen at rate $O(1)$ .
Slow jumps into transient states happen at rate $O(1)$ , and a fraction of them end in a different group; dropping them loses real group changes.
Because $\tilde Q_*$ would no longer be invertible.Deleting the states removes $\tilde Q_*$ altogether. The problem is what happens to $\bar Q$ .

In the explorer, all slow movement from

\mathcal M_2

back to

\mathcal M_1

goes through

s^*_2

. Delete the transient states and

\bar q_{21}

would be 0 instead of

7/26

.

For $t\gg\varepsilon$ , the total probability of the transient states is:

exactly 0.The slow part keeps feeding the transient states.
of order $\varepsilon$ : fed at rate $O(1)$ by slow jumps, drained at rate $O(1/\varepsilon)$ .
decaying like $e^{-t}$ .The drain is fast, on the scale $t/\varepsilon$ , but the feed never stops.

A balance of flows: inflow

\approx c

per unit time, mean stay

\approx\varepsilon/\kappa

, so the mass is

\approx c\varepsilon/\kappa

. The explorer's "max P(transient) ÷ ε" stays roughly constant as you change ε.

Practice

Two absorbing statesA slow leak back

States 1 and 2 are absorbing for the fast part; state 3 is transient with fast rates 1 into state 1 and 2 into state 2. Slowly, state 1 jumps to 3 at rate $b_1$ and state 2 jumps to 3 at rate $b_2$ . Find $\bar Q_*$ and the long-run fraction of time in state 1 as $\varepsilon\to0$ .

Solution

$a=(\tfrac13,\tfrac23)$ . From state 1, a slow jump leads to state 2 with probability $\tfrac23$ : $\bar q_{12}=\tfrac23b_1$ . Similarly $\bar q_{21}=\tfrac13b_2$ . The stationary weight of state 1 is $\bar q_{21}/(\bar q_{12}+\bar q_{21})=b_2/(2b_1+b_2)$ .

The exit probabilities add upUsing the row sums

Show that $\sum_ka_k=\mathbf 1_{m_*}$ .

Solution

The transient rows of $\tilde Q\mathbf 1=0$ say $\sum_k\tilde Q^k_*\mathbf 1_{m_k}+\tilde Q_*\mathbf 1_{m_*}=0$ . Multiply by $-\tilde Q_*^{-1}$ : $\sum_ka_k=-\tilde Q_*^{-1}\big(-\tilde Q_*\mathbf 1_{m_*}\big)=\mathbf 1_{m_*}$ . Starting from a transient state, the fast chain leaves $\mathcal M_*$ with probability 1.

SimulationYin–Zhang, Ch. 5: Occupation Measures: Asymptotic Properties and Ramification

Occupation measures and their fluctuations

The forward equation tells us $P(\alpha^\varepsilon(t)=i)$ . Applications usually need something about one path: how much time it actually spent in state $i$ . That time, minus its expected share, goes to zero like $\sqrt\varepsilon$ , and divided by $\sqrt\varepsilon$ it becomes Gaussian.

The idea

Occupation measure. For

\alpha^\varepsilon

generated by

Q(t)/\varepsilon

(weakly irreducible), set

O^\varepsilon_i(t)=\int_0^t\big(I_{\{\alpha^\varepsilon(s)=i\}}-\nu_i(s)\big)\,ds,

or more generally with a bounded weight

\beta_i(s)

inside the integral, as needed when a control multiplies the indicator.

Size.

\sup_{0\le t\le T}\mathbb E|O^\varepsilon(t)|^2=O(\varepsilon)

. There is also an exponential bound: for a small enough constant

\theta

(depending on

T

),

\mathbb E\exp\big(\theta\sup_{t\le T}|O^\varepsilon(t)|/\sqrt\varepsilon\big)

stays bounded as

\varepsilon\to0

. So large values of

O^\varepsilon/\sqrt\varepsilon

are exponentially rare.

Central limit theorem.

n^\varepsilon(t)=O^\varepsilon(t)/\sqrt\varepsilon

converges weakly to a Gaussian process with independent increments, mean 0 and covariance

\int_0^tA(s)\,ds

, where

A_{ij}(t)=\nu_i(t)\int_0^\infty q^0_{ij}(r,t)\,dr+\nu_j(t)\int_0^\infty q^0_{ji}(r,t)\,dr,\qquad Q^0(r,t)=e^{Q(t)r}-\mathbf 1\nu(t).

For a constant $Q$ this is $A_{ii}=2\nu_iZ_{ii}$ , with the deviation matrix $Z=\int_0^\infty(e^{Qr}-\mathbf 1\nu)\,dr=(\mathbf 1\nu-Q)^{-1}-\mathbf 1\nu$ .

Why $\sqrt\varepsilon$ . In time $t$ the fast chain goes through about $t/\varepsilon$ nearly independent cycles, each contributing a centered amount of size $\varepsilon$ to $O^\varepsilon$ . A sum of $t/\varepsilon$ such terms has standard deviation $\varepsilon\sqrt{t/\varepsilon}=\sqrt{\varepsilon t}$ .

Explore

State i

ε

Generator $Q/\varepsilon$ with $Q=\begin{pmatrix}-2&1&1\\2&-3&1\\1&2&-3\end{pmatrix}$ , so $\nu=(\tfrac7{16},\tfrac5{16},\tfrac14)$ . Time runs over $[0,1]$ .

24 paths of

O^\varepsilon_i(t)

. Note the vertical scale: it is set by

\sqrt\varepsilon

.

400 values of

n^\varepsilon(1)

normal density, variance

A

Worked example

A two-state chain has rates $\lambda/\varepsilon$ ( $1\to2$ ) and $\mu/\varepsilon$ ( $2\to1$ ). Find the limiting variance of $n^\varepsilon_1(t)$ , and use it to estimate the spread of the time spent in state 1 during $[0,1]$ when $\lambda=\mu=1$ and $\varepsilon=0.01$ .

What is $p_{11}(r)-\nu_1$ for the unscaled generator?

From Module 1, $p_{11}(r)=\nu_1+\frac{\lambda}{\lambda+\mu}e^{-(\lambda+\mu)r}$ , with $\nu_1=\frac\mu{\lambda+\mu}$ .
Compute $Z_{11}$ .

$Z_{11}=\int_0^\infty\frac{\lambda}{\lambda+\mu}e^{-(\lambda+\mu)r}\,dr=\frac{\lambda}{(\lambda+\mu)^2}.$
Find $A_{11}$ .

$A_{11}=2\nu_1Z_{11}=\frac{2\lambda\mu}{(\lambda+\mu)^3}.$

Check against renewal theory: on and off periods with means $1/\lambda$ , $1/\mu$ and variances $1/\lambda^2$ , $1/\mu^2$ give variance per unit time $\frac{m_{\text{off}}^2\sigma_{\text{on}}^2+m_{\text{on}}^2\sigma_{\text{off}}^2}{(m_{\text{on}}+m_{\text{off}})^3}$ , which is the same expression.
Now $\lambda=\mu=1$ , $\varepsilon=0.01$ . How spread out is the time spent in state 1 during $[0,1]$ ?

$A=\tfrac14$ , so $\mathrm{Var}\,O^\varepsilon_1(1)\approx\varepsilon A=0.0025$ , a standard deviation of 0.05. The time in state 1 is $0.5\pm0.05$ or so, even though the chain makes about 100 jumps. Fluctuations of order $\sqrt\varepsilon$ are much larger than the $O(\varepsilon)$ errors of the forward equation.

Check your understanding

You divide $\varepsilon$ by 4. The typical size of $O^\varepsilon_i(1)$ is:

divided by 4.That is how the mean $\mathbb EO^\varepsilon$ , of order $\varepsilon$ , scales. The fluctuations are larger.
divided by 2.
unchanged.That is true of $n^\varepsilon=O^\varepsilon/\sqrt\varepsilon$ , not of $O^\varepsilon$ .

The standard deviation is

\approx\sqrt{\varepsilon A}

. Step ε from 0.01 to 0.003 or so in the explorer and watch the vertical scale of the left panel.

Which statement about $\mathbb EO^\varepsilon_i(t)$ and the typical value of $|O^\varepsilon_i(t)|$ is right?

Both are $O(\varepsilon)$ .The mean is. A random quantity can be much larger than its mean in absolute value.
The mean is $O(\varepsilon)$ ; the typical size is of order $\sqrt\varepsilon$ .
The mean is of order $\sqrt\varepsilon$ . $\mathbb EO^\varepsilon_i(t)=\int_0^t(p^\varepsilon_i(s)-\nu_i(s))\,ds$ , and Module 4 bounds the integrand by $|\psi_0(s/\varepsilon)|+O(\varepsilon)$ .

The mean comes from the forward equation: an initial layer of width

\varepsilon

and an

O(\varepsilon)

outer error. The spread comes from the randomness of the path and is of order

\sqrt\varepsilon

. The explorer's sample mean of

n^\varepsilon(1)

is close to 0.

The explorer starts each chain from $\nu$ . If instead every chain started in state 1, would the limit of $n^\varepsilon(1)$ change?

Yes: the mean would shift by an amount of order 1.How long does the effect of the starting state last, and how much does it contribute to $O^\varepsilon$ ?
No: the starting state matters only for a time of order $\varepsilon$ , which changes $O^\varepsilon$ by $O(\varepsilon)$ and $n^\varepsilon$ by $O(\sqrt\varepsilon)$ .

This is why the theorem needs no assumption on

p^0

. The initial layer is invisible after the

\sqrt\varepsilon

scaling.

Practice

The deviation matrixA closed form

Let $Q$ be irreducible and $\Pi=\mathbf 1\nu$ . Show that $\int_0^\infty(e^{Qr}-\Pi)\,dr=(\Pi-Q)^{-1}-\Pi$ , and that this matrix $Z$ satisfies $Z\mathbf 1=0$ and $\nu Z=0$ .

Hint

$Q\Pi=\Pi Q=0$ and $\Pi^2=\Pi$ . Compute $e^{(Q-\Pi)r}$ .

Solution

Since $Q$ and $\Pi$ commute, $e^{(Q-\Pi)r}=e^{Qr}e^{-\Pi r}=e^{Qr}\big(I+(e^{-r}-1)\Pi\big)=e^{Qr}+(e^{-r}-1)\Pi$ , using $e^{Qr}\Pi=\Pi$ . The eigenvalues of $Q-\Pi$ are $-1$ (on $\mathbf 1$ ) and the nonzero eigenvalues of $Q$ , all with negative real part, so $\int_0^\infty e^{(Q-\Pi)r}dr=(\Pi-Q)^{-1}$ . Hence $\int_0^\infty(e^{Qr}-\Pi)\,dr=\int_0^\infty\big(e^{(Q-\Pi)r}-e^{-r}\Pi\big)dr=(\Pi-Q)^{-1}-\Pi$ . Finally $(e^{Qr}-\Pi)\mathbf 1=0$ and $\nu(e^{Qr}-\Pi)=0$ for every $r$ .

A tail estimateUsing the Gaussian limit

With $\lambda=\mu=1$ and $\varepsilon=0.01$ as in the worked example, estimate the probability that the chain spends more than 0.55 time units in state 1 during $[0,1]$ .

Solution

We need $O^\varepsilon_1(1)>0.05$ , i.e. $n^\varepsilon_1(1)>0.5$ . With limiting standard deviation $\sqrt{1/4}=0.5$ , that is one standard deviation: probability $\approx0.16$ .

SimulationYin–Zhang, Ch. 5: Occupation Measures: Asymptotic Properties and Ramification

The aggregated process

With weak and strong interactions, a single sample path splits into two parts. Which group the chain is in behaves more and more like a Markov chain with generator $\bar Q$ . Where it is inside the group becomes a fast blur, whose time fractions are given by $\nu^k$ .

The idea

Aggregated process.

\bar\alpha^\varepsilon(t)=k

when

\alpha^\varepsilon(t)\in\mathcal M_k

. For fixed

\varepsilon

it is usually not Markov: how soon it leaves a group depends on the state inside the group.

Weak convergence. Under the assumptions of Module 5,

\bar\alpha^\varepsilon(\cdot)

converges weakly (in the space of right-continuous paths with left limits) to a Markov chain

\bar\alpha(\cdot)

generated by

\bar Q(t)

.

Occupation inside a group. For every state

s_{kj}

, uniformly in

t\le T

,

\mathbb E\Big(\int_0^t\big(I_{\{\alpha^\varepsilon(s)=s_{kj}\}}-\nu^k_j(s)\,I_{\{\bar\alpha^\varepsilon(s)=k\}}\big)\,ds\Big)^2=O(\varepsilon).

Time in

s_{kj}

is

\nu^k_j

times time in

\mathcal M_k

, along each path, up to

O(\sqrt\varepsilon)

.

The centering uses the random indicator $I_{\{\bar\alpha^\varepsilon=k\}}$ , not a deterministic function: the time spent in a group stays random in the limit, because $\bar\alpha$ itself is random. Scaled by $1/\sqrt\varepsilon$ , these occupation errors converge to a switching diffusion, a Brownian motion whose variance rate depends on the current value of $\bar\alpha(t)$ ; its one-time distributions are Gaussian mixtures rather than Gaussian.

Explore

Time-scale parameter ε

The four-state chain of Module 5: $\bar q_{12}=0.5$ , $\bar q_{21}=0.6$ , $\nu^1=(\tfrac13,\tfrac23)$ , $\nu^2=(\tfrac34,\tfrac14)$ .

Top: which state

\alpha^\varepsilon(t)

is in (blue band: group

\mathcal M_1

; gold band:

\mathcal M_2

). Bottom: the group

\bar\alpha^\varepsilon(t)

.

Worked example

In the explorer's chain, show that $\bar\alpha^\varepsilon$ is not a Markov chain for fixed $\varepsilon$ , and explain why it becomes one in the limit.

Through which states can the chain enter $\mathcal M_1$ ?

The slow moves into $\mathcal M_1$ are $s_{21}\to s_{12}$ (rate 0.4) and $s_{22}\to s_{11}$ (rate 1.2). So it enters at $s_{12}$ or at $s_{11}$ , depending on where it was in $\mathcal M_2$ .
Just after entering, what is the rate of leaving $\mathcal M_1$ ?

0.9 if it entered at $s_{11}$ , 0.3 if at $s_{12}$ . The future of $\bar\alpha^\varepsilon$ depends on more than its present value 1, so $\bar\alpha^\varepsilon$ is not Markov.
How long does the chain remember its entry point?

The block $\tilde Q^1/\varepsilon$ has nonzero eigenvalue $-3/\varepsilon$ . After a time of a few times $\varepsilon/3$ the position inside $\mathcal M_1$ has distribution close to $\nu^1$ , whatever the entry point.
What does a stay in $\mathcal M_1$ look like in the limit?

Over a stay of order 1, the exit rate is averaged against $\nu^1$ and equals $\bar q_{12}=\tfrac13\cdot0.9+\tfrac23\cdot0.3=0.5$ . The memory lasts only $O(\varepsilon)$ , so stays become exponential with mean 2, independent of the past: the Markov property returns in the limit.

Check your understanding

On a fixed interval $[0,T]$ , as $\varepsilon\to0$ :

both the number of jumps and the number of group switches grow like $1/\varepsilon$ .Group switches are driven by $\hat Q$ , whose rates do not depend on $\varepsilon$ .
the number of jumps grows like $1/\varepsilon$ but the number of group switches stays of order 1.
both stay of order 1.Inside a group the rates are of order $1/\varepsilon$ .

Compare the first two statistics in the explorer as you move ε: jumps explode, switches do not.

Why is time in $s_{kj}$ compared with $\nu^k_j\,I_{\{\bar\alpha^\varepsilon(s)=k\}}$ rather than with $\nu^k_j\vartheta_k(s)$ ?

The two give the same limit. $\int_0^t(I_{\{\alpha^\varepsilon=s_{kj}\}}-\nu^k_j\vartheta_k)\,ds$ does not go to 0: the time spent in $\mathcal M_k$ is random even in the limit.
The time the path spends in $\mathcal M_k$ stays random as $\varepsilon\to0$ , so only the split inside the group becomes deterministic.
Because $\vartheta_k$ is not known explicitly.It is: it solves $\vartheta'=\vartheta\bar Q$ . The issue is randomness, not computability.

Contrast Module 7: with a single irreducible group there is nothing slow left to be random, and the centering is deterministic.

When is $\bar\alpha^\varepsilon$ exactly a Markov chain for every $\varepsilon$ ?

Never.It can be, under a condition on $\hat Q$ .
When, for each pair of groups $k\neq j$ , every state of $\mathcal M_k$ has the same total slow rate into $\mathcal M_j$ .
When every $\tilde Q^k$ is irreducible.The explorer's blocks are irreducible, and $\bar\alpha^\varepsilon$ still is not Markov there.

This is lumpability. Then the exit rate does not depend on the position inside the group, and

\bar\alpha^\varepsilon

is Markov with generator

\bar Q

, for every

\varepsilon

.

Practice

Long-run time in one stateCombining the two scales

For the explorer's chain, find the limit as $\varepsilon\to0$ of the long-run fraction of time spent in $s_{12}$ .

Solution

The aggregated chain spends a fraction $\bar q_{21}/(\bar q_{12}+\bar q_{21})=6/11$ of the time in $\mathcal M_1$ , and inside $\mathcal M_1$ a fraction $\nu^1_2=2/3$ is in $s_{12}$ . Answer: $4/11\approx0.364$ , the same as in the worked example of Module 5.

Lumpability exactlyChecking the Markov property

Change the explorer's slow rates so that $\bar\alpha^\varepsilon$ becomes Markov for every $\varepsilon$ while $\bar Q$ stays the same.

Solution

Give both states of $\mathcal M_1$ slow rate 0.5 into $\mathcal M_2$ (say $s_{11}\to s_{21}$ at 0.5 and $s_{12}\to s_{22}$ at 0.5), and both states of $\mathcal M_2$ rate 0.6 into $\mathcal M_1$ . The averages with $\nu^1$ and $\nu^2$ are still 0.5 and 0.6, and now the exit rate does not depend on the state inside the group.

Drag the payoffYin–Zhang, Ch. 6: Asymptotic Expansions of Solutions for Backward Equations

Backward equations and the terminal layer

An expected payoff collected at a final time $T$ solves the backward equation, which runs from $T$ back toward 0. The two-time-scale picture holds in reverse: a thin layer just before $T$ where values inside a group merge, then a slow evolution governed by the same $\bar Q$ as before.

The idea

Backward equation. For a payoff vector

f

(a column), the conditional expectations

u^\varepsilon(t,i)=\mathbb E\big[f(\alpha^\varepsilon(T))\mid\alpha^\varepsilon(t)=i\big]

form the column vector

u^\varepsilon(t)=P^\varepsilon(T,t)f

, which solves

\frac{du^\varepsilon(t)}{dt}=-Q^\varepsilon(t)\,u^\varepsilon(t),\qquad u^\varepsilon(T)=f.

A running reward

g

adds

-g(t)

to the right side.

Expansion.

u^\varepsilon(t)=\sum_{i\le n}\varepsilon^i\varphi_i(t)+\sum_{i\le n}\varepsilon^i\psi_i\big((T-t)/\varepsilon\big)+O(\varepsilon^{n+1})

: an outer part, and a terminal layer in the stretched backward time

\tau=(T-t)/\varepsilon

.

Leading terms (block-diagonal

\tilde Q

, groups

\mathcal M_k

).

\varphi_0(t)=\tilde{\mathbf 1}\,\bar u(t)

: constant on each group, with

\frac{d\bar u}{dt}=-\bar Q(t)\,\bar u(t),\qquad \bar u_k(T)=\nu^k(T)f^k=\sum_j\nu^k_j(T)f(s_{kj}),

and the terminal layer

\psi_0(\tau)=e^{\tilde Q(T)\tau}\big(f-\tilde{\mathbf 1}\bar u(T)\big)

, which decays exponentially in

\tau

.

Duality. Forward equations act on row vectors (distributions), and the leading term has the form $\vartheta\,\mathrm{diag}(\nu^k)$ . Backward equations act on column vectors (functions), and the leading term has the form $\tilde{\mathbf 1}\bar u$ . In both, $\bar Q=\mathrm{diag}(\nu^k)\hat Q\tilde{\mathbf 1}$ governs the slow evolution. Value functions of control problems are backward objects, so Chapters 7 to 10 rest on this expansion.

Explore

Time-scale parameter ε

Payoff presets:

Drag the four dots at

t=T

to set

f

.

s_{11}

s_{12}

s_{21}

s_{22}

; bands:

\bar u_1,\bar u_2

The terminal layer: the last

2.5\,\varepsilon

before

T

, in units of

\varepsilon

. Dashed:

\tilde{\mathbf 1}\bar u+\psi_0

.

Worked example

Use the four-state chain of Module 5 with $f=I_{\{s_{11}\}}$ , so $u^\varepsilon(t,i)=P(\alpha^\varepsilon(T)=s_{11}\mid\alpha^\varepsilon(t)=i)$ . Find the leading-order answer for $t$ well before $T$ .

What does the terminal layer do to $f=(1,0,0,0)$ ?

Inside $\mathcal M_1$ , the values 1 and 0 merge into their $\nu^1$ -average $\tfrac13\cdot1+\tfrac23\cdot0=\tfrac13$ ; inside $\mathcal M_2$ both are 0. So $\bar u(T)=(\tfrac13,0)$ .
Solve $d\bar u/dt=-\bar Q\bar u$ backward from $T$ , writing $s=T-t$ .

$\bar u(t)=e^{\bar Qs}\bar u(T)$ . With $\bar q_{12}=0.5$ and $\bar q_{21}=0.6$ , the first column of $e^{\bar Qs}$ is $\frac1{1.1}\big(0.6+0.5e^{-1.1s},\ 0.6-0.6e^{-1.1s}\big)$ , so

$\bar u_1=\tfrac13\cdot\frac{0.6+0.5e^{-1.1s}}{1.1},\qquad \bar u_2=\tfrac13\cdot\frac{0.6-0.6e^{-1.1s}}{1.1}.$
What do both tend to as $T-t\to\infty$ ? Interpret.

$\tfrac13\cdot\tfrac{0.6}{1.1}=\tfrac2{11}\approx0.182$ , the long-run probability of $s_{11}$ found in Module 5. Far from $T$ , the starting state is forgotten.
Is there a terminal layer when $f$ is constant on each group, say $f=(0,0,1,1)$ ?

No: $f-\tilde{\mathbf 1}\bar u(T)=0$ , so $\psi_0=0$ . The answer is $P(\bar\alpha(T)=2\mid\bar\alpha(t)=k)$ to leading order. Try the "indicator of $\mathcal M_2$ " preset: the curves leave $T$ without any fast transient.

Check your understanding

Near $t=T$ the values $u^\varepsilon(t,s_{11})$ and $u^\varepsilon(t,s_{12})$ differ. Going back from $T$ , how long until they agree to within 1%?

A time of order $\varepsilon$ (here about $\tfrac{\varepsilon}{3}\ln100$ ).
A time of order 1.The merging is done by the fast part $\tilde Q/\varepsilon$ , not the slow part.
They never agree: their payoffs differ.The payoffs differ only at $T$ . Starting slightly earlier, the chain has time to mix inside its group.

The difference is multiplied by

e^{-3(T-t)/\varepsilon}

, 3 being the mixing rate of block

\tilde Q^1

. Look at the zoomed panel.

Which average of the payoff over a group enters $\bar u(T)$ ?

The plain average over the states of the group.States are not visited equally often by the fast chain.
The $\nu^k$ -weighted average $\sum_j\nu^k_jf(s_{kj})$ .
The value at the state the chain is in at time $t$ .That state is forgotten within a time of order $\varepsilon$ .

At time

T

the chain is, given its group, distributed close to

\nu^k

, whatever state it was in a little earlier.

Why is the boundary layer at $T$ for the backward equation, but at 0 for the forward equation?

Each equation has a layer where its data are given: $p^0$ at time 0 for the forward equation, $f$ at time $T$ for the backward one.
Because $Q^\varepsilon$ changes fastest near $T$ .The layer appears even for constant generators. It is about the side condition.
The backward equation has layers at both ends.Its condition is imposed only at $T$ ; at $t=0$ nothing has to be matched.

In both cases the data need not have the shape of the leading term (

\vartheta\,\mathrm{diag}(\nu^k)

or

\tilde{\mathbf 1}\bar u

), and the layer is where the mismatch is absorbed.

Practice

A maximum principleBounds on expected payoffs

Show that $\min_if(i)\le u^\varepsilon(t,j)\le\max_if(i)$ for all $t\le T$ and all $j$ , in two ways: probabilistically, and from the matrix form $u^\varepsilon(t)=P^\varepsilon(T,t)f$ .

Solution

Probabilistically, $u^\varepsilon(t,j)$ is an expectation of $f(\alpha^\varepsilon(T))$ , which lies between the min and max. In matrix form, each row of $P^\varepsilon(T,t)$ is a probability vector, so $u^\varepsilon(t,j)$ is a convex combination of the $f(i)$ . In the explorer, no curve ever leaves the range of the four dots.

Probability of ending in a groupNo layer at all

For $f=(0,0,1,1)$ in the explorer's chain, find $\bar u(t)$ and the leading-order value of $u^\varepsilon(0,s_{11})$ when $T=3$ .

Solution

$\bar u(T)=(0,1)$ , and $\bar u=e^{\bar Qs}(0,1)^{\mathsf T}$ with $s=T-t$ : $\bar u_1=\frac{0.5-0.5e^{-1.1s}}{1.1}$ , $\bar u_2=\frac{0.5+0.6e^{-1.1s}}{1.1}$ . At $t=0$ , $s=3$ and $e^{-3.3}\approx0.0369$ , so $u^\varepsilon(0,s_{11})\approx\frac{0.5\cdot0.9631}{1.1}\approx0.438$ .

CalculatorYin–Zhang, Ch. 7: Markov Decision Problems

Markov decision problems on two time scales

When a controller can change the rates, the optimal cost solves a dynamic-programming equation on the full state space. On two time scales, a smaller problem with one state per group gets the value right to within $O(\varepsilon)$ , and its decisions give a policy that is nearly optimal for the original problem.

The idea

The problem. The generator depends on an action:

Q^\varepsilon(u)=\tilde Q(u)/\varepsilon+\hat Q(u)

, with

\tilde Q(u)

block diagonal as before. A stationary feedback policy picks an action

u(i)\in\Gamma

in each state. Its cost is

J^\varepsilon(i,u)=\mathbb E\int_0^\infty e^{-\rho t}\,G\big(\alpha^\varepsilon(t),u(\alpha^\varepsilon(t))\big)\,dt,\qquad \alpha^\varepsilon(0)=i,

and

v^\varepsilon(i)=\min_uJ^\varepsilon(i,u)

. For a fixed policy,

J^\varepsilon=(\rho I-Q^\varepsilon(u))^{-1}G(u)

.

Dynamic programming.

\displaystyle\rho\,v^\varepsilon(i)=\min_{a\in\Gamma}\Big\{G(i,a)+\sum_{j\neq i}q^\varepsilon_{ij}(a)\big(v^\varepsilon(j)-v^\varepsilon(i)\big)\Big\}.

Limit problem. Its states are the groups. In group

k

a decision is a whole tuple

U^k=(u^{k1},\dots,u^{km_k})

, one action per state of the group. It fixes the fast block

\tilde Q^k(U^k)

and hence its quasi-stationary distribution

\nu^k(U^k)

, which gives the averaged cost and rates

\bar G(k,U^k)=\sum_j\nu^k_j(U^k)\,G(s_{kj},u^{kj}),\qquad \bar Q(U)=\mathrm{diag}\big(\nu^1(U^1),\dots,\nu^l(U^l)\big)\,\hat Q(U)\,\tilde{\mathbf 1}.

Let

v(k)

be the optimal value of this

l

-state problem.

Asymptotic optimality.

v^\varepsilon(i)\to v(k)

for

i\in\mathcal M_k

. With finitely many actions,

v^\varepsilon(i)-v(k)=O(\varepsilon)

, and the policy

u_c(s_{kj})=u^{kj}_*

read off an optimal limit decision satisfies

J^\varepsilon(i,u_c)-v^\varepsilon(i)=O(\varepsilon)

.

Explore

s₁₁fast (÷ε): →s₁₂: 2 / 1slow: →s₂₁: 0.2 / 0.05cost rate: 0 / 2

Your action

s₁₂fast (÷ε): →s₁₁: 1 / 3slow: →s₂₂: 1 / 0.5cost rate: 1 / 4

Your action

s₂₁fast (÷ε): →s₂₂: 1 / 1slow: →s₁₁: 0.1 / 1cost rate: 4 / 6.5

Your action

s₂₂fast (÷ε): →s₂₁: 2 / 4slow: →s₁₂: 0.3 / 1.5cost rate: 6 / 7.5

Your action

Each entry reads action 0 / action 1. Group $\mathcal M_1=\{s_{11},s_{12}\}$ is the cheap regime, $\mathcal M_2$ the expensive one. Action 1 ("maintain") costs more; it slows the drift into $\mathcal M_2$ or speeds the return. Discount rate $\rho=0.5$ .

Time-scale parameter ε

Solid: optimal

v^\varepsilon

for

s_{11}

s_{12}

s_{21}

s_{22}

. Dashed: your policy. Dotted lines: the limit values

v(1),v(2)

. ε is on a log scale. Every value is computed exactly by a linear solve.

Worked example

Solve the explorer's limit problem for the decision $U^1=(0,1)$ , $U^2=(0,1)$ , and compare it with $U^1=(0,0)$ , $U^2=(0,1)$ .

With $U^1=(0,1)$ , what are $\nu^1$ , $\bar G(1)$ and $\bar q_{12}$ ?

Fast rates: $s_{11}\to s_{12}$ at 2, $s_{12}\to s_{11}$ at 3, so $\nu^1=(\tfrac35,\tfrac25)$ . Then $\bar G(1)=\tfrac35\cdot0+\tfrac25\cdot4=1.6$ and $\bar q_{12}=\tfrac35\cdot0.2+\tfrac25\cdot0.5=0.32$ .
Same for group 2 with $U^2=(0,1)$ .

$s_{21}\to s_{22}$ at 1, $s_{22}\to s_{21}$ at 4, so $\nu^2=(\tfrac45,\tfrac15)$ , $\bar G(2)=\tfrac45\cdot4+\tfrac15\cdot7.5=4.7$ and $\bar q_{21}=\tfrac45\cdot0.1+\tfrac15\cdot1.5=0.38$ .
Solve $\rho v=\bar G+\bar Qv$ .

$0.82v_1-0.32v_2=1.6$ and $-0.38v_1+0.88v_2=4.7$ . The determinant is $0.82\cdot0.88-0.32\cdot0.38=0.6$ , so

$v_1=\frac{1.6\cdot0.88+0.32\cdot4.7}{0.6}\approx4.853,\qquad v_2=\frac{0.82\cdot4.7+0.38\cdot1.6}{0.6}\approx7.437.$
Now try $U^1=(0,0)$ with the same $U^2$ . Is it better?

$\nu^1=(\tfrac13,\tfrac23)$ , $\bar G(1)=\tfrac23$ , $\bar q_{12}=\tfrac13\cdot0.2+\tfrac23\cdot1=\tfrac{11}{15}\approx0.733$ . Solving gives $v=(5.00,7.50)$ : cheaper to run, but the chain drifts into the expensive group more than twice as fast. Checking all 16 decisions shows $(0,1),(0,1)$ is optimal, so $u_c=(0,1,0,1)$ , and the explorer confirms it is exactly optimal for the full problem once $\varepsilon\lesssim0.4$ .

Check your understanding

Why is a decision in the limit problem a whole tuple $U^k$ , one action per state of the group, rather than a single action for the group?

Different states of the group can use different actions, and the mix determines $\nu^k$ , hence both $\bar G$ and $\bar Q$ .
Because the limit problem has more states than the original.It has fewer: one per group.
A single action per group would give the same value.Not here: the optimal $u_c=(0,1,0,1)$ uses both actions in each group.

The controller still sees the fast state, so it can act differently in each one. In the limit the effect of that local policy is summarised by

\nu^k(U^k)

.

At $\varepsilon=1$ the optimal policy is $(0,0,1,1)$ , not $u_c=(0,1,0,1)$ . Does this contradict the theorem?

Yes: $u_c$ should be optimal.The theorem promises near-optimality as $\varepsilon\to0$ , not optimality at a given $\varepsilon$ .
No: the theorem says $J^\varepsilon(u_c)-v^\varepsilon=O(\varepsilon)$ , which at $\varepsilon=1$ allows a visible gap.

At

\varepsilon=1

the "fast" moves are not fast and the chain does not mix inside a group, so the averaging behind the limit problem is rough. In the explorer the gap vanishes once

\varepsilon

is below about 0.4.

As $\varepsilon\to0$ , $v^\varepsilon(s_{11})-v^\varepsilon(s_{12})\to0$ , although the running costs in those states differ. Why?

Because the optimal actions in the two states are equal.In the limit policy they are not: 0 in $s_{11}$ , 1 in $s_{12}$ .
The chain moves between $s_{11}$ and $s_{12}$ within a time of order $\varepsilon$ , so the starting state affects the cost only over a time of order $\varepsilon$ .
Because discounting makes all values equal.The values of the two groups stay different, about 4.85 and 7.44.

So

v^\varepsilon(i)

depends at leading order only on the group of

i

, which is why the limit problem has one state per group.

Practice

Discounted costs are averagesThe resolvent

Show that $(\rho I-Q)^{-1}=\int_0^\infty e^{-\rho t}e^{Qt}\,dt$ for a generator $Q$ and $\rho>0$ . Conclude that it has nonnegative entries with row sums $1/\rho$ , so $\rho J$ is a weighted average of running costs.

Solution

$e^{-\rho t}e^{Qt}=e^{(Q-\rho I)t}$ , and $\frac{d}{dt}e^{(Q-\rho I)t}=(Q-\rho I)e^{(Q-\rho I)t}$ . Integrating over $[0,\infty)$ (the integrand decays since $e^{Qt}$ is stochastic) gives $(Q-\rho I)\int_0^\infty e^{(Q-\rho I)t}dt=-I$ . Entries are nonnegative because $e^{Qt}$ is, and row sums are $\int_0^\infty e^{-\rho t}dt=1/\rho$ .

The long-run average cost of

u_c

Letting ρ → 0

In the limit, what is the long-run average cost per unit time of the policy $u_c$ in the explorer?

Solution

The aggregated chain under $u_c$ has $\bar q_{12}=0.32$ , $\bar q_{21}=0.38$ , so it spends fractions $\tfrac{0.38}{0.70}$ and $\tfrac{0.32}{0.70}$ of the time in the two groups. The average cost is $\tfrac{0.38\cdot1.6+0.32\cdot4.7}{0.70}=\tfrac{2.112}{0.70}\approx3.017$ .

SimulationYin–Zhang, Ch. 8: Stochastic Control of Dynamical Systems

Production planning with an unreliable machine

A machine breaks down and is repaired at random, and its output must meet a steady demand. A good policy keeps a safety stock, the hedging point. When breakdowns and repairs are fast compared with the time scale of the costs, the right safety stock is of order $\varepsilon$ , and the problem reduces to a deterministic one in which the machine has its average capacity.

The idea

Controlled system. A state

x^\varepsilon(t)\in\R^n

driven by a two-time-scale chain:

\dot x^\varepsilon=f(x^\varepsilon,\alpha^\varepsilon(t),u(t))

, with

u(t)\in\Gamma

chosen from the past of

\alpha^\varepsilon

, and cost

J^\varepsilon=\mathbb E\int_0^\infty e^{-\rho t}G(x^\varepsilon(t),\alpha^\varepsilon(t),u(t))\,dt

(or over a finite horizon). The book assumes

f

affine in

u

,

G

convex, and

\tilde Q

block diagonal with irreducible blocks.

Limit problem. Take one irreducible fast block, for simplicity. A limit control is a tuple

U=(u^1,\dots,u^m)

, one control per machine state, and everything is averaged against

\nu

:

\dot x=\sum_i\nu_i\,f(x,i,u^i),\qquad \text{cost }\int_0^\infty e^{-\rho t}\sum_i\nu_i\,G(x,i,u^i)\,dt.

With several groups, the averaging is done within each group, and the limit system still switches between groups according to

\bar Q

.

Near optimality. The value functions converge:

|v^\varepsilon(x,\alpha)-v^0(x)|\le K(1+|x|^\kappa)\sqrt\varepsilon

for

\alpha

in the group. A control built from a limit control

U(t)

as

u^\varepsilon(t)=\sum_iI_{\{\alpha^\varepsilon(t)=i\}}u^i(t)

is asymptotically optimal. The occupation-measure bounds of Module 7 are what make the comparison work.

The example. One machine, up or down: breakdown rate $\lambda/\varepsilon$ , repair rate $\mu/\varepsilon$ . Surplus (inventory minus backlog) $\dot x=u-z$ , with $0\le u\le k$ when up and $u=0$ when down; cost $c^+x^++c^-x^-$ per unit time, long-run average. A hedging-point policy with level $\theta$ produces at full rate below $\theta$ , exactly $z$ at $\theta$ , and nothing above. This is the classical setting where such policies are known to be optimal.

Average capacity is $k\mu/(\lambda+\mu)$ , and the limit problem is a reliable machine of that capacity. If that exceeds $z$ , the limit plan keeps $x=0$ forever at zero cost.

Explore

Time-scale parameter ε Hedging point θ

$\lambda=1,\ \mu=2$ (up two thirds of the time), capacity $k=3$ , demand $z=1.5$ , inventory cost $c^+=1$ , backlog cost $c^-=5$ .

The surplus

x(t)

, starting from 0. Shaded: machine down.

exact

J(\theta)

your

\theta

best

\theta^*

. Ring:

J(\theta)

estimated by simulation.

Worked example

Find the long-run behaviour of the surplus under a hedging-point policy, the optimal $\theta$ , and how it scales with $\varepsilon$ . Assume $k\mu/(\lambda+\mu)>z$ .

Which levels does $x$ visit in the long run, and with which velocities?

Once below $\theta$ , $x$ never goes above it. Up and below $\theta$ : $\dot x=k-z>0$ . Up at $\theta$ : $\dot x=0$ . Down: $\dot x=-z$ . So the stationary law lives on $(-\infty,\theta]$ , with an atom at $\theta$ (machine up and waiting).
Let $f_{\rm up}(x),f_{\rm down}(x)$ be the stationary densities below $\theta$ . Why is $(k-z)f_{\rm up}=z\,f_{\rm down}$ , and what form must they have?

In equilibrium the net flow of probability across any level $x<\theta$ is zero: mass crosses upward at rate $(k-z)f_{\rm up}(x)$ and downward at rate $z f_{\rm down}(x)$ . The balance equation for the up state, $\frac{d}{dx}\big[(k-z)f_{\rm up}\big]=-\frac\lambda\varepsilon f_{\rm up}+\frac\mu\varepsilon f_{\rm down}$ , then gives $f\propto e^{\beta(x-\theta)}$ with

$\beta=\frac1\varepsilon\Big(\frac\mu z-\frac\lambda{k-z}\Big)>0.$

So $x-\theta$ has a law that does not depend on $\theta$ : an atom at 0 and an exponential tail of mean depth of order $1/\beta$ , proportional to $\varepsilon$ .
Differentiate $J(\theta)=\mathbb E\,h(\theta+Y)$ , where $Y=x-\theta$ and $h(x)=c^+x^++c^-x^-$ . What is the optimality condition?

Since the law of $Y$ is free of $\theta$ , $J'(\theta)=\mathbb E\,h'(\theta+Y)=c^+P(x>0)-c^-P(x<0)$ ( $x$ has no atom at 0 when $\theta>0$ ). So at the optimum

$P(x<0)=\frac{c^+}{c^++c^-},$

a "critical fractile": backlog should occur a fraction $1/6$ of the time with $c^+=1,\ c^-=5$ .
With the explorer's numbers, find $\theta^*$ .

$\beta=\frac1\varepsilon\big(\tfrac2{1.5}-\tfrac1{1.5}\big)=\frac2{3\varepsilon}$ . Normalising the densities and the atom shows $P(x<\theta)=\tfrac23$ for every $\theta$ and $\varepsilon$ , so $P(x<0)=\tfrac23e^{-\beta\theta}$ . Setting this to $\tfrac16$ : $e^{-\beta\theta}=\tfrac14$ , so

$\theta^*=\frac{\ln4}{\beta}=1.5\,\varepsilon\ln4\approx2.08\,\varepsilon.$

Both $\theta^*$ and the optimal cost are proportional to $\varepsilon$ , and tend to the limit problem's answers, 0 and 0.

Check your understanding

Keep $\theta=1$ fixed and let $\varepsilon\to0$ . What does the long-run cost $J(1)$ tend to?

0, the value of the limit problem.Only if $\theta$ shrinks too. With $\theta=1$ you keep one unit of stock.
$c^+\cdot1=1$ : excursions below $\theta$ are only $O(\varepsilon)$ deep, so $x\approx1$ nearly all the time.
$\infty$ .The average capacity exceeds demand, so the surplus is stable for every $\theta$ .

Set

\theta=1

in the explorer and slide ε down: the cost at the gold dot approaches 1 while the minimum of the curve slides toward

\theta=0

.

In the limit problem, why does the capacity enter only through $k\nu_{\rm up}=k\mu/(\lambda+\mu)$ ?

Because the planner should prepare for the worst case, capacity 0.Worst-case planning would ignore that the machine is up most of the time.
The machine state switches so fast that production averages out: over any interval of order 1 it is up a fraction $\nu_{\rm up}$ of the time.
Because the policy produces at full speed $k$ all the time.It cannot produce when the machine is down.

The dynamics are linear in

u

, so

\int(u-z)\,dt

depends on the path of

\alpha^\varepsilon

only through occupation times, which concentrate at

\nu

(Module 7).

In the explorer, doubling $\varepsilon$ doubles $\theta^*$ . Why?

Down periods last $O(\varepsilon)$ , so the dips below $\theta$ are $O(\varepsilon)$ deep, and the safety stock that guards against them scales the same way.
Because the cost rates $c^\pm$ are proportional to $\varepsilon$ .The cost rates do not change with $\varepsilon$ .
It is a coincidence of these parameter values.It holds for any parameters: the whole stationary law of $x-\theta$ is a rescaling by $\varepsilon$ .

Rescaling time by

\varepsilon

maps the problem with parameter

\varepsilon

onto the one with

\varepsilon=1

, with

x

scaled by

\varepsilon

. Watch "θ* ÷ ε" in the stats: it is always

1.5\ln4\approx2.079

.

Practice

When is the surplus stable?The sign of β

Show that $\beta>0$ exactly when the average capacity exceeds demand, $k\mu/(\lambda+\mu)>z$ . What happens to the surplus otherwise?

Solution

$\varepsilon\beta=\frac{\mu(k-z)-\lambda z}{z(k-z)}$ , and $\mu(k-z)-\lambda z>0\iff\mu k>(\lambda+\mu)z$ . If the inequality fails, the average drift below $\theta$ , $\nu_{\rm up}(k-z)-\nu_{\rm down}z=\nu_{\rm up}k-z$ , is $\le0$ ; the surplus drifts to $-\infty$ (or is null recurrent at equality) and the long-run cost is infinite.

The limit problem from a backlogA deterministic plan

In the limit problem with the explorer's numbers, start from $x(0)=-2$ and ignore discounting over the transient. What is the optimal plan, how long does it take to clear the backlog, and what does it cost?

Solution

Produce at the full average rate $k\nu_{\rm up}=2$ until $x=0$ : $\dot x=0.5$ , so it takes 4 time units. Then hold $x=0$ , producing $z/\nu_{\rm up}=2.25\le3$ whenever the machine is up. The cost is $c^-\int_0^4(2-t/2)\,dt=5\cdot4=20$ .

SimulationYin–Zhang, Ch. 9: Numerical Methods for Control and Optimization

Tuning a threshold by stochastic approximation

When the cost of a threshold policy has no formula, simulate it. A stochastic-approximation recursion uses noisy estimates of the slope of the cost, taken from short simulated runs, and walks the threshold toward the optimum. Here it is on the machine of Module 11, where the exact answer is known, so you can see how well it does.

The idea

Solving the HJB equation numerically. One route to an optimal control is to discretise the dynamic-programming (HJB) equation on a grid of mesh

\Delta

, replacing each derivative by a one-sided difference taken in the direction of the drift. The discrete equations are themselves the dynamic-programming equations of a controlled Markov chain on the grid, and their solutions converge to the value function as

\Delta\to0

(in the viscosity sense).

Threshold policies. Often the control is of threshold type with parameters

\theta

: the action depends only on which side of some levels the state lies, as with the hedging point. The long-run average cost

J(\theta)

becomes a function of finitely many numbers, to be minimised.

Stochastic approximation. Update

\theta_{n+1}=\theta_n-\eta\,g_n

, where

g_n

is an estimate of

\nabla J(\theta_n)

computed from the simulated system over the

n

-th time window. Here, with window length

W

and a small

\delta

,

g_n=\frac{\hat J_n(\theta_n+\delta)-\hat J_n(\theta_n-\delta)}{2\delta},

where

\hat J_n

is the time-average cost over the window. Both runs use the same machine history (common random numbers): it does not depend on

\theta

, so only the effect of the policy differs.

ODE method. Interpolate the iterates in time steps of

\eta

. As

\eta\to0

the interpolated path converges weakly to the solution of

\dot\theta=-\nabla J(\theta)

(up to the bias of the gradient estimate), and over long times it concentrates near a stable point

\theta^*

of that ODE. With a constant step the iterates keep fluctuating around

\theta^*

; the fluctuations shrink as

\eta

does.

Explore

Step size η

Window W

Common random numbers

Start θ₀

The machine of Module 11 with $\varepsilon=1$ ; finite-difference step $\delta=0.25$ . Each step simulates the next window of length $W$ twice, at $\theta_n+\delta$ and at $\theta_n-\delta$ , continuing from where the previous window ended.

\theta_n

exact

\theta^*

Exact

J(\theta)

, with the last 30 iterates (darker is more recent).

Press Run. With η = 0.3 and W = 20 the iterates reach θ* = 2.08 within a few dozen steps and then jitter around it (on average slightly to the right, where J is flatter; that offset shrinks with η). Turn common random numbers off: the two cost estimates now have independent noise, the slope estimates become much noisier, and the iterates wander, mostly into the flat region to the right of θ*. A short window W = 5 also makes each slope estimate noisier. A smaller η is steadier but slower.

Worked example

The finite difference has its own bias: even with perfect cost estimates, the recursion settles where the difference quotient vanishes, not where $J'$ does. Compute that point for the explorer's machine ( $\varepsilon=1$ ).

From Module 11, write $J'(\theta)$ for $\theta>0$ .

$J'(\theta)=c^+-(c^++c^-)P(x<0)=1-6\cdot\tfrac23e^{-\beta\theta}=1-4e^{-\beta\theta}$ , with $\beta=\tfrac23$ .
For $\theta\ge\delta$ , compute the exact central difference $\frac{J(\theta+\delta)-J(\theta-\delta)}{2\delta}$ .

Integrate $J'$ over $[\theta-\delta,\theta+\delta]$ and divide by $2\delta$ :

$\frac{J(\theta+\delta)-J(\theta-\delta)}{2\delta}=1-4e^{-\beta\theta}\,\frac{\sinh\beta\delta}{\beta\delta}.$
Where does it vanish?

At $e^{-\beta\theta}=\tfrac14\cdot\frac{\beta\delta}{\sinh\beta\delta}$ , i.e.

$\theta_\delta=\theta^*+\frac1\beta\ln\frac{\sinh\beta\delta}{\beta\delta}\approx\theta^*+\frac{\beta\delta^2}{6}.$
Evaluate the bias for $\delta=0.25$ and $\delta=1$ .

$\delta=0.25$ : $\beta\delta=\tfrac16$ , bias $\approx0.0069$ , invisible next to the noise. $\delta=1$ : $\beta\delta=\tfrac23$ , $\sinh(\tfrac23)/\tfrac23\approx1.0757$ , bias $\approx1.5\ln1.0757\approx0.109$ . A large $\delta$ gives quieter estimates but a biased answer: the usual trade-off of finite differences.

Check your understanding

Why run $\theta_n+\delta$ and $\theta_n-\delta$ on the same machine history?

Otherwise the gradient estimate would be biased.Independent runs give an estimate with the same mean. The difference is in its variance.
The noise in the two cost estimates is then largely shared and cancels in the difference, so the slope estimate has much smaller variance.
To save simulation time.Both runs still have to be simulated; only the machine history is shared.

\mathrm{Var}(\hat J_+-\hat J_-)=\mathrm{Var}\hat J_++\mathrm{Var}\hat J_--2\,\mathrm{Cov}

, and sharing the randomness makes the covariance large. Toggle it off in the explorer and compare.

With a constant step size $\eta$ , what do the iterates do in the long run?

Converge to $\theta^*$ exactly.Each step adds noise of size $\eta\,\times$ (noise in $g_n$ ), which never stops.
Fluctuate in a neighbourhood of $\theta^*$ whose size shrinks as $\eta\to0$ .
Drift away, because the noise accumulates.Near $\theta^*$ the mean step $-\eta J'(\theta)$ pulls them back.

To converge exactly you need decreasing steps, typically

\sum_n\eta_n=\infty

and

\sum_n\eta_n^2<\infty

. A constant step can track a slowly changing optimum, which decreasing steps cannot.

Which ODE describes the iterates for small $\eta$ , ignoring the finite-difference bias?

$\dot\theta=-J'(\theta)$ , the gradient flow of the long-run cost.
$\dot\theta=-\hat J_n'(\theta)$ , with the noisy estimates.In the limit the noise averages out; the ODE involves only its mean.
$\dot\theta=-J(\theta)$ .The recursion moves against the slope, not against the cost itself.

Over a stretch of

1/\eta

iterations, many independent windows are averaged. This is the same averaging idea as in the two-time-scale chains: here the "fast" variable is the simulated system and the "slow" one is

\theta

.

Practice

A pathwise derivativeNo finite differences needed

Start the surplus at $x(0)=\theta$ and keep the machine history fixed. Show that $x_\theta(t)-\theta$ does not depend on $\theta$ , and derive the unbiased slope estimate $\frac1W\int_0^W\big(c^+I_{\{x(t)>0\}}-c^-I_{\{x(t)<0\}}\big)\,dt$ for $\frac{d}{d\theta}\hat J$ .

Solution

The dynamics of $y=x-\theta$ involve only $y$ and the machine state (full speed when $y<0$ , match demand at $y=0$ ), and $y(0)=0$ . So $x_\theta(t)=\theta+y(t)$ with $y$ free of $\theta$ . Then $\hat J(\theta)=\frac1W\int_0^Wh(\theta+y(t))\,dt$ , and differentiating under the integral gives $\frac1W\int_0^Wh'(x(t))\,dt$ , which is the formula. This is infinitesimal perturbation analysis; it uses one run instead of two.

A Markov chain approximationDiscretising a drift

For $\dot x=b$ with constant $b$ , build a continuous-time chain on the grid $\Delta\Z$ with the right mean drift and a variance per unit time that vanishes as $\Delta\to0$ .

Solution

Jump $x\to x+\Delta$ at rate $b^+/\Delta$ and $x\to x-\Delta$ at rate $b^-/\Delta$ . The mean drift is $\Delta(b^+-b^-)/\Delta=b$ and the variance per unit time is $\Delta^2|b|/\Delta=|b|\Delta\to0$ . With the machine state added, the grid chain jumps also as $\alpha$ does, and its DP equation is the upwind discretisation of the HJB equation.

InteractiveYin–Zhang, Ch. 10: Hybrid LQG Problems

LQ control with Markov switching

A linear system whose coefficients switch with a Markov chain is best controlled by a linear feedback whose gain depends on the current regime. The gains come from coupled Riccati equations, one per state of the chain. On two time scales, they collapse to one averaged Riccati equation per group.

The idea

Hybrid LQG problem.

dx=\big(A(\alpha(t))x+B(\alpha(t))u\big)dt+\sigma\,dw

on

[s,T]

, with

\alpha

independent of the Brownian motion

w

, and cost

J=\mathbb E\Big[\int_s^T\big(x^{\mathsf T}M(\alpha)x+u^{\mathsf T}N(\alpha)u\big)\,dt+x(T)^{\mathsf T}Dx(T)\Big],

with

M(i)\ge0

and

N(i),D>0

.

Coupled Riccati equations. The value is

v(s,x,i)=x^{\mathsf T}K(s,i)x+q(s,i)

, where for each state

i

\dot K_i=-K_iA_i-A_i^{\mathsf T}K_i-M_i+K_iB_iN_i^{-1}B_i^{\mathsf T}K_i-\sum_{j}q_{ij}K_j,\qquad K_i(T)=D,

and the optimal control is

u=-N_i^{-1}B_i^{\mathsf T}K_i\,x

in regime

i

. The term

\sum_jq_{ij}K_j=\sum_{j\neq i}q_{ij}(K_j-K_i)

is the generator acting on the value function: it accounts for future switches.

Two-time-scale limit. With

q_{ij}=q^\varepsilon_{ij}

from

\tilde Q/\varepsilon+\hat Q

,

K^\varepsilon(s,s_{kj})\to\bar K(s,k)

for every state of group

k

, where

\bar K

solves the Riccati system of the averaged problem:

\dot{\bar K}_k=-\bar K_k\bar A_k-\bar A_k^{\mathsf T}\bar K_k-\bar M_k+\bar K_k\,\overline{BN^{-1}B^{\mathsf T}}_k\,\bar K_k-\sum_l\bar q_{kl}\bar K_l,

each bar being a

\nu^k

-average over the group. The control

u^\varepsilon=-N^{-1}(\alpha^\varepsilon)B^{\mathsf T}(\alpha^\varepsilon)\bar K(t,\bar\alpha^\varepsilon(t))\,x

is nearly optimal: its cost approaches the optimal one as

\varepsilon\to0

.

Only the Riccati part is shared within a group. The feedback still uses the $B$ and $N$ of the current state.

Explore

Time-scale parameter ε

Scalar system, the four-state chain of Module 5. $a=(1.2,-0.6,0.5,-1)$ , $b=(1,0.5,1.5,1)$ , $m=n=1$ , $D=0.5$ , $T=3$ . Averages: $\bar a=(0,\ 0.125)$ , $\overline{b^2}=(0.5,\ 1.9375)$ .

Thin lines:

K^\varepsilon(t,i)

for the four states. Wide bands:

\bar K(t,1)

(blue) and

\bar K(t,2)

(gold).

\max_{t,i}|K^\varepsilon-\bar K|

against

\varepsilon

, log–log; dashed: slope 1. Gold dot: the current ε.

The four Riccati equations are integrated backward from

K(T)=D

. With ε large they spread out, each regime pricing its own dynamics; as ε shrinks the two states of each group are pulled together by the

1/\varepsilon

coupling, and the pair follows the averaged equation. The error falls along the slope-1 line: it is

O(\varepsilon)

. Note state

s_{11}

alone is unstable (

a=1.2

), but its group averages to

\bar a=0

.

Worked example

A scalar system switches very fast between two regimes, $a_1=2$ and $a_2=-1$ , spending fractions $\nu=(\tfrac13,\tfrac23)$ of the time in them; $b=m=n=1$ in both, $D=0$ . Find the limiting Riccati solution and its long-horizon gain, and check the closed loop.

Write the scalar Riccati equations for the two regimes.

$\dot K_i=-2a_iK_i-1+K_i^2-\frac1\varepsilon\sum_{j}\tilde q_{ij}K_j$ , $K_i(T)=0$ .
Why are $K_1$ and $K_2$ nearly equal for small $\varepsilon$ ?

The coupling term $\frac1\varepsilon\tilde q_{12}(K_2-K_1)$ is huge unless $K_1\approx K_2$ . Any difference is damped at rate $O(1/\varepsilon)$ , so $K_1-K_2=O(\varepsilon)$ .
Average the two equations with weights $\nu$ to get an equation for the common value $\bar K$ .

The $1/\varepsilon$ terms cancel because $\nu\tilde Q=0$ . What remains is $\dot{\bar K}=-2\bar a\bar K-1+\bar K^2$ with $\bar a=\tfrac13\cdot2+\tfrac23\cdot(-1)=0$ : $\dot{\bar K}=\bar K^2-1$ , $\bar K(T)=0$ . Its solution is $\bar K(t)=\tanh(T-t)$ .
Far from $T$ the gain is $\bar K\approx1$ , so $u\approx-x$ . Is the closed loop stable?

In regime 1 the closed loop is $\dot x=(2-1)x=x$ , unstable on its own. In regime 2 it is $\dot x=-2x$ . With fast switching, $x(t)\approx x(0)\exp\int_0^t(a(\alpha^\varepsilon(s))-1)\,ds\approx x(0)e^{(\bar a-1)t}=x(0)e^{-t}$ : stable. The averaged design stabilises the switched system though it would not stabilise regime 1 alone.

Check your understanding

Where does the term $\sum_jq_{ij}K_j$ in the Riccati equation for regime $i$ come from?

From the Brownian noise $\sigma\,dw$ .The noise only enters the equation for $q(s,i)$ , through $\mathrm{tr}(\sigma\sigma^{\mathsf T}K_i)$ .
From the chance of switching to regime $j$ : it is the generator of $\alpha$ acting on the value function $x^{\mathsf T}K(s,\cdot)x$ .
From the terminal cost $D$ . $D$ enters only through the terminal condition.

In the HJB equation the chain contributes

\sum_{j\neq i}q_{ij}\big(v(s,x,j)-v(s,x,i)\big)

. With

v=x^{\mathsf T}K x+q

this is

x^{\mathsf T}\big(\sum_jq_{ij}K_j\big)x+\dots

, since

\sum_jq_{ij}=0

.

In the nearly optimal feedback $u^\varepsilon=-N^{-1}(\alpha^\varepsilon)B^{\mathsf T}(\alpha^\varepsilon)\bar K(t,\bar\alpha^\varepsilon)x$ , which part no longer depends on the state inside the group?

$\bar K$ , the Riccati solution; $B$ and $N$ are still those of the current state.
All of it: the feedback depends only on the group.The controller still knows the current state and uses its $B$ and $N$ .
None of it.The Riccati part is shared within a group, as the explorer shows.

If

B

and

N

happen to be the same within each group, the whole feedback depends only on the group. Otherwise only the value-function part is shared.

To form the averaged equation, a student averages $b$ first and then squares: $(\sum_j\nu_jb_j)^2/n$ . What is wrong?

Nothing; it gives the same answer.In the explorer's group 2, $(0.75\cdot1.5+0.25\cdot1)^2=1.89$ , while $0.75\cdot2.25+0.25\cdot1=1.94$ .
The control enters the cost quadratically, so the right average is $\overline{b^2/n}=\sum_j\nu_jb_j^2/n_j$ , which is larger by Jensen's inequality.
One should average $K$ instead of the coefficients. $K$ is the unknown being solved for; the coefficients are what get averaged.

In each state the best control removes

b_j^2K^2/n_j

from the cost rate, and the fast chain averages those reductions. Equivalently, the limit control is a tuple

U=(u^{k1},\dots)

and the averaged control effect is

\sum_j\nu_jB_ju^{kj}

.

Practice

From HJB to RiccatiThe scalar case

For a scalar system with one regime, substitute $v=K(s)x^2+q(s)$ into $0=v_s+\min_u\{(ax+bu)v_x+mx^2+nu^2\}+\tfrac12\sigma^2v_{xx}$ and derive the equations for $K$ and $q$ .

Solution

The minimum over $u$ is at $u=-bv_x/(2n)=-bKx/n$ , with minimum value $axv_x+mx^2-b^2K^2x^2/n$ . Matching the $x^2$ terms: $\dot K+2aK+m-b^2K^2/n=0$ ; the constant terms: $\dot q+\sigma^2K=0$ . Terminal conditions $K(T)=D$ , $q(T)=0$ .

The long-horizon gainAn algebraic Riccati equation

For one scalar regime, find the limit of $K(t)$ as $T-t\to\infty$ , and check it against the worked example.

Solution

Set $\dot K=0$ : $\frac{b^2}nK^2-2aK-m=0$ , whose positive root is $K=\frac{n}{b^2}\Big(a+\sqrt{a^2+b^2m/n}\Big)$ . With $a=0$ and $b=m=n=1$ , $K=1$ , the limit of $\tanh(T-t)$ . The closed loop rate is $a-b^2K/n=-\sqrt{a^2+b^2m/n}<0$ : always stable.