M4-T4 — The Linear-Quadratic Regulator from Dynamic Programming

Learning objectives

After this lesson, you should be able to:

• Derive the Riccati equation and the optimal gain from dynamic programming with a quadratic value function.
• Choose the weights of a regulator from the largest acceptable value of each state and of the input (Bryson's rule).
• Prove that an LQR loop has |1 + L| ≥ 1, and deduce its gain and phase margins.
• Say when the guarantees hold and when they do not.

The intuition

UAV M5-T7: LQR control introduced the regulator for a quadrotor hover, with the weights found by trial, and Navigation M5-T6 used it to track a path. Both treat the Riccati equation as a tool. This lesson looks inside: where the equation comes from, how to set the weights from physical limits instead of by trial, and why the resulting loop is unusually robust.

What it optimizes

A regulator brings the state to zero. The faster it does so, the larger the actuator effort. LQR writes the compromise as a cost, accumulated over the whole motion:

$J = \int_0^\infty \bigl(\mathbf{x}^{\mathsf{T}}\mathbf{Q}\mathbf{x} + u^{\mathsf{T}}R\,u\bigr)\,dt$. The first term charges for every state error, the second for every unit of control effort. The two matrices are the designer's statement of what matters. The remarkable fact is that the minimizing input is a linear feedback of the state, $u = -\mathbf{K}\mathbf{x}$, whose gain comes from one matrix equation.

Where the equation comes from

The argument is dynamic programming (the same principle that the Robot Learning track uses for value iteration). Define the value of a state as the smallest cost that can still be incurred starting from it. The principle of optimality says that whatever the first moment's input, the rest of the motion must be optimal from where it lands. Written for a small time step, that gives an equation in the value alone. For a linear model and a quadratic cost the value is a quadratic $V = \mathbf{x}^{\mathsf{T}}\mathbf{P}\mathbf{x}$, and the equation for $\mathbf{P}$ is the algebraic Riccati equation.

Weights from limits

The weights are scaled by Bryson's rule: take the reciprocal of the square of the largest acceptable value. If the position may be off by 1 m, the weight on it is $1/1^2$; if the gimbal may swing 0.15 rad, the weight on the input is $1/0.15^2 = 44.4$. Then every term of the cost equals 1 when its quantity is at its limit, and the weights say, in common units, how much each limit matters. For Hopper: position 1 m, velocity 1 m/s, pitch 0.2 rad, pitch rate 0.5 rad/s, gimbal 0.15 rad. The gain comes out as $\mathbf{K} = (0.150,\ 0.272,\ 1.658,\ 0.386)$ with closed-loop poles at −5.87, −1.61 ± 1.25j and −1.16, found by no one placing anything.

The gift: margins

The loop gain of an LQR design, cut at the plant input, is $L(s) = \mathbf{K}(s\mathbf{I} - \mathbf{A})^{-1}\mathbf{B}$. The Riccati equation forces it to satisfy $|1 + L(j\omega)| \ge 1$ at every frequency. That inequality has a picture: the Nyquist curve never enters the disk of radius 1 centred on the critical point −1. Everything else follows from the geometry: the gain can be raised without limit and lowered to half, and the phase margin is at least 60°. No such guarantee comes with a design by pole placement.

LQR starts from what you care about and returns a gain whose loop cannot have small margins, provided the model it was designed on is the loop that is closed. The guarantee has conditions, and the last section of the lesson is about them.

Definitions

Cost $J$

The integral over time of the weighted state error and weighted input: $J = \int_0^\infty (\mathbf{x}^{\mathsf{T}}\mathbf{Q}\mathbf{x} + u^{\mathsf{T}}\mathbf{R}u)\,dt$.

Weights $\mathbf{Q}$ and $\mathbf{R}$

$\mathbf{Q} \succeq 0$ (positive semidefinite) and $\mathbf{R} \succ 0$ (positive definite): the price of each state error and of the input.

Value function $V(\mathbf{x})$

The smallest cost still to be paid starting from state $\mathbf{x}$; for LQR, $V = \mathbf{x}^{\mathsf{T}}\mathbf{P}\mathbf{x}$.

Return difference $1 + L$

With $L = \mathbf{K}(s\mathbf{I} - \mathbf{A})^{-1}\mathbf{B}$ the loop gain cut at the plant input. Its magnitude is the distance from the Nyquist curve to −1, the inverse of the sensitivity.

For a stabilizable model with the weights above, the optimal input and its gain are:

Bryson's rule for the weights, with the largest acceptable value of each quantity written $x_{i,\max}$ and $u_{\max}$:

The guarantee of the margins comes from an identity, the return-difference equality of Kalman, for a single input:

Derivation: the Riccati equation, and why |1 + L| ≥ 1

The derivation has two parts: the optimal gain from dynamic programming, then the margin identity from the Riccati equation.

1. $$V(\mathbf{x}) = \min_u\left[\,\mathbf{x}^{\mathsf{T}}\mathbf{Q}\mathbf{x}\,dt + u^{\mathsf{T}}\mathbf{R}u\,dt + V(\mathbf{x} + \dot{\mathbf{x}}\,dt)\,\right]$$

   The principle of optimality: the cost from $\mathbf{x}$ is the cost of the next instant plus the optimal cost from where the state then is.

2. $$0 = \min_u\left[\mathbf{x}^{\mathsf{T}}\mathbf{Q}\mathbf{x} + u^{\mathsf{T}}\mathbf{R}u + \nabla V^{\mathsf{T}}(\mathbf{A}\mathbf{x} + \mathbf{B}u)\right]$$

   Expand $V(\mathbf{x} + \dot{\mathbf{x}}dt)$ to first order in $dt$ and divide by $dt$. In steady state the value does not change with time, which gives the Hamilton-Jacobi-Bellman equation.

3. $$V = \mathbf{x}^{\mathsf{T}}\mathbf{P}\mathbf{x} \;\Rightarrow\; \nabla V = 2\mathbf{P}\mathbf{x}$$

   Guess that the value is a quadratic form. A linear model with a quadratic cost makes this consistent, as the next steps confirm.

4. $$2\mathbf{R}u + 2\mathbf{B}^{\mathsf{T}}\mathbf{P}\mathbf{x} = 0 \;\Rightarrow\; u = -\mathbf{R}^{-1}\mathbf{B}^{\mathsf{T}}\mathbf{P}\mathbf{x}$$

   Set the derivative of the bracket with respect to $u$ to zero. The optimal input is linear in the state.

5. $$\mathbf{x}^{\mathsf{T}}\bigl(\mathbf{Q} + \mathbf{A}^{\mathsf{T}}\mathbf{P} + \mathbf{P}\mathbf{A} - \mathbf{P}\mathbf{B}\mathbf{R}^{-1}\mathbf{B}^{\mathsf{T}}\mathbf{P}\bigr)\mathbf{x} = 0$$

   Put the optimal input back into the bracket. For this to hold for every state, the matrix in brackets must vanish: the Riccati equation.

The margin identity

1. $$\mathbf{P}(j\omega\mathbf{I} - \mathbf{A}) + (-j\omega\mathbf{I} - \mathbf{A}^{\mathsf{T}})\mathbf{P} = -(\mathbf{A}^{\mathsf{T}}\mathbf{P} + \mathbf{P}\mathbf{A}) = \mathbf{Q} - \mathbf{P}\mathbf{B}\mathbf{R}^{-1}\mathbf{B}^{\mathsf{T}}\mathbf{P}$$

   The $j\omega\mathbf{P}$ terms cancel, and the Riccati equation gives the right side.

2. $$\mathbf{B}^{\mathsf{T}}(-j\omega\mathbf{I} - \mathbf{A}^{\mathsf{T}})^{-1}\bigl[\cdots\bigr](j\omega\mathbf{I} - \mathbf{A})^{-1}\mathbf{B}$$

   Multiply by $\mathbf{B}^{\mathsf{T}}(-j\omega\mathbf{I} - \mathbf{A}^{\mathsf{T}})^{-1}$ on the left and $(j\omega\mathbf{I} - \mathbf{A})^{-1}\mathbf{B}$ on the right. The terms with $\mathbf{P}$ become $\mathbf{B}^{\mathsf{T}}\mathbf{P}(j\omega\mathbf{I} - \mathbf{A})^{-1}\mathbf{B} = R\,L(j\omega)$, its conjugate, and $L^*RL$.

3. $$R + R\,L + L^{*}R + L^{*}R\,L = R + \mathbf{B}^{\mathsf{T}}(\cdots)^{-1}\mathbf{Q}(\cdots)^{-1}\mathbf{B}$$

   Collect and add $R$ to both sides. The left side factors as $R\,(1 + L)^{*}(1 + L)$.

4. $$|1 + L(j\omega)|^2 \ge 1$$

   The extra term on the right is $\mathbf{v}^{*}\mathbf{Q}\mathbf{v} \ge 0$ for a vector $\mathbf{v}$, so the right side is at least $R$.

The margins that follow

1. $$L \text{ real and negative} \;\Rightarrow\; |1 + L| \ge 1 \;\Rightarrow\; L \le -2$$

   On the negative real axis the curve cannot enter the interval (−2, 0): the points there are less than 1 from −1.

2. $$k\,L = -1 \text{ needs } k = -1/L \le 1/2 \text{ or no crossing at all}$$

   Scaling the gain by $k$ moves the curve radially. A crossing of −1 can happen only for $k \le 1/2$: the gain can be raised without limit and cut to one half.

3. $$|1 + e^{j(\mathrm{PM} - 180^\circ)}| = 2\sin\tfrac{\mathrm{PM}}{2} \ge 1 \;\Rightarrow\; \mathrm{PM} \ge 60^\circ$$

   At the crossover the curve is on the unit circle at an angle PM from −1, at distance $2\sin(\mathrm{PM}/2)$ from it. That must be at least 1.

The Riccati equation forces the return difference to be at least 1 everywhere, and geometry turns that into at least 60° of phase margin and a gain margin from 1/2 to infinity. The guarantee is for the loop cut at the plant input, with the model the weights were designed for.

A worked example

Design the regulator for Hopper's lateral channel and examine it.

1. Weights. Bryson's rule with limits of 1 m, 1 m/s, 0.2 rad and 0.5 rad/s for the states and 0.15 rad for the gimbal gives $\mathbf{Q} = \operatorname{diag}(1, 1, 25, 4)$ and $R = 44.4$.
2. Solve. The Riccati equation (the library solves it by the matrix sign function, checked against SciPy to 10⁻¹³) gives $\mathbf{K} = (0.150,\ 0.272,\ 1.658,\ 0.386)$ and closed-loop poles at −5.87, −1.61 ± 1.25j and −1.16. A 1 m position error alone calls for 0.15 rad of gimbal, exactly the limit: Bryson's rule did what it says.
3. The identity. The smallest value of $|1 + L(j\omega)|$ over all frequencies is 1.000018: at the bound, as the theorem says.
4. Margins. The loop crosses 0 dB at 10.22 rad/s with a phase margin of 71.3°, above the guaranteed 60°, and a delay margin of 0.122 s. It stays stable for any increase of gain and for a decrease down to 0.19 of the design gain, better than the guaranteed 0.5.
5. Add the servo. The design did not include the 40 rad/s gimbal servo. With it the crossover is 9.93 rad/s, the phase margin 56.8° (below 60°: the guarantee does not cover the servo), the delay margin 0.100 s and the peak of the sensitivity 1.22. The servo costs $\arctan(10.2/40) = 14^\circ$ at the crossover, and $71.3 - 14 = 57$ predicts the result.
6. Optimality. The cost of starting 1 m off is $P_{11} = 1.81$ for LQR (and the Lyapunov equation of the closed loop gives the same, 1.81). The pole placement of the last lesson with poles at (−1, −1.5, −3, −4) costs 1.90, 5% more, and the faster placement (−2, −3, −6, −8) costs 5.93, more than three times as much.

import cmath, math
import gnckit as g

gr, lam, b = 9.81, 0.7, 19.62
A = [[0, 1, 0, 0], [0, 0, gr, 0], [0, 0, 0, 1], [0, 0, lam ** 2, 0]]
B = [0, gr, 0, b]
x_max, v_max, th_max, q_max, d_max = 1.0, 1.0, 0.2, 0.5, 0.15       # Bryson's rule: the largest acceptable value of each quantity
Q = [[1 / x_max ** 2, 0, 0, 0], [0, 1 / v_max ** 2, 0, 0], [0, 0, 1 / th_max ** 2, 0], [0, 0, 0, 1 / q_max ** 2]]
R = [[1 / d_max ** 2]]
K, P = g.lqr(A, [[v] for v in B], Q, R)
K = K[0]
print('K =', [round(k, 4) for k in K])
poles = g.mat_eigvals(g.ss_closed_loop(A, B, K))
print('closed-loop poles:', sorted((round(complex(p).real, 2), round(complex(p).imag, 2)) for p in poles))

L = lambda w: g.ss_freq(A, B, K, 0.0, w)                                # loop broken at the plant input
ws = g.logspace(-2, 3, 4000)
print(f'smallest |1 + L| over all frequencies: {min(abs(1 + L(w)) for w in ws):.6f} (the theorem says at least 1)')
wc = [ws[i] for i in range(len(ws) - 1) if (abs(L(ws[i])) - 1) * (abs(L(ws[i + 1])) - 1) < 0][0]
pm = 180 + math.degrees(cmath.phase(L(wc)))
print(f'ideal LQR loop: crossover {wc:.2f} rad/s, phase margin {pm:.1f} deg (at least 60), delay margin {math.radians(pm) / wc:.3f} s')
lo = min(k for k in [10 ** (i / 100) for i in range(-300, 100)] if max(complex(p).real for p in g.mat_eigvals(g.ss_closed_loop(A, B, [k * c for c in K]))) < 0)
print(f'it stays stable for any gain increase and for a gain decrease down to x{lo:.2f} (the theorem guarantees x0.5)')

Ls = lambda w: L(w) / (1 + 1j * w / 40)                                  # the same loop with the 40 rad/s gimbal servo
wc = [ws[i] for i in range(len(ws) - 1) if (abs(Ls(ws[i])) - 1) * (abs(Ls(ws[i + 1])) - 1) < 0][0]
pm = 180 + math.degrees(cmath.phase(Ls(wc)))
print(f'with the servo: crossover {wc:.2f} rad/s, phase margin {pm:.1f} deg, delay margin {math.radians(pm) / wc:.3f} s, peak |S| {max(1 / abs(1 + Ls(w)) for w in ws):.2f}')

# optimality: cost of a 1 m position error, J = x0' P x0, for LQR and for two pole-placement designs
def cost(Kv):
    Acl = g.ss_closed_loop(A, B, Kv)
    Qc = [[Q[i][j] + Kv[i] * R[0][0] * Kv[j] for j in range(4)] for i in range(4)]
    return g.lyap(Acl, Qc)[0][0]
print(f'cost of a 1 m error: LQR {P[0][0]:.2f} (check {cost(K):.2f}); poles (-1,-1.5,-3,-4): {cost(g.acker(A, B, [-1, -1.5, -3, -4])):.2f}; poles (-2,-3,-6,-8): {cost(g.acker(A, B, [-2, -3, -6, -8])):.2f}')

K = [0.15, 0.2715, 1.6581, 0.3861]
closed-loop poles: [(-5.87, 0.0), (-1.61, -1.25), (-1.61, 1.25), (-1.16, 0.0)]
smallest |1 + L| over all frequencies: 1.000018 (the theorem says at least 1)
ideal LQR loop: crossover 10.22 rad/s, phase margin 71.3 deg (at least 60), delay margin 0.122 s
it stays stable for any gain increase and for a gain decrease down to x0.19 (the theorem guarantees x0.5)
with the servo: crossover 9.93 rad/s, phase margin 56.8 deg, delay margin 0.100 s, peak |S| 1.22
cost of a 1 m error: LQR 1.81 (check 1.81); poles (-1,-1.5,-3,-4): 1.90; poles (-2,-3,-6,-8): 5.93

In robotics

LQR is the standard first design for multivariable regulation, and its margins are why it is trusted.

• Autopilots. Attitude and velocity regulators of aircraft, missiles and spacecraft are designed by LQR as often as by classical loop shaping. The weights express the requirements in the units of the problem, and an engineer can read the rule behind each gain.
• Weights from requirements. Bryson's rule connects a requirement ("keep the pitch within 0.2 rad, the gimbal within 0.15 rad") to a design directly. Tuning then moves the weights by looking at the closed loop, as the lab M4-L1: LQR for Hopper and Its Margins does.
• Checking a design. The identity $|1 + L| \ge 1$ is a diagnostic: an LQR loop that violates it has a bug, an unmodelled actuator or a wrong weight. The worked example finds the 56.8° of the servo loop by this reasoning.
• Beyond the regulator. The same machinery gives the Kalman gain of M4-T3: Observers and the Separation Principle (the dual problem), the linear-quadratic-Gaussian controller of the next lesson, and, over a finite horizon with constraints, the model predictive control of module 9.
• Limits. The solution needs the whole state (the observer of the last lesson) and a model that is good in the region used. It does not know about actuator saturation, so the weight on the input has to be large enough to keep the input within its limits for the errors expected, as Bryson's rule does.

For the quadrotor version see UAV M5-L2: Implement LQR; this track's labs for the same ideas are in module 4. The math of positive-definite matrices that the cost relies on is in Linear Algebra M6-T4: Positive Definite Matrices.

Common misconceptions

Key takeaways

1. LQR minimizes $\int(\mathbf{x}^{\mathsf{T}}\mathbf{Q}\mathbf{x} + u^{\mathsf{T}}\mathbf{R}u)\,dt$ by linear state feedback $u = -\mathbf{R}^{-1}\mathbf{B}^{\mathsf{T}}\mathbf{P}\mathbf{x}$, where $\mathbf{P}$ solves the Riccati equation.
2. The equation comes from dynamic programming with a quadratic value function $V = \mathbf{x}^{\mathsf{T}}\mathbf{P}\mathbf{x}$; the cost of a start $\mathbf{x}_0$ is $\mathbf{x}_0^{\mathsf{T}}\mathbf{P}\mathbf{x}_0$.
3. Bryson's rule sets each weight to one over the square of the largest acceptable value; only the ratios of the weights matter.
4. The Riccati equation implies $|1 + L(j\omega)| \ge 1$, hence a gain margin from 1/2 to infinity and a phase margin of at least 60°.
5. The guarantee holds at the plant input, with full state and the design model; Hopper's servo takes the margin from 71.3° to 56.8°.
6. For Hopper, K = (0.150, 0.272, 1.658, 0.386), with the cost of a 1 m error 1.81, below any pole placement tried.

References and further reading

Each source is marked free or paid. Free to read does not mean free to reuse: the course resources page lists every source with its licence.

• K. J. Åström and R. M. Murray, Feedback Systems, 2nd edition, Princeton University Press, 2021 , the chapter on optimal control (the linear-quadratic regulator). Free to read as the authors' PDF; the book is copyrighted.
• MIT 16.323 Principles of Optimal Control (OpenCourseWare) , the lectures on dynamic programming and the LQR. Free; licence CC BY-NC-SA 4.0.
• Z. Manchester, CMU 16-745 Optimal Control and Reinforcement Learning , with notes and notebooks on the LQR and its derivation. Free to read; no licence statement was found.
• B. D. O. Anderson and J. B. Moore, Optimal Control: Linear Quadratic Methods, Prentice Hall, 1989, the chapters on the regulator and its robustness. The classic text (a scan is on an author's page, with no licence statement). Free to read there.

Knowledge check

Three questions. Write an answer, then reveal the explanation. Your answers stay in your browser.

What’s next

Related content

Nothing below is required to finish this lesson.

• Before
  M4-T3: Observers and the Separation Principle

  Design a Luenberger observer, show that controller and observer poles can be chosen separately, and identify the Kalman filter as the optimal observer.

• Lab
  M4-L1: LQR for Hopper and Its Margins

  Design an LQR for Hopper's lateral channel (position, velocity, pitch, pitch rate), then break the loop at the gimbal and measure the margins it actually has with the actuator lag included.

• Lab
  M12-L2: Capstone: Land Hopper

  Write the navigation filter, the guidance law and the controller of Hopper, and pass an acceptance campaign of dispersed flights with wind, sensor errors, mass uncertainty and one beacon fault.

Other courses on this site treat neighbouring topics from their own angle:

• UAV M5-T7: LQR control
• UAV M5-L2: Implement LQR
• Navigation M5-T6: LQR and MPC for Path Tracking
• UAV M3-T4: Kalman Filter
• UAV M3-T12: Observability and consistency
• Linear Algebra M6-T1: Eigenvalues and Eigenvectors
• Linear Algebra M6-T4: Positive Definite Matrices
• Linear Algebra M9-T7: Linear Kalman Filter Mathematics
• Robot Learning M4-T4: Dynamic Programming and Value Iteration
• Probability for Robotics (module 7 derives the Kalman filter)

This course leaves the following to the pages that already teach them:

• The derivation of the Kalman filter
• LQR applied to a linearized quadrotor
• LQR applied to path tracking of a ground robot