2 Linear Quadratic Regulators

2.1Introduction¶

Up to this point, we have considered decision problems with finitely many states and actions. However, in many applications, states and actions may take on continuous values. For example, consider autonomous driving, controlling a robot’s joints, and automated manufacturing. How can we teach computers to solve these kinds of problems? This is the task of continuous control.

Aside from the change in the state and action spaces, the general problem setup remains the same: we seek to construct an optimal policy that outputs actions to solve the desired task. We will see that many key ideas and algorithms, in particular dynamic programming algorithms, carry over to this new setting.

This chapter introduces a fundamental tool to solve a simple class of continuous control problems: the linear quadratic regulator. We will then extend this basic method to more complex settings.

2.2Optimal control¶

Recall that an MDP is defined by its state space $\mathcal{S}$, action space $\mathcal{A}$, state transitions $P$, reward function $r$, and discount factor γ or time horizon $\hor$. These have equivalents in the control setting:

• The state and action spaces are continuous rather than finite. That is, $\mathcal{S} \subseteq \mathbb{R}^{n_\st}$ and $\mathcal{A} \subseteq \mathbb{R}^{n_\act}$, where $n_\st$ and $n_\act$ are the corresponding dimensions of these spaces, i.e. the number of coordinates to specify a single state or action respectively.

• We call the state transitions the dynamics of the system. In the most general case, these might change across timesteps and also include some stochastic noise $w_\hi$ at each timestep. We denote these dynamics as the function $f_\hi$ such that $\st_{\hi+1} = f_\hi(\st_\hi, \act_\hi, w_\hi)$. Of course, we can simplify to cases where the dynamics are deterministic/noise-free (no $w_\hi$ term) and/or time-homogeneous (the same function $f$ across timesteps).

• Instead of maximizing the reward function, we seek to minimize the cost function $c_\hi: \mathcal{S} \times \mathcal{A} \to \mathbb{R}$. Often, the cost function describes how far away we are from a target state-action pair $(\st^\star, \act^\star)$. An important special case is when the cost is time-homogeneous; that is, it remains the same function $c$ at each timestep $h$.

• We seek to minimize the undiscounted cost within a finite time horizon $\hor$. Note that we end an episode at the final state $\st_\hor$ -- there is no $\act_\hor$, and so we denote the cost for the final state as $c_\hor(\st_\hor)$.

With all of these components, we can now formulate the optimal control problem: compute a policy to minimize the expected undiscounted cost over $\hor$ timesteps. In this chapter, we will only consider deterministic, time-dependent policies $\pi = (\pi_0, \dots, \pi_{H-1})$ where $\pi_h : \mathcal{S} \to \mathcal{A}$ for each $\hi \in [\hor]$.

2.2.1A first attempt: Discretization¶

Can we solve this problem using tools from the finite MDP setting? If $\mathcal{S}$ and $\mathcal{A}$ were finite, then we’d be able to work backwards using the DP algorithm for computing the optimal policy in an MDP (Definition 1.11). This inspires us to try discretizing the problem.

Suppose $\mathcal{S}$ and $\mathcal{A}$ are bounded, that is, $\max_{\st \in \mathcal{S}} \|\st\| \le B_\st$ and $\max_{\act \in \mathcal{A}} \|\act\| \le B_\act$. To make $\mathcal{S}$ and $\mathcal{A}$ finite, let’s choose some small positive ε, and simply round each coordinate to the nearest multiple of ε. For example, if $\epsilon = 0.01$, then we round each element of $\st$ and $\act$ to two decimal spaces.

However, the discretized $\widetilde{\mathcal{S}}$ and $\widetilde{\mathcal{A}}$ may be finite, but they may be infeasibly large: we must divide each dimension into intervals of length $\varepsilon$, resulting in $|\widetilde{\mathcal{S}}| = (B_\st/\varepsilon)^{n_\st}$ and $|\widetilde{\mathcal{A}}| = (B_\act/\varepsilon)^{n_\act}$. To get a sense of how quickly this grows, consider $\varepsilon = 0.01, n_\st = n_\act = 10$. Then the number of elements in the transition matrix would be $|\widetilde{\mathcal{S}}|^2 |\widetilde{\mathcal{A}}| = (100^{10})^2 (100^{10}) = 10^{60}$! (That’s a trillion trillion trillion trillion trillion.)

What properties of the problem could we instead make use of? Note that by discretizing the state and action spaces, we implicitly assumed that rounding each state or action vector by some tiny amount $\varepsilon$ wouldn’t change the behavior of the system by much; namely, that the cost and dynamics were relatively continuous. Can we use this continuous structure in other ways? This leads us to the linear quadratic regulator.

