Concrete examples

Grocery shopping. Bob really likes to eat cheerios for breakfast every morning. Let the state X_t denotes the amount of Cheerios in Bob’s kitchen on day t. The action U_t denotes the amount of Cheerios Bob buys on day t and W_t denotes the amount of Cheerios he eats that day. The random variable W_t varies with Bob’s hunger. This yields the dynamical system X_{t+1} = X_t + U_t - W_t\,. While this example is a bit cartoonish, it turns out that such simple models are commonly used in managing large production supply chains. Any system where resources are stochastically depleted and must be replenished can be modeled comparably. If Bob had a rough model for how much he eats in a given day, he could forecast when his supply would be depleted. And he could then minimize the number of trips he’d need to make to the grocer using optimal control.

Moving objects. Consider a physical model of a flying object. The simplest model of the dynamics of this system are given by Newton’s laws of mechanics. Let Z_t denote the position of the vehicle (this is a three-dimensional vector). The derivative of position is velocity V_t = \frac{\partial Z_t}{\partial t}\,, and the derivative of velocity is acceleration A_t = \frac{\partial V_t}{\partial t}\,. Now, we can approximate the rules in discrete time with the simple Taylor approximations \begin{aligned} Z_{t+1} &= Z_t + \Delta V_t\\ V_{t+1} &= V_t + \Delta A_t \end{aligned} We also know by Newton’s second law, that acceleration is equal to the total applied force divided by the mass of the object: F=mA. The flying object will be subject to external forces such as gravity and wind W_t and it will also receive forces from its propellers U_t. Then we can add this to the equations to yield a model \begin{aligned} Z_{t+1} &= Z_t + \Delta V_t\\ V_{t+1} &= V_t + \frac{\Delta}{m} (W_t+U_t)\,. \end{aligned}

Oftentimes, we only observe the acceleration through an accelerometer. Then estimating the position and velocity becomes a filtering problem. Optimal control problems for this model include flying the object along a given trajectory or flying to a desired location in a minimal amount of time.

Markov decision processes

Our definition in terms of structural equations is not the only dynamical model used in machine learning. Some people prefer to work directly with probabilistic transition models and conditional probabilities. In a Markov Decision Process, we again have a state X_t and input U_t, and they are linked by a probabilistic model \mathop\mathbb{P}[X_{t+1}\mid X_t,U_t]\,. This is effectively the same as the structural equation model above except we hide the randomness in this probabilistic notation.

Example: machine repair. The following example illustrates the elegance of the conditional probability models for dynamical systems. This example comes from BertsekasBertsekas, Dynamic Programming and Optimal Control, 4th ed., vol. 1 (Athena Scientific, 2017).. Suppose we have a machine with ten states of repair. State 10 denotes excellent condition and 1 denotes the inability to function. Every time one uses the machine in state j, it has a probability of falling into disrepair, given by the probabilities \mathop\mathbb{P}[X_{t+1} = i\mid X_t=j]\,, where \mathop\mathbb{P}[X_{t+1}=i\mid X_t=j]=0 if i>j. The action a one can take at any time is to repair the machine, resetting the system state to 10. Hence \mathop\mathbb{P}[X_{t+1}=i\mid X_{t}=j,U_t=0] = \mathop\mathbb{P}[X_{t+1}=i\mid X_{t}=j] and \mathop\mathbb{P}[X_{t+1}=i\mid X_{t}=j,U_t=1] = \mathbb{1}\left\{ i=10 \right\} \,. While we could write this dynamical system as a structural equation model, it is more conveniently expressed by these probability tables.

Infinite time horizons and stationary policies

The Q-functions we derived for these finite time horizons are time varying. One applies a different policy for each step in time. However, on long horizons with time invariant dynamics and costs, we can get a simpler formula. First, for example, consider the limit: \begin{array}{ll} \text{maximize} & \lim_{N\rightarrow \infty} \mathop\mathbb{E}_{W_t}[ \frac{1}{N} \sum_{t=0}^N R(X_t,U_t,W_t) ]\\ \text{subject to} & X_{t+1} = f(X_t, U_t, W_t),~U_t=\pi_t(X_t)\\ & \text{($x_0$ given).} \end{array} Such infinite time horizon problems are referred to as average cost dynamic programs. Note that there are no subscripts on the rewards or transition functions in this model.

