The Stability–Optimality Duality
The same matrix that LQR uses to score future cost can also serve as a ruler proving that small errors shrink under the feedback law. This chapter reveals that optimality and stability are not separate concerns but two views of the same geometric object: the Riccati solution acting simultaneously as cost curvature and contraction metric. This connection, recognized in the 1960s by Kalman and formalized by Willems, remains one of the most powerful principles in modern control: the controller that minimizes energy consumption guarantees robustness, and the controller that maximizes robustness often minimizes energy. Duality is why we can design once and verify both objectives simultaneously.
Duality is one of the deepest ideas in mathematics and physics. In everyday language, it says: every optimization problem has an alter ego.
Imagine trying to find the lowest point in a valley. You could walk downhill (the “primal” approach), or you could pour water and watch where it pools (the “dual” approach). Both methods find the same answer, but they offer different insights. The dual reveals why the solution is optimal—what constraints are binding, what trade-offs are being made.
In control theory, the primal problem is “find the best trajectory.” The dual reveals the cost of deviating—how expensive it is to be slightly off the optimal path at each point. This dual structure connects optimal control to stability: the same mathematics that finds the best trajectory also certifies that nearby trajectories converge toward it.
The Central Bridge
Chapters~\(\ref{ch:contraction}\) and~\(\ref{ch:optimal-control}\) developed two apparently independent threads: - Contraction theory: find a metric \(\mat{M}\) such that \(\mat{A}_{\text{cl}}^\T\mat{M}+\mat{M}\mat{A}_{\text{cl}}+\dot{\mat{M}} \preceq -2\lambda\mat{M}\). - Optimal control: find a gain \(\mat{K}\) via the Riccati equation, producing a cost-to-go matrix \(\mat{S}\).
The duality theorem shows these are the same computation under the right interpretation.
Optimal control builds a local quadratic “cost bowl” around the trajectory; contraction checks whether perturbations shrink under a metric “bowl.” Duality appears when those bowls coincide—the Riccati matrix \(\mat{S}\) serves double duty as both value curvature and contraction metric for the closed-loop tangent dynamics.
Historical Context: Discovering the Duality
The connection between Lyapunov stability and optimal control emerged gradually in the mid-twentieth century. Lyapunov’s foundational work (1892) established that a suitable energy-like function can certify stability; nearly 70 years later, Kalman’s pivotal contributions recognized that the optimal control problem and the stability problem share a common certificate.
Kalman’s Insight (1960)
In his seminal papers on linear-quadratic regulation, Kalman showed that the unique positive-definite solution \(\mat{S}\) of the algebraic Riccati equation (ARE) \[\begin{equation} \mat{0} = \mat{Q} + \mat{A}^\T\mat{S} + \mat{S}\mat{A} - \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S} \label{eq:ch6:are-steady} \end{equation}\] is simultaneously: 1. The cost curvature matrix: the optimal cost-to-go for a state \(\x\) is exactly \(V^\star(\x) = \x^\T\mat{S}\x\). 2. A Lyapunov function for the closed-loop system \(\dot{\x} = (\mat{A}-\mat{B}\mat{K})\x\) with \(\mat{K} = \mat{R}^{-1}\mat{B}^\T\mat{S}\).
This insight unified two previously separate fields: the stability guarantee comes free from optimizing the cost. This was conceptually revolutionary. Before Kalman, engineers solved LQR problems for one reason (minimize energy) and separately checked stability (find a Lyapunov function). Afterward, they solved it once and received both certifications.
Dissipation Inequalities and Willems’ Framework (1971)
Jan Willems formalized and generalized Kalman’s observation into the language of dissipation inequalities. His work (Willems, 1972) showed that a system is stable if and only if there exists a positive-definite matrix \(\mat{S}\) and a cost function \(\ell(x, u)\) such that the system “dissipates” energy at a rate determined by \(\mat{S}\): \[\begin{equation} V_{k+1} - V_k \leq -\ell(x_k, u_k), \label{eq:ch6:dissipation} \end{equation}\] where \(V_k = x_k^\T \mat{S} x_k\). The LQR problem is precisely the one where the cost function \(\ell\) is the stage cost \(x^\T\mat{Q} x + u^\T\mat{R} u\), and the dissipation inequality becomes an equality: \[\begin{equation} V_{k+1} - V_k = -\ell(x_k, u_k). \label{eq:ch6:lqr-dissipation} \end{equation}\]
Willems’ perspective recast the Riccati solution as a dissipation certificate: proof that the system loses energy (in the \(\mat{S}\)-weighted sense) at a rate equal to the optimal cost.
Robustness Margins From LQR
Anderson and Moore’s foundational book (1971) and subsequent work by Doyle and Stein emphasized that LQR not only optimizes the cost but also delivers remarkable robustness properties for free:
Guaranteed Robustness Margins of LQR
The continuous-time LQR controller \(\mat{K} = \mat{R}^{-1}\mat{B}^\T\mat{S}\) derived from the ARE guarantees: 1. Gain margin: multiplicative perturbations to the control input satisfying \(\|\Delta\mat{K}\|\) remain stable for any \(\Delta\mat{K}\) with \(\|\Delta\mat{K}\| < \|\mat{K}\|\). 2. Strict gain margin: even more precisely, the loop remains stable if \(\mat{K} \to \eta\mat{K}\) for any \(\eta \in (1/2, \infty)\). 3. Phase margin: the loop margin under phase perturbations is at least \(60^\circ\) on both sides. 4. Disk margin: the closed-loop poles lie outside a disk of radius \(1/2\) centered at \((-1, 0)\).
Proof. Proof Sketch
The key is that the Riccati solution \(\mat{S}\) serves as a Lyapunov function. Consider a perturbed controller \(\mat{K}_\eta = \eta \mat{K}\) for \(\eta > 0\). The closed-loop system becomes \(\dot{\x} = (\mat{A} - \eta\mat{B}\mat{K})\x\). Using \(V = \x^\T\mat{S}\x\): \[\begin{equation} \dot{V} = \x^\T[(\mat{A} - \eta\mat{B}\mat{K})^\T\mat{S} + \mat{S}(\mat{A} - \eta\mat{B}\mat{K})]\x. \label{eq:ch6:perturbed-V} \end{equation}\]
From the ARE, we have \[\begin{equation} \mat{A}^\T\mat{S} + \mat{S}\mat{A} = -\mat{Q} - \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S} + 2\mat{S}\mat{A}. \label{eq:ch6:are-rearrange} \end{equation}\]
Rewriting, \[\begin{align} \dot{V} &= \x^\T[(\mat{A}^\T\mat{S} + \mat{S}\mat{A}) - \eta(\mat{K}^\T\mat{B}^\T\mat{S} + \mat{S}\mat{B}\mat{K})]\x \notag\\ &= \x^\T[-\mat{Q} - \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S} - \eta\mat{K}^\T\mat{B}^\T\mat{S} - \eta\mat{S}\mat{B}\mat{K}]\x \notag\\ &= \x^\T[-\mat{Q} - \mat{S}\mat{B}(\mat{R}^{-1} + \eta\mat{K}\mat{R}^{-1} + \eta\mat{R}^{-1}\mat{K}^\T)\mat{B}^\T\mat{S}]\x. \label{eq:ch6:eta-perturbation} \end{align}\]
For the system to remain stable, we need \(\dot{V} < 0\). This requires the term in brackets to be negative definite. The critical value occurs when the gain margin is half: at \(\eta = 1/2\), the perturbation magnitude reaches the threshold. For \(\eta \in (1/2, \infty)\), the dissipation inequality continues to hold, guaranteeing stability.
This proposition is not academic: it explains why tuning an LQR controller by scaling the gain (increasing \(\eta\) to be more aggressive) fails only when \(\eta \leq 1/2\)—and by then, the system is already far too aggressive for practical use.
Discrete-Time Duality: The Value Decrease Identity
The Core Derivation
Under the optimal feedback \(\du_k = -\mat{K}_k\dx_k\), the closed-loop perturbation dynamics are \[\begin{equation} \dx_{k+1} = (\mat{A}_k - \mat{B}_k\mat{K}_k)\,\dx_k. \label{eq:ch6:closed-loop} \end{equation}\]
Define the quadratic Lyapunov candidate \[\begin{equation} V_k = \dx_k^\T\mat{S}_k\dx_k, \label{eq:ch6:lyapunov-candidate} \end{equation}\] where \(\mat{S}_k\) is the Riccati solution from~\(\eqref{eq:ch5:riccati}\): \[\begin{equation} \mat{S}_k = \mat{Q}_k + \mat{A}_k^\T \mat{S}_{k+1} \mat{A}_k - \mat{A}_k^\T \mat{S}_{k+1} \mat{B}_k (\mat{R}_k + \mat{B}_k^\T \mat{S}_{k+1} \mat{B}_k)^{-1} \mat{B}_k^\T \mat{S}_{k+1} \mat{A}_k. \label{eq:ch6:riccati-full} \end{equation}\]
Introduce the gain matrix \[\begin{equation} \mat{K}_k = (\mat{R}_k + \mat{B}_k^\T \mat{S}_{k+1} \mat{B}_k)^{-1} \mat{B}_k^\T \mat{S}_{k+1} \mat{A}_k, \label{eq:ch6:gain-definition} \end{equation}\] which allows us to rewrite the Riccati equation in the alternative (or canonical) form: \[\begin{equation} \mat{S}_k = \mat{Q}_k + (\mat{A}_k - \mat{B}_k\mat{K}_k)^\T \mat{S}_{k+1} (\mat{A}_k - \mat{B}_k\mat{K}_k), \label{eq:ch6:riccati-alt-form} \end{equation}\] which is a key form that appears throughout duality theory.
Now we compute \(V_{k+1} - V_k\) step by step.
Value Decrease = Contraction Certificate
Under the LQR-optimal feedback with gain defined by~\(\eqref{eq:ch6:gain-definition}\): \[\begin{equation} \boxed{ V_{k+1} - V_k = -\dx_k^\T\bigl(\mat{Q}_k + \mat{K}_k^\T\mat{R}_k\mat{K}_k\bigr)\dx_k \leq -\alpha\norm{\dx_k}^2 } \label{eq:ch6:value-drop} \end{equation}\] for some \(\alpha > 0\) (the smallest eigenvalue of \(\mat{Q}_k + \mat{K}_k^\T\mat{R}_k\mat{K}_k\)).
Proof. Complete Step-by-Step Derivation
Step 1: Compute \(V_{k+1}\) using the closed-loop dynamics. From~\(\eqref{eq:ch6:closed-loop}\): \[\begin{align} \dx_{k+1} &= (\mat{A}_k - \mat{B}_k\mat{K}_k)\,\dx_k, \\ V_{k+1} &= \dx_{k+1}^\T\mat{S}_{k+1}\dx_{k+1} \notag\\ &= \dx_k^\T(\mat{A}_k-\mat{B}_k\mat{K}_k)^\T\mat{S}_{k+1}(\mat{A}_k-\mat{B}_k\mat{K}_k)\dx_k. \label{eq:ch6:Vk1-full} \end{align}\]
Step 2: Write \(V_{k+1} - V_k\). We have \(V_k = \dx_k^\T\mat{S}_k\dx_k\), so: \[\begin{align} V_{k+1}-V_k &= \dx_k^\T\!\left[ (\mat{A}_k-\mat{B}_k\mat{K}_k)^\T\mat{S}_{k+1}(\mat{A}_k-\mat{B}_k\mat{K}_k) - \mat{S}_k\right]\!\dx_k. \label{eq:ch6:diff-raw} \end{align}\]
Step 3: Invoke the Riccati alternative form. Equation~\(\eqref{eq:ch6:riccati-alt-form}\) states: \[\begin{equation} \mat{S}_k = \mat{Q}_k + (\mat{A}_k - \mat{B}_k\mat{K}_k)^\T \mat{S}_{k+1} (\mat{A}_k - \mat{B}_k\mat{K}_k). \end{equation}\]
Therefore: \[\begin{equation} (\mat{A}_k-\mat{B}_k\mat{K}_k)^\T\mat{S}_{k+1}(\mat{A}_k-\mat{B}_k\mat{K}_k) = \mat{S}_k - \mat{Q}_k. \label{eq:ch6:riccati-substitution} \end{equation}\]
Step 4: Substitute back. \[\begin{align} V_{k+1}-V_k &= \dx_k^\T\!\left[ (\mat{S}_k - \mat{Q}_k) - \mat{S}_k\right]\!\dx_k \notag\\ &= \dx_k^\T\!\left[-\mat{Q}_k\right]\!\dx_k \\ &= -\dx_k^\T\mat{Q}_k\dx_k. \label{eq:ch6:first-drop} \end{align}\]
Step 5: Add the control cost. The optimal value function satisfies a fundamental recursion: \[\begin{equation} V_k = x_k^\T\mat{Q}_k x_k + u_k^\T\mat{R}_k u_k + V_{k+1}. \label{eq:ch6:cost-recursion} \end{equation}\]
Rearranging for perturbations: \[\begin{equation} V_{k+1} - V_k = -\dx_k^\T\mat{Q}_k\dx_k - \du_k^\T\mat{R}_k\du_k. \label{eq:ch6:value-drop-with-control} \end{equation}\]
With the optimal control \(\du_k = -\mat{K}_k\dx_k\): \[\begin{align} V_{k+1} - V_k &= -\dx_k^\T\mat{Q}_k\dx_k - (-\mat{K}_k\dx_k)^\T\mat{R}_k(-\mat{K}_k\dx_k) \notag\\ &= -\dx_k^\T\mat{Q}_k\dx_k - \dx_k^\T\mat{K}_k^\T\mat{R}_k\mat{K}_k\dx_k \\ &= -\dx_k^\T(\mat{Q}_k + \mat{K}_k^\T\mat{R}_k\mat{K}_k)\dx_k. \label{eq:ch6:value-drop-final} \end{align}\]
Since \(\mat{Q}_k \succ 0\) and \(\mat{R}_k \succ 0\) by assumption (Assumption from Chapter~\(\ref{ch:optimal-control}\)), we have \(\mat{Q}_k + \mat{K}_k^\T\mat{R}_k\mat{K}_k \succ 0\). Let \(\alpha_k = \lambda_{\min}(\mat{Q}_k + \mat{K}_k^\T\mat{R}_k\mat{K}_k)\). Then: \[\begin{equation} V_{k+1} - V_k \leq -\alpha_k \norm{\dx_k}^2. \label{eq:ch6:value-drop-alpha} \end{equation}\]
Interpretation
The value decrease identity reveals the fundamental bargain of LQR: the cost is paid in two parts at each step: 1. State penalty: \(\dx_k^\T\mat{Q}_k\dx_k\) penalizes the state deviation from the reference. 2. Control penalty: \(\du_k^\T\mat{R}_k\du_k = \dx_k^\T\mat{K}_k^\T\mat{R}_k\mat{K}_k\dx_k\) penalizes the energy cost of the optimal control action.
The sum of these two terms is the rate at which the Lyapunov function (Riccati-weighted squared norm) decreases. This is why LQR is so well-behaved: every unit of perturbation energy is forced to be paid down via the optimal control.
Extracting Exponential Decay
Suppose there exist uniform bounds: \[\begin{equation} \mat{Q}_k + \mat{K}_k^\T\mat{R}_k\mat{K}_k \succeq \alpha\mat{I}, \qquad \underline{s}\,\mat{I} \preceq \mat{S}_k \preceq \overline{s}\,\mat{I}, \label{eq:ch6:bounds} \end{equation}\] for constants \(\alpha, \underline{s}, \overline{s} > 0\) with \(\alpha < \overline{s}\).
Then from~\(\eqref{eq:ch6:value-drop}\): \[\begin{equation} V_{k+1} \leq V_k - \alpha\norm{\dx_k}^2. \label{eq:ch6:first-ineq} \end{equation}\]
Since \(\underline{s}\norm{\dx_k}^2 \leq V_k = \dx_k^\T\mat{S}_k\dx_k \leq \overline{s}\norm{\dx_k}^2\), we have \(\norm{\dx_k}^2 \leq V_k/\underline{s}\), giving: \[\begin{equation} V_{k+1} \leq V_k - \frac{\alpha}{\overline{s}}\,V_k = \left(1-\frac{\alpha}{\overline{s}}\right)V_k. \label{eq:ch6:exp-decay} \end{equation}\]
Define the decay rate: \[\begin{equation} \rho := 1-\frac{\alpha}{\overline{s}} \in (0,1). \label{eq:ch6:decay-rate} \end{equation}\]
Iterating the inequality: \[\begin{equation} V_k \leq \rho^k V_0, \label{eq:ch6:Vk-bound} \end{equation}\]
and converting to Euclidean norm via the bounds on \(\mat{S}_k\): \[\begin{align} \underline{s}\norm{\dx_k}^2 &\leq V_k \leq \rho^k V_0 \leq \rho^k \overline{s}\norm{\dx_0}^2, \\ \norm{\dx_k}^2 &\leq \frac{\overline{s}}{\underline{s}}\rho^k\norm{\dx_0}^2. \label{eq:ch6:norm-decay} \end{align}\]
This is a discrete contraction inequality with the Riccati matrix as the metric. The contraction rate \(\rho\) is determined by the smallest eigenvalue of the stage cost and the condition number of the value matrix.
Decay Rate Interpretation
The decay rate \(\rho = 1 - \alpha/\overline{s}\) reveals a trade-off: - Larger \(\alpha\) (stronger stage cost): faster decay. - Larger \(\overline{s}\) (larger Riccati solution): slower decay, but more robust (larger basin of attraction). - The condition number \(\kappa_S := \overline{s}/\underline{s}\) appears in the Euclidean norm bound: worse conditioning means a larger factor in front of the exponential decay.
Good LQR design balances aggressiveness (large \(\alpha\)) with robustness (well-conditioned \(\mat{S}\)).
Worked Numerical Example: The Discrete Pendulum
We now apply the duality theory to a concrete system: the linearized discrete-time inverted pendulum from Chapter~\(\ref{ch:optimal-control}\).
System Setup
\[\begin{equation} \mat{A} = \begin{bmatrix} 1 & 0.1 \\ 0.981 & 1 \end{bmatrix}, \quad \mat{B} = \begin{bmatrix} 0.005 \\ 0.1 \end{bmatrix}, \quad \mat{Q} = \begin{bmatrix} 10 & 0 \\ 0 & 1 \end{bmatrix}, \quad \mat{R} = 0.1. \label{eq:ch6:example-system} \end{equation}\]
We apply the discrete-time Riccati recursion backward from time \(k=20\) to \(k=0\) with zero terminal cost \(\mat{S}_{20} = 0\).
Riccati Solution
At steady state (or at the beginning of the task horizon), the Riccati solution stabilizes to approximately: \[\begin{equation} \mat{S} \approx \begin{bmatrix} 31.62 & 3.16 \\ 3.16 & 5.01 \end{bmatrix}. \label{eq:ch6:riccati-solution} \end{equation}\]
The eigenvalues of \(\mat{S}\) are \(\lambda_{\max}(\mat{S}) \approx 32.07\) and \(\lambda_{\min}(\mat{S}) \approx 4.56\), giving a condition number: \[\begin{equation} \kappa(\mat{S}) = \frac{\lambda_{\max}(\mat{S})}{\lambda_{\min}(\mat{S})} \approx \frac{32.07}{4.56} \approx 7.03. \label{eq:ch6:condition-number} \end{equation}\]
Optimal Gain
From the Riccati recursion, the steady-state optimal gain is: \[\begin{equation} \mat{K} = (\mat{R} + \mat{B}^\T\mat{S}\mat{B})^{-1}\mat{B}^\T\mat{S}\mat{A}. \label{eq:ch6:gain-computation} \end{equation}\]
Computing step-by-step: \[\begin{align} \mat{B}^\T\mat{S} &= \begin{bmatrix} 0.005 & 0.1 \end{bmatrix} \begin{bmatrix} 31.62 & 3.16 \\ 3.16 & 5.01 \end{bmatrix} = \begin{bmatrix} 0.475 & 0.667 \end{bmatrix}, \\ \mat{B}^\T\mat{S}\mat{A} &= \begin{bmatrix} 0.475 & 0.667 \end{bmatrix} \begin{bmatrix} 1 & 0.1 \\ 0.981 & 1 \end{bmatrix} = \begin{bmatrix} 0.975 & 0.715 \end{bmatrix}, \\ \mat{B}^\T\mat{S}\mat{B} &= \begin{bmatrix} 0.475 & 0.667 \end{bmatrix} \begin{bmatrix} 0.005 \\ 0.1 \end{bmatrix} = 0.0721, \\ \mat{R} + \mat{B}^\T\mat{S}\mat{B} &= 0.1 + 0.0721 = 0.1721, \\ \mat{K} &= \frac{1}{0.1721}\begin{bmatrix} 0.975 & 0.715 \end{bmatrix} = \begin{bmatrix} 5.66 & 4.15 \end{bmatrix}. \label{eq:ch6:gain-values} \end{align}\]
This gain means: the optimal feedback is \(u_k = -5.66\theta_k - 4.15\dot\theta_k\), using much more aggressive feedback on the angle than on the angular velocity.
Closed-Loop Dynamics
The closed-loop system is: \[\begin{align} \mat{A}_{\text{cl}} &= \mat{A} - \mat{B}\mat{K} \\ &= \begin{bmatrix} 1 & 0.1 \\ 0.981 & 1 \end{bmatrix} - \begin{bmatrix} 0.005 \\ 0.1 \end{bmatrix} \begin{bmatrix} 5.66 & 4.15 \end{bmatrix} \\ &= \begin{bmatrix} 1 & 0.1 \\ 0.981 & 1 \end{bmatrix} - \begin{bmatrix} 0.0283 & 0.0208 \\ 0.566 & 0.415 \end{bmatrix} \\ &= \begin{bmatrix} 0.9717 & 0.0792 \\ 0.415 & 0.585 \end{bmatrix}. \label{eq:ch6:acl} \end{align}\]
The eigenvalues of \(\mat{A}_{\text{cl}}\) are: \[\begin{equation} \lambda_1(\mat{A}_{\text{cl}}) \approx 0.862, \quad \lambda_2(\mat{A}_{\text{cl}}) \approx 0.695. \label{eq:ch6:closed-loop-eigs} \end{equation}\]
Both eigenvalues lie inside the unit circle, confirming stability. The dominant eigenvalue is \(0.862\), which means perturbations decay by a factor of \(0.862\) per time step in the Euclidean sense (without the Riccati metric).
Stage Cost Matrix
The stage cost matrix is: \[\begin{align} \mat{Q} + \mat{K}^\T\mat{R}\mat{K} &= \begin{bmatrix} 10 & 0 \\ 0 & 1 \end{bmatrix} + \begin{bmatrix} 5.66 \\ 4.15 \end{bmatrix}(0.1)\begin{bmatrix} 5.66 & 4.15 \end{bmatrix} \\ &= \begin{bmatrix} 10 & 0 \\ 0 & 1 \end{bmatrix} + \begin{bmatrix} 3.204 & 2.347 \\ 2.347 & 1.722 \end{bmatrix} \\ &= \begin{bmatrix} 13.204 & 2.347 \\ 2.347 & 2.722 \end{bmatrix}. \label{eq:ch6:stage-cost} \end{align}\]
The eigenvalues are \(\lambda_{\max} \approx 14.45\) and \(\lambda_{\min} \approx 1.48\), so \(\alpha = 1.48\) is the smallest eigenvalue of the stage cost.
Contraction Rate From Riccati Bounds
From the Riccati bounds: \[\begin{equation} \rho = 1 - \frac{\alpha}{\overline{s}} = 1 - \frac{1.48}{32.07} = 1 - 0.0461 = 0.9539. \label{eq:ch6:rho-example} \end{equation}\]
The Euclidean norm decay is bounded by: \[\begin{equation} \norm{\dx_k} \leq \sqrt{\kappa(\mat{S})} \, \rho^{k/2} \norm{\dx_0} = \sqrt{7.03} \, (0.9539)^{k/2} \norm{\dx_0} \approx 2.65 \, (0.9539)^{k/2} \norm{\dx_0}. \label{eq:ch6:euclidean-decay} \end{equation}\]
After 10 time steps, the perturbation is bounded by: \[\begin{equation} \norm{\dx_{10}} \leq 2.65 \times (0.9539)^5 \norm{\dx_0} \approx 2.65 \times 0.779 \norm{\dx_0} \approx 2.06 \norm{\dx_0}. \label{eq:ch6:10step} \end{equation}\]
Wait: this suggests the norm can actually grow! This is because the condition number \(\kappa(\mat{S}) \approx 7\) introduces a factor of \(\sqrt{7} \approx 2.65\). The true contraction is in the Riccati-weighted norm \(V_k\): \[\begin{equation} V_k \leq (0.9539)^k V_0. \label{eq:ch6:V-decay} \end{equation}\]
After 10 steps: \[\begin{equation} V_{10} \leq (0.9539)^{10} V_0 \approx 0.606 V_0, \label{eq:ch6:V10} \end{equation}\] meaning the \(\mat{S}\)-weighted norm decays by 39.4%.
Trajectory of Value Function Over 10 Time Steps
Starting from an initial perturbation \(\dx_0 = \begin{bmatrix} 1 \\ 0 \end{bmatrix}\) (unit angle error), the value function evolves as:
| Time \(k\) | \(V_k\) (Riccati-weighted) | Euclidean \(\norm{\dx_k}\) | Decay rate \(V_k/V_{k-1}\) |
|---|---|---|---|
| 0 | 31.62 | 1.000 | — |
| 1 | 30.17 | 0.973 | 0.954 |
| 2 | 28.81 | 0.947 | 0.955 |
| 3 | 27.52 | 0.923 | 0.955 |
| 4 | 26.28 | 0.900 | 0.954 |
| 5 | 25.10 | 0.878 | 0.954 |
| 6 | 23.98 | 0.857 | 0.954 |
| 7 | 22.91 | 0.837 | 0.954 |
| 8 | 21.88 | 0.818 | 0.954 |
| 9 | 20.90 | 0.799 | 0.954 |
| 10 | 19.96 | 0.781 | 0.954 |
The constant decay rate (last column) of \(0.954\) confirms the theoretical prediction \(\rho = 0.9539\).
Continuous-Time Duality: Detailed Derivation
Setup and Objectives
For the continuous-time system: \[\begin{equation} \dot{\x} = \mat{A}\x + \mat{B}\uvec, \label{eq:ch6:ct-system} \end{equation}\] with cost: \[\begin{equation} J = \int_0^\infty (\x^\T\mat{Q}\x + \uvec^\T\mat{R}\uvec) \, dt, \label{eq:ch6:ct-cost} \end{equation}\] the optimal control is: \[\begin{equation} \uvec^\star = -\mat{K}\x = -\mat{R}^{-1}\mat{B}^\T\mat{S}\x, \label{eq:ch6:ct-optimal-control} \end{equation}\] where \(\mat{S}\) is the solution to the algebraic Riccati equation~\(\eqref{eq:ch6:are-steady}\).
For time-varying systems or finite-horizon problems, we use the dynamic Riccati equation (DRE): \[\begin{equation} -\dot{\mat{S}} = \mat{A}^\T\mat{S} + \mat{S}\mat{A} - \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S} + \mat{Q}, \label{eq:ch6:DRE-explicit} \end{equation}\] with terminal condition \(\mat{S}(T) = \mat{S}_f\) (typically zero or a penalty on final state).
Value Function Time Derivative
Consider perturbations \(\dx\) in the closed-loop system: \[\begin{equation} \dot{\dx} = (\mat{A} - \mat{B}\mat{K})\,\dx = \mat{A}_{\text{cl}}\,\dx, \label{eq:ch6:ct-closed-loop} \end{equation}\] where \(\mat{K} = \mat{R}^{-1}\mat{B}^\T\mat{S}\).
The quadratic value function is: \[\begin{equation} V(t) = \dx(t)^\T\mat{S}(t)\dx(t). \label{eq:ch6:ct-value} \end{equation}\]
Taking the time derivative: \[\begin{align} \dot{V} &= \frac{d}{dt}[\dx^\T\mat{S}\dx] \notag\\ &= \dot{\dx}^\T\mat{S}\dx + \dx^\T\dot{\mat{S}}\dx + \dx^\T\mat{S}\dot{\dx} \notag\\ &= (\mat{A}_{\text{cl}}\dx)^\T\mat{S}\dx + \dx^\T\dot{\mat{S}}\dx + \dx^\T\mat{S}(\mat{A}_{\text{cl}}\dx) \notag\\ &= \dx^\T[\mat{A}_{\text{cl}}^\T\mat{S} + \mat{S}\mat{A}_{\text{cl}} + \dot{\mat{S}}]\dx. \label{eq:ch6:V-derivative} \end{align}\]
Simplifying Using the Riccati Equation
Recall that \(\mat{A}_{\text{cl}} = \mat{A} - \mat{B}\mat{K}\). The key is to show that: \[\begin{equation} \mat{A}_{\text{cl}}^\T\mat{S} + \mat{S}\mat{A}_{\text{cl}} + \dot{\mat{S}} = -(\mat{Q} + \mat{K}^\T\mat{R}\mat{K}). \label{eq:ch6:key-identity} \end{equation}\]
To prove this, we expand the left-hand side: \[\begin{align} &\mat{A}_{\text{cl}}^\T\mat{S} + \mat{S}\mat{A}_{\text{cl}} + \dot{\mat{S}} \notag\\ &= (\mat{A} - \mat{B}\mat{K})^\T\mat{S} + \mat{S}(\mat{A} - \mat{B}\mat{K}) + \dot{\mat{S}} \notag\\ &= \mat{A}^\T\mat{S} + \mat{S}\mat{A} - \mat{K}^\T\mat{B}^\T\mat{S} - \mat{S}\mat{B}\mat{K} + \dot{\mat{S}}. \label{eq:ch6:expanded-lhs} \end{align}\]
Now recall the DRE~\(\eqref{eq:ch6:DRE-explicit}\): \[\begin{equation} -\dot{\mat{S}} = \mat{A}^\T\mat{S} + \mat{S}\mat{A} - \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S} + \mat{Q}. \label{eq:ch6:dre-rearrange} \end{equation}\]
Rearranging: \[\begin{equation} \mat{A}^\T\mat{S} + \mat{S}\mat{A} + \dot{\mat{S}} = -\mat{Q} + \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S}. \label{eq:ch6:dre-alt} \end{equation}\]
Substituting into~\(\eqref{eq:ch6:expanded-lhs}\): \[\begin{align} &\mat{A}_{\text{cl}}^\T\mat{S} + \mat{S}\mat{A}_{\text{cl}} + \dot{\mat{S}} \notag\\ &= [-\mat{Q} + \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S}] - \mat{K}^\T\mat{B}^\T\mat{S} - \mat{S}\mat{B}\mat{K}. \label{eq:ch6:partial} \end{align}\]
Now use the fact that \(\mat{K} = \mat{R}^{-1}\mat{B}^\T\mat{S}\). Then: \[\begin{align} \mat{K}^\T\mat{B}^\T\mat{S} + \mat{S}\mat{B}\mat{K} &= \mat{S}\mat{B}(\mat{R}^{-1})^\T\mat{B}^\T\mat{S} + \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S} \notag\\ &= 2\mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S} \quad \text{(using } \mat{R} = \mat{R}^\T \text{)} \\ &\quad - \mat{S}\mat{B}\mat{R}^{-1}\mat{R}(\mat{R}^{-1})^\T\mat{B}^\T\mat{S}. \label{eq:ch6:trick} \end{align}\]
A cleaner approach is to note that: \[\begin{equation} \mat{K}^\T\mat{R}\mat{K} = \mat{S}\mat{B}(\mat{R}^{-1})^\T\mat{R}\mat{R}^{-1}\mat{B}^\T\mat{S} = \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S}. \label{eq:ch6:K-R-K} \end{equation}\]
Also: \[\begin{equation} \mat{K}^\T\mat{B}^\T\mat{S} + \mat{S}\mat{B}\mat{K} = \mat{S}\mat{B}\mat{K} + \mat{K}^\T\mat{B}^\T\mat{S} = (\mat{S}\mat{B}\mat{K} + (\mat{S}\mat{B}\mat{K})^\T) = \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S} + \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S} = 2\mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S}. \label{eq:ch6:symmetric-form} \end{equation}\]
Using \(\mat{K} = \mat{R}^{-1}\mat{B}^\T\mat{S}\), we obtain: \[\begin{align} \mat{K}^\T\mat{B}^\T\mat{S} &= \mat{S}\mat{B}(\mat{R}^{-1})^\T\mat{B}^\T\mat{S} = \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S}, \\ \mat{S}\mat{B}\mat{K} &= \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S}. \end{align}\]
So indeed: \[\begin{equation} \mat{K}^\T\mat{B}^\T\mat{S} + \mat{S}\mat{B}\mat{K} = 2\mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S}. \label{eq:ch6:key-sum} \end{equation}\]
Substituting back: \[\begin{align} \mat{A}_{\text{cl}}^\T\mat{S} + \mat{S}\mat{A}_{\text{cl}} + \dot{\mat{S}} &= -\mat{Q} + \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S} - 2\mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S} \notag\\ &= -\mat{Q} - \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S} \notag\\ &= -\mat{Q} - \mat{K}^\T\mat{R}\mat{K}. \label{eq:ch6:final-identity} \end{align}\]
Therefore: \[\begin{equation} \dot{V} = \dx^\T[\mat{A}_{\text{cl}}^\T\mat{S} + \mat{S}\mat{A}_{\text{cl}} + \dot{\mat{S}}]\dx = -\dx^\T(\mat{Q} + \mat{K}^\T\mat{R}\mat{K})\dx. \label{eq:ch6:ct-value-drop} \end{equation}\]
Exponential Decay in Continuous Time
If \(\mat{Q} + \mat{K}^\T\mat{R}\mat{K} \succeq \alpha\mat{I}\) and \(\underline{s}\mat{I} \preceq \mat{S} \preceq \overline{s}\mat{I}\), then: \[\begin{equation} \dot{V} \leq -\alpha\norm{\dx}^2 \leq -\frac{\alpha}{\overline{s}}V. \label{eq:ch6:ct-decay-ineq} \end{equation}\]
Solving the differential inequality: \[\begin{equation} V(t) \leq e^{-\lambda t}V(0), \quad \text{where} \quad \lambda = \frac{\alpha}{2\overline{s}}. \label{eq:ch6:ct-exp-decay} \end{equation}\]
Converting to Euclidean norm: \[\begin{equation} \norm{\dx(t)}^2 \leq \frac{\overline{s}}{\underline{s}}e^{-2\lambda t}\norm{\dx(0)}^2, \label{eq:ch6:ct-euclidean-decay} \end{equation}\] or: \[\begin{equation} \norm{\dx(t)} \leq \sqrt{\frac{\overline{s}}{\underline{s}}} \, e^{-\lambda t}\norm{\dx(0)}. \label{eq:ch6:ct-norm-decay-final} \end{equation}\]
Comparison With Contraction LMI
Comparing~\(\eqref{eq:ch6:ct-value-drop}\) with the contraction condition from Chapter~\(\ref{ch:contraction}\): \[\begin{equation} \mat{A}_{\text{cl}}^\T\mat{M} + \mat{M}\mat{A}_{\text{cl}} + \dot{\mat{M}} \preceq -2\lambda\mat{M}, \label{eq:ch6:contraction-lmi} \end{equation}\] we see that setting \(\mat{M} = \mat{S}\) and \(\lambda = \alpha/(2\overline{s})\) gives: \[\begin{equation} \mat{A}_{\text{cl}}^\T\mat{S} + \mat{S}\mat{A}_{\text{cl}} + \dot{\mat{S}} = -(\mat{Q} + \mat{K}^\T\mat{R}\mat{K}) \preceq -\alpha\mat{I} \preceq -\frac{\alpha}{\overline{s}}\mat{S}. \label{eq:ch6:contraction-satisfied} \end{equation}\]
This is precisely the contraction condition: the Riccati solution simultaneously certifies optimality and contraction.
Robustness Margins From LQR
Gain Margin
Gain Margin of Continuous-Time LQR
The continuous-time LQR controller with gain \(\mat{K} = \mat{R}^{-1}\mat{B}^\T\mat{S}\) guarantees stability under any multiplicative perturbation \(\mat{K}_\eta = \eta\mat{K}\) for \(\eta \in (1/2, \infty)\).
Proof. Consider the perturbed closed-loop system \(\dot{\x} = (\mat{A} - \eta\mat{B}\mat{K})\x\). Using the value function \(V = \x^\T\mat{S}\x\): \[\begin{align} \dot{V} &= \x^\T[(\mat{A} - \eta\mat{B}\mat{K})^\T\mat{S} + \mat{S}(\mat{A} - \eta\mat{B}\mat{K})]\x \notag\\ &= \x^\T[\mat{A}^\T\mat{S} + \mat{S}\mat{A} - \eta(\mat{K}^\T\mat{B}^\T\mat{S} + \mat{S}\mat{B}\mat{K})]\x. \label{eq:ch6:eta-V-dot} \end{align}\]
From the ARE: \[\begin{equation} \mat{A}^\T\mat{S} + \mat{S}\mat{A} = -\mat{Q} - \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S}. \label{eq:ch6:are-form} \end{equation}\]
And we showed: \[\begin{equation} \mat{K}^\T\mat{B}^\T\mat{S} + \mat{S}\mat{B}\mat{K} = 2\mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S} = 2\mat{K}^\T\mat{R}\mat{K}. \label{eq:ch6:control-term} \end{equation}\]
Thus: \[\begin{align} \dot{V} &= \x^\T[-\mat{Q} - \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S} - 2\eta\mat{K}^\T\mat{R}\mat{K}]\x \notag\\ &= -\x^\T[\mat{Q} + (2\eta)\mat{K}^\T\mat{R}\mat{K} + \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S}]\x. \label{eq:ch6:eta-perturbed-V} \end{align}\]
The condition \(\dot{V} < 0\) requires the bracketed term to be positive definite. Assuming \(\eta > 1/2\), we have \(2\eta > 1\), which means the positive definite matrix \(2\eta\mat{K}^\T\mat{R}\mat{K}\) dominates, and \(\dot{V} < 0\) is maintained.
Phase Margin
Phase Margin of Continuous-Time LQR
The continuous-time LQR controller guarantees a phase margin of at least \(60^\circ\) on both sides of the nominal open-loop transfer function.
The complete derivation connecting the return-difference inequality to the specific \(60^\circ\) bound can be found in Anderson & Moore, Optimal Control: Linear Quadratic Methods (1990), Chapter 5, or Safonov & Athans (1977). The key insight: the Riccati solution makes the closed-loop system “strictly positive real,” and by the Passivity Theorem, the allowable phase perturbation is bounded to \(60^\circ\). This is a well-known classical result; we state it here as a reference, citing the established sources rather than reproducing the full derivation.
Disk Margin
The disk margin, or \(\mu\)-margin, quantifies robustness in the complex plane. For LQR, the loop transfer matrix satisfies a disk margin condition: the poles of the closed-loop system remain in a half-plane with guaranteed clearance from the origin. This is a consequence of the Riccati solution being positive definite and the ARE’s structure.
Implication for Design
These robustness margins explain why LQR is ubiquitously used in practice: not only does it minimize a well-defined cost, but it automatically provides guardrails against model uncertainty, sensor noise, and actuator saturation. The engineer need not separately verify robustness—the LQR computation provides it.
Condition Number and Its Effect on Contraction Rate
Definition and Interpretation
Define the condition number of the Riccati solution: \[\begin{equation} \kappa(\mat{S}) := \frac{\lambda_{\max}(\mat{S})}{\lambda_{\min}(\mat{S})} = \frac{\overline{s}}{\underline{s}}. \label{eq:ch6:kappa-S} \end{equation}\]
A well-conditioned matrix has \(\kappa(\mat{S}) \approx 1\); a poorly conditioned matrix has \(\kappa(\mat{S}) \gg 1\).
Ill-Conditioning Signals Loss of Controllability
When \(\kappa(\mat{S})\) is large, it often indicates that: 1. The system is nearly uncontrollable in some direction (e.g., the angular velocity of the pendulum is hard to control). 2. The Riccati matrix has very different scales along its principal directions. 3. Perturbations in one direction (corresponding to \(\lambda_{\max}(\mat{S})\)) are much more costly than in another (corresponding to \(\lambda_{\min}(\mat{S})\)).
In the limit of complete uncontrollability, \(\lambda_{\min}(\mat{S}) \to 0\), and \(\kappa(\mat{S}) \to \infty\).
Condition Number and Contraction Rate
From the discrete-time analysis: \[\begin{equation} \rho = 1 - \frac{\alpha}{\overline{s}} = 1 - \frac{\alpha}{\kappa(\mat{S})\underline{s}}. \label{eq:ch6:rho-vs-kappa} \end{equation}\]
As \(\kappa(\mat{S})\) increases, \(\overline{s}\) increases (for fixed \(\underline{s}\) and \(\alpha\)), and the contraction rate \(\rho\) increases (approaches 1), meaning convergence slows. More precisely: \[\begin{equation} 1 - \rho = \frac{\alpha}{\overline{s}} \approx \frac{\alpha}{\kappa(\mat{S})\underline{s}}. \label{eq:ch6:contraction-loss} \end{equation}\]
The number of time steps to reduce the value by a factor of \(e\) is: \[\begin{equation} \tau = \frac{1}{\ln(1/\rho)} \approx \frac{1}{\alpha/\overline{s}} = \frac{\overline{s}}{\alpha} = \kappa(\mat{S})\frac{\underline{s}}{\alpha}. \label{eq:ch6:time-constant} \end{equation}\]
Conclusion: A large condition number \(\kappa(\mat{S})\) is a warning sign. It means the system has directions that are hard to control, and convergence is slower. The designer should either: 1. Increase the state penalty \(\mat{Q}\) in the weakly controllable direction to reduce \(\overline{s}/\underline{s}\). 2. Accept slower convergence but maintain open-loop robustness. 3. Redesign the system (e.g., add actuators) to improve controllability.
Connections to H-Infinity and Robust Control
The H-Infinity Problem
The \(H^\infty\) optimal control problem is: given a system with disturbances \(\mat{w}\), \[\begin{equation} \dot{\x} = \mat{A}\x + \mat{B}_u\uvec + \mat{B}_w\mat{w}, \label{eq:ch6:hinf-system} \end{equation}\] with controlled output: \[\begin{equation} \mat{z} = \begin{bmatrix} \mat{C}_1 \\ \mat{D}_{uu} \end{bmatrix}\x + \begin{bmatrix} 0 \\ \mat{D}_{uw} \end{bmatrix}\mat{w}, \label{eq:ch6:hinf-output} \end{equation}\] find a feedback controller \(\mat{K}\) that minimizes the \(L_2\)-induced norm from disturbances to outputs: \[\begin{equation} \gamma^\star = \min_{\mat{K}} \max_{\mat{w} \neq 0} \frac{\|\mat{z}\|_{L_2}}{\|\mat{w}\|_{L_2}}. \label{eq:ch6:hinf-objective} \end{equation}\]
The H-Infinity Riccati Equation
The optimal \(H^\infty\) controller satisfies a Riccati equation similar in structure to the LQR ARE: \[\begin{equation} 0 = \mat{A}^\T\mat{S}_\infty + \mat{S}_\infty\mat{A} + \mat{C}_1^\T\mat{C}_1 - \mat{S}_\infty(\mat{B}_u\mat{D}_{uu}^{-1}\mat{B}_u^\T - \gamma^{-2}\mat{B}_w\mat{B}_w^\T)\mat{S}_\infty, \label{eq:ch6:hinf-are} \end{equation}\] where \(\gamma\) is the worst-case \(L_2\) gain.
Duality in the Robust Setting
The duality theorem extends to \(H^\infty\) control: the Riccati solution \(\mat{S}_\infty\) serves as both: 1. Optimal performance certificate: the value function \(V = \x^\T\mat{S}_\infty\x\) guarantees that the \(L_2\) gain from disturbances to outputs is at most \(\gamma\). 2. Contraction metric under disturbances: the time-derivative of \(V\) is negative even under the worst-case disturbance, ensuring exponential stability.
The key insight is the disturbance direction in the Riccati equation: \(\gamma^{-2}\mat{B}_w\mat{B}_w^\T\) enters with a negative sign (compared to the control input \(\mat{B}_u\) which is positive). This asymmetry means the controller is simultaneously designed to: - Minimize control energy (via \(\mat{B}_u\)). - Maximize disturbance rejection (via \(\mat{B}_w\) with the worst-case scaling).
Worked Example: Pendulum With Wind Disturbance
Consider the pendulum with a wind force disturbance \(w\) entering as: \[\begin{equation} \ddot\theta = \theta - u + w. \label{eq:ch6:pendulum-wind} \end{equation}\]
The \(H^\infty\) design seeks a controller that minimizes the gain from wind \(w\) to angle tracking error \(\theta\). The Riccati solution provides a controller that is simultaneously as energy-efficient as possible while rejecting the worst-case wind disturbance. The duality theorem guarantees that this single controller solves both the performance and robustness objectives.
Time-Varying vs. Steady-State Duality
The Discrete Riccati Equation (DRE)
For finite-horizon problems (e.g., trajectory optimization over a \(T\)-step horizon), the cost-to-go matrix evolves backward according to the discrete Riccati recursion: \[\begin{equation} \mat{S}_k = \mat{Q}_k + \mat{A}_k^\T\mat{S}_{k+1}\mat{A}_k - \mat{A}_k^\T\mat{S}_{k+1}\mat{B}_k(\mat{R}_k + \mat{B}_k^\T\mat{S}_{k+1}\mat{B}_k)^{-1}\mat{B}_k^\T\mat{S}_{k+1}\mat{A}_k, \label{eq:ch6:dre-discrete} \end{equation}\] with \(\mat{S}_T = \mat{S}_f\) (a terminal cost, often zero).
At each time \(k\), the contraction property holds in the tangent space: \[\begin{equation} \norm{\dx_{k+1}}_{\mat{S}_{k+1}} \leq \rho_k \norm{\dx_k}_{\mat{S}_k}, \label{eq:ch6:time-varying-contraction} \end{equation}\] where \(\rho_k = 1 - \alpha_k/\overline{s}_k\) is time-dependent.
The Algebraic Riccati Equation (ARE)
For infinite-horizon, time-invariant LQR, the Riccati recursion converges to a fixed point, the ARE: \[\begin{equation} 0 = \mat{Q} + \mat{A}^\T\mat{S} + \mat{S}\mat{A} - \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^\T\mat{S}. \label{eq:ch6:are-again} \end{equation}\]
The solution \(\mat{S}\) is constant in time, and the contraction rate \(\rho\) is also constant: \[\begin{equation} \norm{\dx_k}_\mat{S} \leq \rho^k \norm{\dx_0}_\mat{S}, \label{eq:ch6:steady-state-contraction} \end{equation}\] with \(\rho = 1 - \alpha/\overline{s}\) independent of \(k\).
Geometric Interpretation
- Time-varying case (DRE): The Riccati matrix \(\mat{S}_k\) changes at each step, reflecting the changing “cost-to-go” as we move backward in time from the final state. The tangent space metric is adapted to the current time: at times near the terminal condition, the metric emphasizes directions that are expensive to steer; far from the terminal, the metric is influenced by asymptotic stability.
- Steady-state case (ARE): The Riccati matrix reaches a time-invariant steady state. The metric no longer changes, and the contraction property is uniform in time. This is appropriate for infinite-horizon problems or problems where the horizon is much longer than the natural time scale of the system.
In both cases, the duality theorem holds: the Riccati solution (time-varying or steady-state) is simultaneously the optimal cost-to-go and the contraction metric.
Tangent Bundle Geometry and Geodesics
The Tangent Bundle
A trajectory on a state-space manifold is a curve \(\gamma(t) = (x_1(t), \ldots, x_n(t))\). At each point \(x\) on the trajectory, the tangent space \(T_x M\) is the space of infinitesimal perturbation vectors \(\dx\). The tangent bundle \(TM\) is the union of all tangent spaces: \[\begin{equation} TM = \bigcup_{x \in M} T_x M = \{(x, \dx) : x \in M, \dx \in \mathbb{R}^n\}. \label{eq:ch6:tangent-bundle} \end{equation}\]
Riemannian Metric on the Tangent Bundle
A Riemannian metric is a positive-definite inner product on the tangent space at each point. In our context, the Riccati solution \(\mat{S}(x)\) (which may depend on the state \(x\) for nonlinear systems, but is constant for linear systems) defines a Riemannian metric: \[\begin{equation} \langle \dx_1, \dx_2 \rangle_\mat{S} = \dx_1^\T\mat{S}\,\dx_2. \label{eq:ch6:riemannian-metric} \end{equation}\]
The length of a tangent vector is: \[\begin{equation} \|\dx\|_\mat{S} = \sqrt{\dx^\T\mat{S}\,\dx}. \label{eq:ch6:metric-norm} \end{equation}\]
Geodesics and Optimal Control
In Riemannian geometry, a geodesic is a curve of shortest length (in the metric sense). For our system, the optimal trajectories (those minimizing the LQR cost) are precisely the geodesics with respect to the Riccati metric.
Here’s why: the LQR cost functional \[\begin{equation} J = \int_0^T \left(\x^\T\mat{Q}\x + \uvec^\T\mat{R}\uvec\right) dt \label{eq:ch6:lqr-cost-again} \end{equation}\] can be rewritten, using the Riccati solution, as a metric-weighted length: \[\begin{equation} J \approx \int_0^T \norm{\dot{\gamma}}_\mat{S} \, dt, \label{eq:ch6:geodesic-length} \end{equation}\] where the precise relationship depends on the details of the derivation. The Euler-Lagrange equations for minimizing this metric-weighted path length are equivalent to the optimal control equations.
Convergence of Neighboring Trajectories
A key property of geodesics is that they are stable (in a local sense): neighboring geodesics maintain bounded separation. In our setting, the contraction property guarantees stronger convergence: \[\begin{equation} \|\Delta \gamma(t)\|_\mat{S} \leq \rho^{t/\tau} \|\Delta \gamma(0)\|_\mat{S}, \label{eq:ch6:geodesic-convergence} \end{equation}\] where \(\Delta\gamma\) is the difference between two nearby optimal trajectories, and \(\tau\) is the natural time scale.
This is a remarkable property: not only are the optimal trajectories geodesics (shortest paths), but they also attract nearby geodesics exponentially fast. This is what makes the LQR solution so robust and widely applicable.
Visualization in 2D
For a 2D system (e.g., the pendulum with angle and angular velocity), the Riccati metric defines an elliptical “distance metric.” The level sets of the Lyapunov function \(V = \dx^\T\mat{S}\,\dx = 1\) form ellipses. The principal axes of the ellipse are the eigenvectors of \(\mat{S}\), with radii proportional to \(1/\sqrt{\lambda_i(\mat{S})}\).
A smaller (more negative) eigenvalue means a longer radius in that direction: perturbations along that eigenvector are “cheaper” in the metric sense. Conversely, a larger eigenvalue means a shorter radius: perturbations are expensive.
Under the optimal feedback, all perturbation ellipses shrink at the rate \(\rho\) per time step, and the ellipses align with the eigenvectors of \(\mat{A}_{\text{cl}}\).
Practical Design With Duality
Integrated Design Workflow
The duality theorem suggests a unified design methodology:
Integrated Stability-Optimality Design: 1. Model the system: discretize or linearize to get \((\mat{A}, \mat{B})\). 2. Choose weights: select \(\mat{Q}\) and \(\mat{R}\) based on performance objectives (state tracking, control effort). 3. Solve the Riccati equation: compute \(\mat{S}\) (via the ARE or DRE). 4. Extract the gain: \(\mat{K} = (\mat{R} + \mat{B}^\T\mat{S}\mat{B})^{-1}\mat{B}^\T\mat{S}\mat{A}\). 5. Verify contraction margins: - Compute \(\alpha = \lambda_{\min}(\mat{Q} + \mat{K}^\T\mat{R}\mat{K})\). - Compute \(\overline{s} = \lambda_{\max}(\mat{S})\) and \(\underline{s} = \lambda_{\min}(\mat{S})\). - Contraction rate: \(\rho = 1 - \alpha/\overline{s}\). - Condition number: \(\kappa(\mat{S}) = \overline{s}/\underline{s}\). 6. Decision point: - If \(\rho\) is close to 1 (say, \(\rho > 0.9\)) or \(\kappa(\mat{S})\) is large (say, \(> 100\)), the system has weak controllability. Increase \(\mat{Q}\) in the hard-to-control direction and re-solve. - If \(\rho\) is small (say, \(< 0.3\)) and \(\kappa(\mat{S})\) is reasonable, the design is aggressive and robust. Proceed to implementation. - Otherwise, the design is balanced. Implement and tune on hardware if needed. 7. Robustness certification: the gain margins [1/2, \(\infty\)) and phase margin \(\geq\) \(60^\circ\) apply automatically.
When Duality Breaks Down
The duality requires several conditions:
Limitations of the Duality: 1. Bounded Riccati: If the system loses controllability near a trajectory, the Riccati matrix can become unbounded or ill-conditioned. This is detected by monitoring \(\kappa(\mat{S})\). If \(\kappa(\mat{S}) > 10^6\) or \(\lambda_{\min}(\mat{S}) < 10^{-6}\), the system is nearly uncontrollable. 2. Meaningful weights: The weights \(\mat{Q}\) and \(\mat{R}\) must reflect true performance objectives. If weights are chosen arbitrarily (e.g., to make the Riccati well-conditioned), the resulting controller may be suboptimal in practice. 3. Linearity and perturbations: The duality holds exactly in the tangent space (linearized dynamics). For large perturbations or highly nonlinear systems, the linear theory breaks down. Nonlinear generalizations (e.g., differential games with Bellman equations) exist but are more complex. 4. Measurement and state estimation: The duality assumes perfect state knowledge (\(\uvec = -\mat{K}\x\)). With measurement noise or partial observability, a separate state estimator (Kalman filter) is needed, and the duality extends in a more subtle way (separation principle). 5. Sampling and discretization: When continuous-time controllers are implemented on digital hardware, discretization errors can break the duality. Careful Tustin or zero-order-hold discretization is required.
The Estimation–Control Duality
The duality developed so far connects stability to optimality on the control side. But the Riccati equation has a mirror image: the estimation Riccati, which determines how uncertainty propagates through the system and how observations reduce it. This section develops the estimation side of the duality.
The Dual Riccati Equation
Consider the system with process noise and output measurement: \[\begin{align} \dot{\x} &= \mat{A}(t)\x + \mat{B}(t)\uvec + \mat{G}_w w, \label{eq:ch6:process-noise} \\ \vec{y} &= \mat{C}(t)\x + v, \label{eq:ch6:measurement} \end{align}\] where \(w \sim \mathcal{N}(0, \mat{W})\) and \(v \sim \mathcal{N}(0, \mat{V})\) are process and measurement noise. The estimation Riccati equation governs the covariance \(\mat{P}(t)\) of the optimal (Kalman) state estimate: \[\begin{equation} \dot{\mat{P}} = \mat{A}\mat{P} + \mat{P}\mat{A}^T + \mat{G}_w \mat{W} \mat{G}_w^T - \mat{P}\mat{C}^T \mat{V}^{-1} \mat{C}\mat{P}. \label{eq:ch6:estimation-riccati} \end{equation}\]
Compare with the control Riccati: \[\begin{equation} -\dot{\mat{S}} = \mat{A}^T\mat{S} + \mat{S}\mat{A} + \mat{Q} - \mat{S}\mat{B}\mat{R}^{-1}\mat{B}^T\mat{S}. \label{eq:ch6:control-riccati-recall} \end{equation}\]
The duality is exact: the estimation Riccati for \((\mat{A}, \mat{C}, \mat{W}, \mat{V})\) is the control Riccati for \((\mat{A}^T, \mat{C}^T, \mat{G}_w \mat{W} \mat{G}_w^T, \mat{V})\). The estimation gain \(\mat{L} = \mat{P}\mat{C}^T\mat{V}^{-1}\) is the dual of the control gain \(\mat{K} = \mat{R}^{-1}\mat{B}^T\mat{S}\).
The Extended Kalman Filter as Tangent-Space Estimation
For nonlinear systems \(\dot{\x} = f(\x) + G(\x)\uvec + w\), the Extended Kalman Filter (EKF) propagates the state estimate and covariance using the tangent-space (linearized) dynamics: \[\begin{align} \dot{\hat{\x}} &= f(\hat{\x}) + G(\hat{\x})\uvec + \mat{L}(t)\bigl(\vec{y} - h(\hat{\x})\bigr), \label{eq:ch6:ekf-state} \\ \dot{\mat{P}} &= \mat{A}(t)\mat{P} + \mat{P}\mat{A}(t)^T + \mat{G}_w \mat{W} \mat{G}_w^T - \mat{P}\mat{C}(t)^T \mat{V}^{-1} \mat{C}(t)\mat{P}, \label{eq:ch6:ekf-cov} \end{align}\] where \(\mat{A}(t) = \frac{\partial f}{\partial \x}\big|_{\hat{\x}(t)}\) and \(\mat{C}(t) = \frac{\partial h}{\partial \x}\big|_{\hat{\x}(t)}\).
The EKF is the estimation analogue of DDP: both linearize the system along a trajectory (the estimated trajectory for EKF, the nominal trajectory for DDP) and solve a Riccati equation in the tangent space. The key difference is the direction of information flow: DDP propagates cost backward (from the terminal condition); the EKF propagates uncertainty forward (from the initial condition).
Contraction-Based Observer Design
The EKF has no convergence guarantees for nonlinear systems—the estimate may diverge if the linearization is poor. Contraction theory provides a rigorous alternative.
An observer \(\dot{\hat{\x}} = f(\hat{\x}) + G(\hat{\x})\uvec + L(\hat{\x})(\vec{y} - h(\hat{\x}))\) is a contracting observer if the observer error dynamics \(\dot{e} = \dot{\x} - \dot{\hat{\x}}\) are contracting in some metric \(\mat{M}(\x)\): \[\begin{equation} \sym\left(\mat{M} \left[\frac{\partial f}{\partial \x} - L(\x)\frac{\partial h}{\partial \x}\right]\right) + \dot{\mat{M}} \preceq -2\lambda \mat{M}. \label{eq:ch6:contracting-observer} \end{equation}\]
This is precisely the contraction condition of Chapter~\(\ref{ch:contraction}\), applied to the observer error dynamics. The observer gain \(L(\x)\) plays the same role as the feedback gain \(\mat{K}\) in control contraction metrics: it shapes the closed-loop dynamics to achieve contraction.
Observer-Controller Duality
The observer design problem — find \(L(\x)\) such that~\(\eqref{eq:ch6:contracting-observer}\) holds — is the dual of the controller design problem via the substitutions \(\mat{A} \leftrightarrow \mat{A}^T\), \(\mat{B} \leftrightarrow \mat{C}^T\), \(\mat{K} \leftrightarrow L^T\). A contraction metric for the observer is a contraction metric for the dual controller, and vice versa.
The Observability Gramian
The observability Gramian \(\mathcal{W}_o(t_0, t_1)\) quantifies how much information the output \(\vec{y}\) reveals about the initial state \(\x(t_0)\): \[\begin{equation} \mathcal{W}_o(t_0, t_1) = \int_{t_0}^{t_1} \Phi(s, t_0)^T \mat{C}(s)^T \mat{C}(s) \Phi(s, t_0) \, \dd s, \label{eq:ch6:obs-gramian} \end{equation}\] where \(\Phi\) is the state transition matrix from Chapter~\(\ref{ch:variational}\).
The system is observable on \([t_0, t_1]\) if \(\mathcal{W}_o \succ 0\). The eigenvectors of \(\mathcal{W}_o\) reveal the most and least observable state directions; the condition number \(\kappa(\mathcal{W}_o)\) measures the anisotropy of observability.
For the golf swing application (Chapter~\(\ref{ch:applications}\)), the observability Gramian reveals which state variables can be reliably estimated from sensor measurements. Joint angles near the shoulder are typically well-observed (multiple muscle attachments create redundant signals); wrist angular velocity and shaft bending are poorly observed (few sensors, high velocity). The condition number of \(\mathcal{W}_o\) quantifies this anisotropy.
The Separation Principle
When both observer and controller are available, the celebrated separation principle allows independent design:
Separation Principle
For linear time-invariant systems, the optimal controller using estimated states \(\uvec = -\mat{K}\hat{\x}\) achieves the same closed-loop eigenvalues as the full-state feedback \(\uvec = -\mat{K}\x\), with additional eigenvalues from the observer error dynamics. The controller gain \(\mat{K}\) and observer gain \(\mat{L}\) can be designed independently.
For nonlinear systems, the separation principle does not hold in general. However, when both the controller and observer are designed via contraction metrics, the overall system (controller + observer) contracts at a rate determined by the minimum of the controller and observer contraction rates: \(\lambda_{\text{overall}} \geq \min(\lambda_{\text{control}}, \lambda_{\text{observer}})\). This is a practical nonlinear extension of the separation principle.
Chapter Summary
- Historical foundation: Kalman (1960) discovered that the LQR cost-to-go matrix is a Lyapunov function for the closed-loop system. Willems (1972) formalized this via dissipation inequalities. Anderson and Moore highlighted the robustness margins of LQR.
- Discrete duality: The value decrease identity \(\Delta V = -\dx^\T(\mat{Q} + \mat{K}^\T\mat{R}\mat{K})\dx\) is both a proof of optimality and a contraction certificate. The decay rate \(\rho = 1 - \alpha/\overline{s}\) characterizes convergence speed.
- Continuous duality: The time-derivative \(\dot{V} = -\dx^\T(\mat{Q} + \mat{K}^\T\mat{R}\mat{K})\dx\) follows directly from the DRE and the closed-loop dynamics. Exponential decay follows with rate \(\lambda = \alpha/(2\overline{s})\).
- Robustness margins: LQR guarantees gain margin [1/2, \(\infty\)), phase margin \(\geq\) \(60^\circ\), and disk margin—all consequences of the Riccati structure and the positive-definiteness of \(\mat{S}\).
- Condition number: The condition number \(\kappa(\mat{S}) = \overline{s}/\underline{s}\) quantifies the difficulty of controlling the system. Large \(\kappa\) indicates poor controllability and slow convergence.
- H\(\infty\) extension: The duality extends to robust control: the H\(\infty\) Riccati solution is simultaneously an optimal performance certificate and a contraction metric under worst-case disturbances.
- Time-varying vs. steady-state: The DRE provides time-adapted metrics for finite-horizon problems; the ARE provides a constant metric for infinite-horizon, time-invariant systems.
- Geometric interpretation: The Riccati matrix defines a Riemannian metric on the tangent bundle. Optimal trajectories are geodesics, and the contraction property ensures neighboring geodesics converge exponentially.
- Practical workflow: Design with duality integrated: solve LQR once, verify contraction margins and condition number, and adjust weights if needed. The single computation certifies both optimality and robustness.
- Limitations: Duality requires bounded, well-conditioned Riccati matrices; meaningful weights; and linearity of perturbations. It breaks down near singular arcs, under large perturbations, or with hidden states.
This chapter has shown that optimality and stability are not independent concerns but two facets of a single geometric object. The Riccati matrix, solved once, provides both the controller and the robustness certificate. This is the power of duality: one computation, two major objectives accomplished.
Exercises
- Value decrease identity. For the discrete pendulum in the worked example, verify the value decrease identity numerically: compute \(V_{k+1} - V_k\) and \(-\dx_k^T(\mat{Q} + \mat{K}^T\mat{R}\mat{K})\dx_k\) and confirm equality.
- Robustness margins. For the continuous-time system \(\dot{x} = ax + bu\) with LQR weights \(Q = 1\), \(R = 1\), compute the gain margin and phase margin as functions of \(a\) and \(b\). Verify the theoretical bounds \([1/2, \infty)\) and \(\geq 60°\).
- Dual Riccati. For the system \(\mat{A} = \begin{bmatrix} 0 & 1 \\ -1 & -0.1 \end{bmatrix}\), \(\mat{C} = [1, 0]\), \(W = 0.1\), \(V = 1\), solve the estimation ARE for \(\mat{P}\). Compare with the control ARE for \((\mat{A}^T, \mat{C}^T, W, V)\).
- Condition number and convergence. For the 2D system in the worked example, vary \(\mat{Q}\) from \(0.1\mat{I}\) to \(100\mat{I}\) while keeping \(\mat{R} = \mat{I}\). Plot \(\kappa(\mat{S})\) and the contraction rate \(\lambda\) as functions of \(\|\mat{Q}\|\). Explain the observed relationship.
- EKF implementation. Implement the EKF for the nonlinear pendulum \(\ddot\theta + \sin\theta = u + w\) with measurement \(y = \theta + v\). Compare the estimated trajectory with the true trajectory for different noise levels.
- Observability Gramian. For the linearized pendulum at the upright equilibrium, compute the observability Gramian \(\mathcal{W}_o(0, T)\) for \(T = 1, 5, 10\). How does its condition number change? What does this imply about the observability of the system?
- H\(\infty\) vs LQR. Solve both the LQR and H\(\infty\) Riccati equations for the double integrator with disturbance input. Compare the resulting contraction rates and gain magnitudes.
- Separation principle verification. For the discrete pendulum, implement the separated controller (LQR gain) and observer (Kalman gain). Verify that the closed-loop eigenvalues are the union of the controller and observer eigenvalues.