2.3The Linear Quadratic Regulator¶

The optimal control problem Definition 2.1 seems highly complex in general. Is there a relevant simplification that we can analyze? The linear quadratic regulator (LQR) is a solvable case and a fundamental tool in control theory.

We will henceforth abbreviate “symmetric positive definite” as s.p.d. and “positive definite” as p.d.

It will be helpful to reintroduce the value function notation for a policy to denote the average cost it incurs. These will be instrumental in constructing the optimal policy via dynamic programming, as we did in Section 1.3.2 for MDPs.

2.4Optimality and the Riccati Equation¶

In this section, we’ll compute the optimal value function $V^\star_h$, Q-function $Q^\star_h$, and policy $\pi^\star_h$ in the linear quadratic regulator using dynamic programming in a very similar way to the DP algorithms in the MDP setting. Recall the definition of the optimal value function:

We will prove the striking fact that the solution has very simple structure: $V_h^\star$ and $Q^\star_h$ are upward-curved quadratics and $\pi_h^\star$ is linear and furthermore does not depend on the noise!

The construction (and inductive proof) proceeds similarly to the one in the MDP setting.

1. We’ll compute $V_\hor^\star$ (at the end of the horizon) as our base case.
2. Then we’ll work step-by-step backwards in time, using $V_{\hi+1}^\star$ to compute $Q_\hi^\star$, $\pi_{\hi}^\star$, and $V_\hi^\star$.

Base case: At the final timestep, there are no possible actions to take, and so $V^\star_\hor(\st) = c(\st) = \st^\top Q \st$. Thus $V_\hor^\star(\st) = \st^\top P_\hor \st + p_\hor$ where $P_\hor = Q$ and $p_\hor = 0$.

Inductive hypothesis: We seek to show that the inductive step holds for both theorems: If $V^\star_{\hi+1}(\st)$ is an upward-curved quadratic, then $V^\star_\hi(\st)$ must also be an upward-curved quadratic, and $\pi^\star_\hi(\st)$ must be linear. We’ll break this down into the following steps:

1. Show that $Q^\star_\hi(\st, \act)$ is an upward-curved quadratic (in both $\st$ and $\act$).
2. Derive the optimal policy $\pi^\star_\hi(\st) = \arg \min_\act Q^\star_\hi(\st, \act)$ and show that it’s linear.
3. Show that $V^\star_\hi(\st)$ is an upward-curved quadratic.

We first assume the inductive hypothesis that our theorems are true at time $\hi+1$. That is,

Now we’ve shown that $V^\star_\hi(\st) = \st^\top P_\hi \st + p_\hi$, where $P_\hi$ is s.p.d., proving the inductive hypothesis and completing the proof of Theorem 2.2 and Theorem 2.1.

In summary, we just demonstrated that at each timestep $\hi \in [\hor]$, the optimal value function $V^\star_\hi$ and optimal Q-function $Q^\star_\hi$ are both upward-curved quadratics and the optimal policy $\pi^\star_\hi$ is linear. We also showed that all of these quantities can be calculated using a sequence of s.p.d. matrices $P_0, \dots, P_H$ that can be defined recursively using the Riccati equation Definition 2.5.

Before we move on to some extensions of LQR, let’s consider how the state at time $\hi$ behaves when we act according to this optimal policy.

2.4.1Expected state at time $\hi$¶

How can we compute the expected state at time $\hi$ when acting according to the optimal policy? Let’s first express $\st_\hi$ in a cleaner way in terms of the history. Note that having linear dynamics makes it easy to expand terms backwards in time:

Let’s consider the average state at this time, given all the past states and actions. Since we assume that $\E [w_\hi] = 0$ (this is the zero vector in $d$ dimensions), when we take an expectation, the $w_\hi$ term vanishes due to linearity, and so we’re left with

This introdces the quantity $A - B K_i$, which shows up frequently in control theory. For example, one important question is: will $\st_\hi$ remain bounded, or will it go to infinity as time goes on? To answer this, let’s imagine for simplicity that these $K_i$s are equal (call this matrix $K$). Then the expression above becomes $(A-BK)^\hi \st_0$. Now consider the maximum eigenvalue $\lambda_{\max}$ of $A - BK$. If $|\lambda_{\max}| > 1$, then there’s some nonzero initial state $\bar \st_0$, the corresponding eigenvector, for which

Otherwise, if $|\lambda_{\max}| < 1$, then it’s impossible for your original state to explode as dramatically.

2.5Extensions¶