Average cost dynamic programming is deceptively difficult. These formulations are not directly amenable to standard dynamic programming techniques except in cases with special structure. A considerably simpler infinite time formulation is known as discounted dynamic programming, and this is the most popular studied formulation. Discounting is a mathematical convenience that dramatically simplifies algorithms and analysis. Consider the SDM problem \begin{array}{ll} \text{maximize} & (1-\gamma) \mathop\mathbb{E}_{W_t}[ \sum_{t=0}^\infty \gamma^t R(X_t,U_t,W_t) ]\\ \text{subject to} & X_{t+1} = f(X_t, U_t, W_t),~U_t=\pi_t(X_t)\\ & \text{($x_0$ given).} \end{array} where \gamma is a scalar in (0,1) called the discount factor. For \gamma close to 1, the discounted reward is approximately equal to the average reward. However, unlike the average cost model, the discounted cost has particularly clean optimality conditions. If we define \mathcal{Q}_\gamma(x,u) to be the Q-function obtained from solving the discounted problem with initial condition x, then we have a discounted version of dynamic programming, now with the same Q-functions on the left and right hand sides: \mathcal{Q}_\gamma(x,u) = \mathop\mathbb{E}_{W} \left[ R(x,u,W) + \gamma \max_{u'} \mathcal{Q}_\gamma(f(x,u,W),u')\right]\,. The optimal policy is now for all times to let u_t= \arg\max_u \mathcal{Q}_\gamma(x_t,u)\,. The policy is time invariant and one can execute it without any knowledge of the reward or dynamics functions. At every stage, one simply has to maximize a function to find the optimal action. Foreshadowing to the next chapter, the formula additionally suggest that the amount that needs to be “learned” in order to “control” is not very large for these infinite time horizon problems.

Tabular MDPs

