Lecture 12: Markov decision processes (MDP)

Goal: mathematical model for optimal control with random dynamics and cost.

References

Chapter 2: Bertsekas, D. P. & Tsitsiklis, J. N. Neuro-dynamic programming. Athena Scientific (1996).
Chapter 3: Sutton, R. S. & Barto, A. G. Reinforcement Learning: An Introduction. (The MIT Press, 1998).

Mathematical model

Analogous to how Markov processes naturally generalize deterministic dynamics, we can try to generalize optimal control problems into problems with random dynamics and cost. Instead of

\[ \min_u c(x,u)\ \text{s.t.}\ x^+ = F(x,u) \]

we might want to write

\[ \min_u c(x,u)\ \text{s.t.}\ x^+ \sim P(x,u) \]

but there’s a problem: in the second equation \(c(x,u)\) does not take on a single value due to the fact that \(x^+\) does not take on a single value. We need a way to evaluate the cost of a distribution of trajectories. But we already have a tool to convert a probability distribution to a single number: statistics.

Code

from imports import *

# Dynamics: stable 2D linear system
A = np.array([[0.9, -0.3],
              [0.2,  0.85]])
x0 = np.array([1.0, 0.0])
T = 30
sigma = 0.15
N = 200  # number of stochastic trajectories

# --- Deterministic trajectory ---
x_det = np.zeros((T + 1, 2))
x_det[0] = x0
for t in range(T):
    x_det[t + 1] = A @ x_det[t]

# --- Stochastic trajectories ---
x_sto = np.zeros((N, T + 1, 2))
x_sto[:, 0, :] = x0
for t in range(T):
    noise = np.random.randn(N, 2) * sigma
    x_sto[:, t + 1, :] = (A @ x_sto[:, t, :].T).T + noise

# Compute median and IQR across trajectories
x_med = np.median(x_sto, axis=0)          # (T+1, 2)
x_q25 = np.percentile(x_sto, 25, axis=0)  # (T+1, 2)
x_q75 = np.percentile(x_sto, 75, axis=0)  # (T+1, 2)

time = np.arange(T + 1)
colors = ['tab:blue', 'tab:orange']
labels = [r'$x_1$', r'$x_2$']

fig, axes = plt.subplots(1, 2, sharey=True)

# --- Left: deterministic ---
ax = axes[0]
for i in range(2):
    ax.plot(time, x_det[:, i], color=colors[i], label=labels[i])
ax.set_title(r'$x^+ = F(x,u)$')
ax.set_xlabel(r'$t$')
ax.set_ylabel(r'$x$')
ax.legend()

# --- Right: stochastic distribution ---
ax = axes[1]
for i in range(2):
    ax.plot(time, x_med[:, i], color=colors[i], label=labels[i])
    ax.fill_between(time, x_q25[:, i], x_q75[:, i],
                    color=colors[i], alpha=0.3)
ax.set_title(r'$x^+ \sim P(x,u)$')
ax.set_xlabel(r'$t$')
ax.legend()

fig.tight_layout()
show_with_alt(fig,
    "Two side-by-side plots of x vs time for a stable 2D linear system. "
    "Left: deterministic trajectory x^+ = F(x,u) showing x1 and x2 converging smoothly to zero. "
    "Right: stochastic distribution x^+ ~ P(x,u) showing the median of 200 sample trajectories "
    "with interquartile range shaded, illustrating the spread introduced by process noise."
)

Q: what statistic(s) should we choose for the cost?

There are many valid choices. By far the most common is the expected cost \(\text{E}[c(x,u)]\). A key motivation for this choice is that it admits many nice theoretical results. If the system will run many times, and if there are no catastrophic outcomes (e.g. break a sensor or actuator), it is a reasonable objective.

However, if you care about the “spread” in the cost, or if some outcomes are catastrophic, considering \(\text{Var}[c(x,u)]\) may be a good idea, or perhaps the interquartile or other confidence intervals.

Expected cost

We’ll focus on minimizing expected cost \(\min_u \text{E}[c(s,u)]\ \text{s.t.}\ x^+\sim P(x,u)\) where

\[\begin{align*} c(x,u) & = \ell(t,x_t) + \sum_{s=0}^{t-1}\mathcal{L}(s,x_s,u_s) \quad \text{-- finite-horizon} \\ \text{or}\quad & = \lim_{t\rightarrow\infty} \frac{1}{t} \sum_{s=0}^{t-1}\mathcal{L}(s,x_s,u_s) \quad \text{-- infinite-horizon, time-averaged} \\ \text{or}\quad & = \lim_{t\rightarrow\infty} \sum_{s=0}^{t}\gamma^s\mathcal{L}(s,x_s,u_s) \quad \text{-- infinite-horizon, exponentially-discounted} \end{align*}\]

where \(\gamma\in(0,1)\) is termed a discount factor.

A bit of notation: we let \(T\) denote a time interval (e.g. \(T = [0,t]\) or \(= [0,\infty)\)) and \(X^T = \{x:T\rightarrow X\}\), \(U^T = \{u:T\rightarrow U\}\) denote the set of functions from times into states or inputs, respectively. (This comports with the notation \(\R^n\), as vectors \(v\in\R^n\) can be identified with functions \(v:\{1,\dots,n\}\rightarrow\R\).)

Definition: a Markov decision process (MDP) is specified by \((X,U,P,c)\) where \(X\) and \(U\) are sets, \(P:X\times U\rightarrow\Delta(X)\), and \(c:X^T\times U^T\rightarrow\R\).