We’ve now formulated an optimal solution for the time-homogeneous LQR and computed the expected state under the optimal policy. However, real world tasks rarely have such simple dynamics, and we may wish to design more complex cost functions. In this section, we’ll consider more general extensions of LQR where some of the assumptions we made above are relaxed. Specifically, we’ll consider:

1. Time-dependency, where the dynamics and cost function might change depending on the timestep.

2. General quadratic cost, where we allow for linear terms and a constant term.

3. Tracking a goal trajectory rather than aiming for a single goal state-action pair.

Combining these will allow us to use the LQR solution to solve more complex setups by taking Taylor approximations of the dynamics and cost functions.

2.5.1Time-dependent dynamics and cost function¶

So far, we’ve considered the time-homogeneous case, where the dynamics and cost function stay the same at every timestep. However, this might not always be the case. As an example, in many sports, the rules and scoring system might change during an overtime period. To address these sorts of problems, we can loosen the time-homogeneous restriction, and consider the case where the dynamics and cost function are time-dependent. Our analysis remains almost identical; in fact, we can simply add a time index to the matrices $A$ and $B$ that determine the dynamics and the matrices $Q$ and $R$ that determine the cost.

The modified problem is now defined as follows:

The derivation of the optimal value functions and the optimal policy remains almost exactly the same, and we can modify the Riccati equation accordingly:

Additionally, by allowing the dynamics to vary across time, we gain the ability to locally approximate nonlinear dynamics at each timestep. We’ll discuss this later in the chapter.

2.5.2More general quadratic cost functions¶

Our original cost function had only second-order terms with respect to the state and action, incentivizing staying as close as possible to $(\st^\star, \act^\star) = (0, 0)$. We can also consider more general quadratic cost functions that also have first-order terms and a constant term. Combining this with time-dependent dynamics results in the following expression, where we introduce a new matrix $M_\hi$ for the cross term, linear coefficients $q_\hi$ and $r_\hi$ for the state and action respectively, and a constant term $c_\hi$:

Similarly, we can also include a constant term $v_\hi \in \mathbb{R}^{n_\st}$ in the dynamics (note that this is deterministic at each timestep, unlike the stochastic noise $w_\hi$):

2.5.3Tracking a predefined trajectory¶

Consider applying LQR to a task like autonomous driving, where the target state-action pair changes over time. We might want the vehicle to follow a predefined trajectory of states and actions $(\st_\hi^\star, \act_\hi^\star)_{\hi=0}^{\hor-1}$. To express this as a control problem, we’ll need a corresponding time-dependent cost function:

Note that this punishes states and actions that are far from the intended trajectory. By expanding out these multiplications, we can see that this is actually a special case of the more general quadratic cost function above (2.38):

2.6Approximating nonlinear dynamics¶

The LQR algorithm solves for the optimal policy when the dynamics are linear and the cost function is an upward-curved quadratic. However, real settings are rarely this simple! Let’s return to the CartPole example from the start of the chapter (Example 2.1). The dynamics (physics) aren’t linear. How can we approximate this by an LQR problem?

Concretely, let’s consider a noise-free problem since, as we saw, the noise doesn’t factor into the optimal policy. Let’s assume the dynamics and cost function are stationary, and ignore the terminal state for simplicity:

This is now only slightly simplified from the general optimal control problem (see Definition 2.1). Here, we don’t know an analytical form for the dynamics $f$ or the cost function $c$, but we assume that we’re able to query/sample/simulate them to get their values at a given state and action. To clarify, consider the case where the dynamics are given by real world physics. We can’t (yet) write down an expression for the dynamics that we can differentiate or integrate analytically. However, we can still simulate the dynamics and cost function by running a real-world experiment and measuring the resulting states and costs. How can we adapt LQR to this more general nonlinear case?

2.6.1Local linearization¶

How can we apply LQR when the dynamics are nonlinear or the cost function is more complex? We’ll exploit the useful fact that we can take a function that’s locally continuous around $(s^\star, a^\star)$ and approximate it nearby with low-order polynomials (i.e. its Taylor approximation). In particular, as long as the dynamics $f$ are differentiable around $(\st^\star, \act^\star)$ and the cost function $c$ is twice differentiable at $(\st^\star, \act^\star)$, we can take a linear approximation of $f$ and a quadratic approximation of $c$ to bring us back to the regime of LQR.

Linearizing the dynamics around $(\st^\star, \act^\star)$ gives:

and quadratizing the cost function around $(\st^\star, \act^\star)$ gives:

where the gradients and Hessians are defined as

Exercise: Note that this cost can be expressed in the general quadratic form seen in (2.38). Derive the corresponding quantities $Q, R, M, q, r, c$.