Tabular MDPs refer to Markov Decision Processes with small number of states and actions. Say that there are S states and A actions. Then the transition rules are given by tables of conditional probabilities \mathop\mathbb{P}[X_{t+1}|X_t,U_t], and the size of such tables are S^2 A. The Q-functions for the tabular case are also tables, each of size SA, enumerating the cost-to-go for all possible state action pairs. In this case, the maximization \max_{u'}\mathcal{Q}_{a\rightarrow b}(x,u') corresponds to looking through all of the actions and choosing the largest entry in the table. Hence, in the case that the rewards are deterministic functions of x and u, Bellman’s equation simplifies to \mathcal{Q}_{t\rightarrow T}(x,u) = R_t(x,u) + \sum_{x'} \mathop\mathbb{P}[X_{t+1}=x'|X_t=x,U_t=u] \max_{u'} \mathcal{Q}_{t+1\rightarrow T}(x',u')\,. This function can be computed by elementary matrix-vector operations: the Q-functions are S \times A arrays of numbers. The “max” operation can be performed by operating over each row in such an array. The summation with respect to x' can be implemented by multiplying a SA \times S array by an S-dimensional vector. We complete the calculation by summing the resulting expression with the S\times A array of rewards. Hence, the total time to compute \mathcal{Q}_{t\rightarrow T} is O(S^2A).

Linear quadratic regulator

The other important problem where dynamic programming is efficiently solvable is the case when the dynamics are linear and the rewards are quadratic. In control design, this class of problems is generally referred to as the problem of the Linear Quadratic Regulator (LQR): \begin{array}{ll} \text{minimize} \, & \mathop\mathbb{E}_{W_t} \left[\frac{1}{2}\sum_{t=0}^T X_t^T \Phi_t X_t + U_t^T \Psi_t U_t\right], \\ \text{subject to} & X_{t+1} = A_t X_t+ B_t U_t + W_t,~U_t=\pi_t(X_t) \\ & \text{($x_0$ given).} \end{array} Here, \Phi_t and \Psi_t are most commonly positive semidefinite matrices. w_t is noise with zero mean and bounded variance, and we assume W_t and W_{t'} are independent when t\neq t'. The state transitions are governed by a linear update rule with A_t and B_t appropriately sized matrices. We also abide by the common convention in control textbooks to pose the problem as a minimization—not maximization—problem.

As we have seen above, many systems can be modeled by linear dynamics in the real world. However, we haven’t yet discussed cost functions. It’s important to emphasize here that cost functions are designed not given. Recall back to supervised learning: though we wanted to minimize the number of errors made on out-of-sample data, on in-sample data we minimized convex surrogate problems. The situation is exactly the same in this more complex world of dynamical decision making. Cost functions are designed by the engineer so that the SDM problems are tractable but also so that the desired outcomes are achieved. Cost function design is part of the toolkit for online decision making, and quadratic costs can often yield surprisingly good performance for complex problems.

Quadratic costs are also attractive for computational reasons. They are convex as long as \Phi_t and \Psi_t are positive definite. Quadratic functions are closed under minimization, maximization, addition. And for zero mean noise W_t with covariance \Sigma, we know that the noise interacts nicely with the cost function. That is, we have \mathbb{E}_W[(x+W)^T M (x+W)] = x^T M x + \operatorname{Tr}(M\Sigma) for any vector x and matrix M. Hence, when we run dynamic programming, every Q-function is necessarily quadratic. Moreover, since the Q-functions are quadratic, the optimal action is a linear function of the state U_t = -K_t X_t for some matrix K_t.

Now consider the case where there are static costs \Phi_t=\Phi and \Psi_t=\Psi, and time invariant dynamics such that A_t=A and B_t=B for all t. One can check that the Q-function on a finite time horizon satisfies a recursion \mathcal{Q}_{t \rightarrow T}(x,u) = x^T \Phi x + u^T \Psi u + (Ax+Bu)^T M_{t+1} (Ax+Bu) + c_t\,. for some positive definite matrix M_{t+1}. In the limit as the time horizon tends to infinity, the optimal control policy is static, linear state feedback: u_t = -K x_t\,. Here the matrix K is defined by K=(\Psi + B^T M B)^{-1} B^T M A and M is a solution to the Discrete Algebraic Riccati Equation M = \Phi + A^T M A - (A^T M B)(\Psi + B^T M B)^{-1} (B^T M A)\,. Here, M is the unique solution of the Riccati equation where all of the eigenvalues of A-BK have magnitude less than 1. Finding this specific solution is relatively easy using standard linear algebraic techniques. It is also the limit of the Q-functions computed above.

Policy and value iteration

Two of the most well studied methods for solving such discounted infinite time horizon problems are value iteration and policy iteration. Value iteration proceeds by the steps \mathcal{Q}_{k+1} (x,u) = \mathop\mathbb{E}_{W} \left[ R(x,u,W) + \gamma \max_{u'} \mathcal{Q}_k (f(x,u,W),u')\right]\,. That is, it simply tries to solve the Bellman equation by running a fixed point operation. This method succeeds when the iteration is a contraction mapping, and this occurs in many contexts.

On the other hand, Policy Iteration is a two step procedure: policy evaluation followed by policy improvement. Given a policy \pi_k, the policy evaluation step is given by \mathcal{Q}_{k+1} (x,u) = \mathop\mathbb{E}\left[ R(x,u,W) + \gamma \mathcal{Q}_k (f(x,\pi_k(x),W),\pi_k(x))\right]\,. And then the policy is updated by the rule \pi_{k+1}(x) = \arg\max_u \mathcal{Q}_{k+1} (x,u)\,. Often times, several steps of policy evaluation are performed before updating the policy.

For both policy and value iteration, we need to be able to compute expectations efficiently and must be able to update all values of x and u in the associated Q functions. This is certainly doable for tabular MDPs. For general low dimensional problems, policy iteration and value iteration can be approximated by gridding state space, and then treating the problem as a tabular one. Then, the resulting Q function can be extended to other (x,u) pairs by interpolation. There are also special cases where the maxima and minima yield closed form solutions and hence these iterations reduce to simpler forms. LQR is a canonical example of such a situation.

Model predictive control

If the Q-functions in value or policy iteration converge quickly, long-term planning might not be necessary, and we can effectively solve infinite horizon problem with short-term planning. This is the key idea behind one of the most powerful techniques for efficiently and effectively finding quality policies for SDM problems called model predictive control.

Suppose that we aim to solve the infinite horizon average reward problem: \begin{array}{ll} \text{maximize} & \lim_{T\rightarrow \infty} \mathop\mathbb{E}_{W_t}[\frac{1}{T}\sum_{t=0}^T R_t(W_t,U_t) ]\\ \text{subject to} & X_{t+1} = f_t(X_t, U_t, W_t)\\ & U_t = \pi_t(X_t)\\ & \text{($x_0$ given).} \end{array} Model Predictive Control computes an open loop policy on a finite horizon H \begin{array}{ll} \text{maximize}_{u_t} & \mathop\mathbb{E}_{W_t}[ \sum_{t=0}^H R_t(X_t,u_t)]\\ \text{subject to} & X_{t+1} = f_t(X_t, u_t, W_t)\\ & \text{($X_0$ = x).} \end{array} This gives a sequence u_0(x),\ldots, u_H(x). The policy is then set to be \pi(x)=u_0(x). After this policy is executed, we observe a new state, x', based on the dynamics. We then recompute the optimization, now using x_0=x' and setting the action to be \pi(x')= u_0(x').

MPC is a rather intuitive decision strategy. The main idea is to plan out a sequence of actions for a given horizon, taking into account as much uncertainty as possible. But rather than executing the entire sequence, we play the first action and then gain information from the environment about the noise. This direct feedback influences the next planning stage. For this reason, model predictive control is often a successful control policy even when implemented with inaccurate or approximate models. Model Predictive Control also allows us to easily add a variety of constraints to our plan at little cost, such as bounds on the magnitude of the actions. We just append these to the optimization formulation and then lean on the computational solver to make us a plan with these constraints.

To concretely see how Model Predictive Control can be effective, it’s helpful to work through an example. Let’s suppose the dynamics and rewards are time invariant. Let’s suppose further that the reward function is bounded above, and there is some state-action pair (x_\star,u_\star) which achieves this maximal reward R_\mathrm{max}.

Suppose we solve the finite time horizon problem where we enforce that (x,u) must be at (x_\star,u_\star) at the end of the time horizon: \begin{array}{ll} \text{maximize} & \mathop\mathbb{E}_{W_t}[ \sum_{t=0}^H R(X_t,u_t)]\\ \text{subject to} & X_{t+1} = f(X_t, u_t, W_t)\\ & (X_H,U_H) = (x_\star,u_\star)\\ & \text{($X_0$ = x).} \end{array} We replan at every time step by solving this optimization problem and taking the first action.

The following proposition summarizes how this policy performs

The proposition asserts that there is a “burn in” cost associated with the initial horizon length. This term goes to zero with T, but will have different values for different H. The policy converges to a residual average cost due to the stochasticity of the problem and the fact that we try to force the system to the state (x_\star,u_\star).

The main caveat with this argument is that there may not exist a policy that drives x to x_\star from an arbitrary initial condition and any realization of the disturbance signal. Much of the analysis of MPC schemes is devoted to guaranteeing that the problems are recursively feasible, meaning that such constraints can be met for all time.

This example also shows how it is often helpful to have some sort of recourse at the end of the planning horizon to mitigate the possibility of being too greedy and driving the system into a bad state. The terminal condition of forcing x_H=0 adds an element of safety to the planning, and ensures stable execution for all time. More general, adding some terminal condition to the planning horizon \mathcal{C}(x_{H}) is part of good Model Predictive Control design and is often a powerful way to balance performance and robustness.

Optimal filtering

Given a model of the state transition function and the observation function, we can attempt to compute the maximum a posteriori estimate of the state from the data. That is, we could compute p(x_t | \tau_t) and then estimate the mode of this density. Here, we show that such an estimator has a relatively simple recursive formula, though it is not always computationally tractable to compute this formula.

To proceed, we first need a calculation that takes advantage of the conditional independence structure of our dynamical system model. Note that \begin{aligned} p(y_t, x_t, & x_{t-1}, u_{t-1} | \tau_{t-1} ) = \\ &p(y_t| x_t, x_{t-1}, u_{t-1}, \tau_{t-1} ) p(x_t | x_{t-1},u_{t-1}, \tau_{t-1} )\\ &\qquad\qquad\qquad\qquad\times p(x_{t-1}|\tau_{t-1} ) p(u_{t-1}|\tau_{t-1} )\\ &= p(y_t| x_t ) p(x_t | x_{t-1},u_{t-1} ) p(x_{t-1}|\tau_{t-1} ) p(u_{t-1}|\tau_{t-1} )\,. \end{aligned} This decomposition is into terms we now recognize. p(x_t | x_{t-1},u_{t-1} ) and p(y_t| x_t ) define the POMDP model and are known. p(u_t|\tau_t) is our policy and it’s what we’re trying to design. The only unknown here is p(x_{t-1}|\tau_{t-1} ), but this expression gives us a recursive formula to p(x_{t}|\tau_{t} ) for all t.

To derive this formula, we apply Bayes rule and then use the above calculation: \begin{aligned} p(x_t | \tau_t) &=\frac{\int_{x_{t-1}} p(x_t , y_t, x_{t-1}, u_{t-1} | \tau_{t-1}) } {\int_{x_t,x_{t-1}} p(x_t , y_t, x_{t-1}, u_{t-1}| \tau_{t-1})}\\ &=\frac{\int_{x_{t-1}} p(y_t| x_t ) p(x_t | x_{t-1},u_{t-1} ) p(x_{t-1}|\tau_{t-1} ) p(u_{t-1}|\tau_{t-1} ) } {\int_{x_{t},x_{t-1}} p(y_t| x_t ) p(x_t | x_{t-1},u_{t-1} ) p(x_{t-1}|\tau_{t-1} ) p(u_{t-1}|\tau_{t-1} )}\\ &=\frac{\int_{x_{t-1}} p(y_t| x_t ) p(x_t | x_{t-1},u_{t-1} ) p(x_{t-1}|\tau_{t-1} ) } {\int_{x_t,x_{t-1}} p(y_t| x_t ) p(x_t| x_{t-1},u_{t-1} ) p(x_{t-1}|\tau_{t-1} )}\,.\qquad (1) \end{aligned} Given a prior for x_0, this now gives us a formula to compute a MAP estimate of x_t for all t, incorporating data in a streaming fashion. For tabular POMDP models with small state spaces, this formula can be computed simply by summing up the conditional probabilities. In POMDPs without inputs—also known as hidden Markov models—this formula gives the forward pass of Viterbi’s decoding algorithm. For models where the dynamics are linear and the noise is Gaussian, these formulas reduce into an elegant closed form solution known as Kalman filtering. In general, this optimal filtering algorithm is called belief propagation and is the basis of a variety of algorithmic techniques in the field of graphical models.

Kalman filtering

For the case of linear dynamical systems, the above calculation has a simple closed form solution that looks similar to the solution of LQR. This estimator is called a Kalman filter, and is one of the most important tools in signal processing and estimation. The Kalman filter assumes a linear dynamical system driven by Gaussian noise with observations corrupted by Gaussian noise \begin{aligned} X_{t+1} &= A X_t+ B U_t + W_t\\ Y_{t} & = C X_t + V_t\,. \end{aligned} Here, assume W_t and V_t are independent for all time and Gaussian with means zero and covariances \Sigma_W and \Sigma_V respectively. Because of the joint Gaussianity, we can compute a closed form formula for the density of X_t conditioned on the past observations and actions, p(x_t | \tau_t). Indeed, X_t is a Gaussian itself.

On an infinite time horizon, the Kalman filter takes a simple and elucidating form: \begin{aligned} \hat{x}_{t+1} &= A \hat{x}_t + B u_t - L (y_t - \hat{y}_t)\\ \hat{y}_t &= C \hat{x}_t \end{aligned} where L=A P C^T (C P C^T + \Sigma_V)^{-1} and P is the positive semidefinite solution to the discrete algebraic Riccati equation P = A P A^T + \Sigma_W - (A P C^T) (CPC^T +\Sigma_V)^{-1} (C\Sigma A^T)\,.

The derivation of this form follows from the calculation in Equation 1. But the explanation of the formulas tend to be more insightful than the derivation. Imagine the case where L=0. Then our estimate \hat{x}_t simulates the same dynamics as our model with no noise corruption. The matrix L computes a correction for this simulation based on the observed y_t. This feedback correction is chosen in such a way such that P is the steady state covariance of the error X_t-\hat{x}_t. P ends up being the minimum variance possible with an estimator that is unbiased in the sense that \mathop\mathbb{E}[\hat{x}_t-X_t]=0.

Another interesting property of this calculation is that the L matrix is the LQR gain associated with the LQR problem \begin{array}{ll} \text{minimize} \, & \lim_{T\rightarrow \infty}\mathop\mathbb{E}_{W_t} \left[\frac{1}{2}\sum_{t=0}^T X_t^T \Sigma_W X_t + U_t^T \Sigma_V U_t\right], \\ \text{subject to} & X_{t+1} = A^T X_t+ B^T U_t + W_t,~U_t=\pi_t(X_t) \\ & \text{($x_0$ given).} \end{array} Control theorists often refer to this pairing as the duality between estimation and control.

Feedforward prediction

While the optimal filter can be computed in simple cases, we often do not have simple computational means to compute the optimal state estimate. That said, the problem of state estimation is necessarily one of prediction, and the first half of this course gave us a general strategy for building such estimators from data. Given many simulations or experimental measurements of our system, we can try to estimate a function h such that X_t \approx h(\tau_t). To make this concrete, we can look at time lags of the history \tau_{t-s\rightarrow t} := (y_t,\ldots, y_{t-s}, u_{t-1}, \ldots, u_{t-s})\,. Such time lags are necessarily all of the same length. Then estimating \begin{array}{ll} \text{minimize}_h & \sum_t \operatorname{loss}(h(\tau_{t-s\rightarrow t}), x_t) \end{array} is a supervised learning problem, and standard tools can be applied to design architectures for and estimate h.