2.6.2Finite differencing¶

To calculate these gradients and Hessians in practice, we use a method known as finite differencing for numerically computing derivatives. Namely, we can simply use the limit definition of the derivative, and see how the function changes as we add or subtract a tiny δ to the input.

Note that this only requires us to be able to query the function, not to have an analytical expression for it, which is why it’s so useful in practice.

2.6.3Local convexification¶

However, simply taking the second-order approximation of the cost function is insufficient, since for the LQR setup we required that the $Q$ and $R$ matrices were positive definite, i.e. that all of their eigenvalues were positive.

One way to naively force some symmetric matrix $D$ to be positive definite is to set any non-positive eigenvalues to some small positive value $\varepsilon > 0$. Recall that any real symmetric matrix $D \in \mathbb{R}^{n \times n}$ has an basis of eigenvectors $u_1, \dots, u_n$ with corresponding eigenvalues $\lambda_1, \dots, \lambda_n$ such that $D u_i = \lambda_i u_i$. Then we can construct the positive definite approximation by

Exercise: Convince yourself that $\widetilde{D}$ is indeed positive definite.

Note that Hessian matrices are generally symmetric, so we can apply this process to $Q$ and $R$ to obtain the positive definite approximations $\widetilde{Q}$ and $\widetilde{R}$. Now that we have an upward-curved quadratic approximation to the cost function, and a linear approximation to the state transitions, we can simply apply the time-homogenous LQR methods from Section 2.4.

But what happens when we enter states far away from $\st^\star$ or want to use actions far from $\act^\star$? A Taylor approximation is only accurate in a local region around the point of linearization, so the performance of our LQR controller will degrade as we move further away. We’ll see how to address this in the next section using the iterative LQR algorithm.

2.6.4Iterative LQR¶

To address these issues with local linearization, we’ll use an iterative approach, where we repeatedly linearize around different points to create a time-dependent approximation of the dynamics, and then solve the resulting time-dependent LQR problem to obtain a better policy. This is known as iterative LQR or iLQR:

Now let’s go through the details of each step. We’ll use superscripts to denote the iteration of the algorithm. We’ll also denote $\bar \st_0 = \E_{\st_0 \sim \mu_0} [\st_0]$ as the expected initial state.

At iteration $i$ of the algorithm, we begin with a candidate trajectory $\bar \tau^i = (\bar \st^i_0, \bar \act^i_0, \dots, \bar \st^i_{\hor-1}, \bar \act^i_{\hor-1})$.

Step 1: Form a time-dependent LQR problem. At each timestep $\hi \in [\hor]$, we use the techniques from Section 2.6 to linearize the dynamics and quadratize the cost function around $(\bar \st^i_\hi, \bar \act^i_\hi)$:

Step 2: Compute the optimal policy. We can now solve the time-dependent LQR problem using the Riccati equation from Section 2.5.1 to compute the optimal policy $\pi^i_0, \dots, \pi^i_{\hor-1}$.

Step 3: Generate a new series of actions. We can then generate a new sample trajectory by taking actions according to this optimal policy:

Note that the states are sampled according to the true dynamics, which we assume we have query access to.

Step 4: Compute a better candidate trajectory., Note that we’ve denoted these actions as $\widetilde \act_\hi$ and aren’t directly using them for the next iteration $\bar \act^{i+1}_\hi$. Rather, we want to interpolate between them and the actions from the previous iteration $\bar \act^i_0, \dots, \bar \act^i_{\hor-1}$. This is so that the cost will increase monotonically, since if the new policy turns out to actually be worse, we can stay closer to the previous trajectory. (Can you think of an intuitive example where this might happen?)

Formally, we want to find $\alpha \in [0, 1]$ to generate the next iteration of actions $\bar \act^{i+1}_0, \dots, \bar \act^{i+1}_{\hor-1}$ such that the cost is minimized:

Note that this optimizes over the closed interval $[0, 1]$, so by the Extreme Value Theorem, it’s guaranteed to have a global maximum.

The final output of this algorithm is a policy $\pi^{n_\text{steps}}$ derived after $n_\text{steps}$ of the algorithm. Though the proof is somewhat complex, one can show that for many nonlinear control problems, this solution converges to a locally optimal solution (in the policy space).

2.7Summary¶

This chapter introduced some approaches to solving different variants of the optimal control problem Definition 2.1. We began with the simple case of linear dynamics and an upward-curved quadratic cost. This model is called the LQR and we solved for the optimal policy using dynamic programming. We then extended these results to the more general nonlinear case via local linearization. We finally saw the iterative LQR algorithm for solving nonlinear control problems.