The Linear-Quadratic Regulator: Differential and Algebraic Riccati Equations

Intuition Beginner

Imagine balancing a tall pole on your hand, or steering a drone back to a hover after a gust. At every instant you are off from where you want to be, and you can push to correct it. Two things annoy you. First, being off-target is bad: a leaning pole, a drifting drone. Second, pushing hard is also bad: it burns energy, strains motors, and feels jerky. The linear-quadratic regulator is the recipe for the best trade-off between these two annoyances, when the system responds in a straight-line way to your pushes and the penalties grow with the square of how far off you are and how hard you push.

Why square the penalties? Squaring makes small errors cheap and large errors very expensive, so the recipe works hard to kill big deviations and tolerates tiny ones. Squaring also has a clean payoff: the best control turns out to be a simple rule. You measure how far the system is from its target, multiply by a fixed table of numbers, and push by that amount. That rule is called feedback, and the table is the gain.

The deep surprise is where the gain comes from. To know the best push right now you would think you need to plan the entire future. It collapses into one object: a matrix that records how costly it is to be in each state, given that you will steer optimally from here on. Computing that matrix is the whole game, and it satisfies a special equation named after Riccati.

The Riccati equation is the bridge from "what do all my future choices cost" to "what single push is best now." Solve it once, read off the gain table, and you have an autopilot.

Visual Beginner

The picture shows the loop the regulator runs and where the gain table sits inside it.

Read the table left to right as one cycle, repeated forever. The system has a current state (how far off you are). The gain table turns that state into a push. The push, plus the system's own drift, produces the next state, and the cycle repeats. The gain table is fixed once you have solved the Riccati equation, so running the autopilot is cheap.

stage	what happens	object
measure	read how far off you are	state $x$
decide	multiply state by the gain table	push $u=-Kx$
act	push moves the system, which also drifts on its own	$x$ changes
repeat	the new state feeds back to measure	loop closes

The takeaway: the optimal controller is just "push in proportion to the error," and the only hard work, finding the proportion, is done in advance by solving the Riccati equation. Once you have the gain, controlling the system costs one multiplication per step.

Worked example Beginner

Take the simplest possible system: a single number $x$ that drifts according to $\dot{x}=u$, where $u$ is your push. There is no built-in motion, so $A=0$ and $B=1$. You pay for being off-target at rate $x$ squared and for pushing at rate $u$ squared, so the state penalty is $Q=1$ and the control penalty is $R=1$. We want the steady, long-run best rule.

Step 1. Guess the cost-to-go has the form $V(x)=p{x}^2$ for some single positive number $p$, recording how expensive it is to sit at state $x$ and steer optimally forever after.

Step 2. The steady Riccati equation in this scalar case reads $2Ap-B^2{p}^2/R+Q=0$. Plug in $A=0$, $B=1$, $R=1$, $Q=1$:

$$0-p^2+1=0,\text{so}{p}^2=1.$$

Step 3. Pick the positive root $p=1$ (the cost-to-go must be a cost, so it cannot be negative).

Step 4. Read off the gain. The rule is $u=-Kx$ with $K=Bp/R=1\cdot 1/1=1$. So the best push is $u=-x$: push back by exactly your current error.

Step 5. Check it stabilises. With $u=-x$ the system becomes $\dot{x}=-x$, which decays smoothly to zero. The error shrinks; the autopilot works.

What this tells us: even with the future stretching to infinity, the answer is one number $p=1$ and one gain $K=1$. The Riccati equation turned an infinite planning problem into solving $p^2=1$, and the resulting controller $u=-x$ drives the state to rest.

Check your understanding Beginner

Formal definition Intermediate+

Let the plant be the linear time-invariant system $\dot{x}=Ax+Bu$ of 02.06.03, with state $x(t)\in {\mathbb{R}}^n$, control $u(t)\in {\mathbb{R}}^m$, $A\in {\mathbb{R}}^{n\times n}$, $B\in {\mathbb{R}}^{n\times m}$. Fix symmetric weight matrices $Q=Q^{\mathsf{T}}⪰0$ on the state and $R=R^{\mathsf{T}}\succ 0$ on the control (the positive-semidefinite and positive-definite cones of 01.01.13), and a terminal weight $S=S^{\mathsf{T}}⪰0$.

Definition (finite-horizon LQR). Over $[t_0,t_f]$ the cost of a state-control trajectory is

$$J(u)=x(t_f{)}^{\mathsf{T}}Sx(t_f)+{\int }_{t_0}^{t_f}(x(t{)}^{\mathsf{T}}Qx(t)+u(t{)}^{\mathsf{T}}Ru(t))dt.$$

The linear-quadratic regulator is the control $u$ minimising $J$ subject to the dynamics and the initial condition $x(t_0)=x_0$.

Definition (value function / cost-to-go). The value function is

$$V(t,x)=\underset{u[t,t_f]}{\min }[x(t_f{)}^{\mathsf{T}}Sx(t_f)+{\int }_t^{t_f}(x^{\mathsf{T}}Qx+u^{\mathsf{T}}Ru)ds],x(t)=x,$$

the least cost achievable from state $x$ at time $t$ steering optimally to $t_f$.

Definition (matrix Riccati differential equation, RDE). The Riccati differential equation for the symmetric matrix $P(t)\in {\mathbb{R}}^{n\times n}$ is the terminal-value problem

$$-\dot{P}=A^{\mathsf{T}}P+PA-PB{R}^{-1}{B}^{\mathsf{T}}P+Q,P(t_f)=S.$$

Its solution $P(t)$, run backward from $t_f$, is the Hessian of the value function: $V(t,x)=x^{\mathsf{T}}P(t)x$. The associated time-varying gain is $K(t)=R^{-1}{B}^{\mathsf{T}}P(t)$, and the optimal control is the state feedback $u(t)=-K(t)x(t)$.

Definition (continuous-time algebraic Riccati equation, CARE). Letting $t_f\to \infty$ with $S=0$ and $Q$, $R$ constant gives the infinite-horizon LQR, minimising ${\int }_0^{\infty }(x^{\mathsf{T}}Qx+u^{\mathsf{T}}Ru)dt$. Its value function is $V(x)=x^{\mathsf{T}}Px$ with $P$ the steady solution $\dot{P}=0$, i.e. the symmetric solution of the continuous-time algebraic Riccati equation

$$A^{\mathsf{T}}P+PA-PB{R}^{-1}{B}^{\mathsf{T}}P+Q=0.$$

The constant optimal gain is $K=R^{-1}{B}^{\mathsf{T}}P$ and the optimal control is the static feedback $u=-Kx$, producing the closed-loop matrix $A_{\mathrm{cl}}=A-BK=A-B{R}^{-1}{B}^{\mathsf{T}}P$.

Definition (stabilisability, detectability). The pair $(A,B)$ is stabilisable when there exists $F$ with $A-BF$ Hurwitz (every uncontrollable mode is already stable); equivalently $\mathrm{rank}[A-\lambda IB]=n$ for every $\lambda$ with $\mathrm{Re}\lambda \ge 0$. The pair $(A,C)$, with $C=Q^{1/2}$, is detectable when $(A^{\mathsf{T}},C^{\mathsf{T}})$ is stabilisable (every unobservable mode is already stable). These are the hypotheses under which the CARE has a unique stabilising solution.

Definition (discrete-time LQR and DARE). For $x_{k+1}=A{x}_k+B{u}_k$ with cost ${\sum }_{k=0}^{\infty }(x_k^{\mathsf{T}}Q{x}_k+u_k^{\mathsf{T}}R{u}_k)$, the steady value Hessian solves the discrete-time algebraic Riccati equation

$$P=A^{\mathsf{T}}PA-A^{\mathsf{T}}PB(R+B^{\mathsf{T}}PB{)}^{-1}{B}^{\mathsf{T}}PA+Q,$$

with optimal gain $K=(R+B^{\mathsf{T}}PB{)}^{-1}{B}^{\mathsf{T}}PA$ and closed loop $A-BK$.

Notation: $A,B$ are the system matrices; $Q⪰0$, $R\succ 0$ the state and control weights; $P(t)$ the Riccati matrix; $K$ the feedback gain; $V$ the value function; $A_{\mathrm{cl}}$ the closed-loop matrix; $\mathcal{H}$ (introduced below) the Hamiltonian matrix; $⪰,\succ$ the Loewner order on symmetric matrices.

Counterexamples to common slips

• Positive-definiteness of $R$ is not optional. If $R$ is only positive-semidefinite, some control directions are free, $R^{-1}$ fails to exist, and the gain $K=R^{-1}{B}^{\mathsf{T}}P$ is undefined; the minimisation can have no finite optimum. The state weight $Q$ may be merely semidefinite, but $R\succ 0$ is structural.

• The CARE has many symmetric solutions; the regulator wants exactly one. The equation $A^{\mathsf{T}}P+PA-PB{R}^{-1}{B}^{\mathsf{T}}P+Q=0$ is quadratic in $P$, so it generally has a whole lattice of symmetric solutions (one per Lagrangian invariant subspace of the Hamiltonian). Only the stabilising solution — the one making $A-B{R}^{-1}{B}^{\mathsf{T}}P$ Hurwitz — is the optimal cost-to-go; picking the wrong root gives an unstable closed loop.

• Stabilisability is weaker than controllability and is the right hypothesis. A system with a stable but uncontrollable mode is not controllable, yet the infinite-horizon regulator is still well posed because the offending mode decays on its own. Demanding full controllability would needlessly exclude such systems.

Key theorem with proof Intermediate+

Theorem (LQR optimality via dynamic programming and the Riccati equation; Bertsekas §4.1 ^{[source pending]}; Anderson-Moore Ch. 2-3 ^{[source pending]}). For the finite-horizon problem with $R\succ 0$, $Q⪰0$, $S⪰0$, let $P(t)$ be the symmetric solution of the Riccati differential equation $-\dot{P}=A^{\mathsf{T}}P+PA-PB{R}^{-1}{B}^{\mathsf{T}}P+Q$ with $P(t_f)=S$. Then $P(t)$ exists, is symmetric and positive-semidefinite on all of $[t_0,t_f]$, the value function is the quadratic form $V(t,x)=x^{\mathsf{T}}P(t)x$, and the unique optimal control is the linear feedback $u^{\star }(t)=-R^{-1}{B}^{\mathsf{T}}P(t)x(t)$.

Proof. The value function obeys the Hamilton-Jacobi-Bellman (HJB) equation: writing the running cost $\ell (x,u)=x^{\mathsf{T}}Qx+u^{\mathsf{T}}Ru$, the principle of optimality gives, for the value $V(t,x)$,

$$-\frac{\partial V}{\partial t}=\underset{u\in {\mathbb{R}}^m}{\min }[\ell (x,u)+{\nabla }_{x}V\cdot (Ax+Bu)],V(t_f,x)=x^{\mathsf{T}}Sx,$$

with ${\nabla }_{x}V$ the gradient in $x$. Adopt the quadratic ansatz $V(t,x)=x^{\mathsf{T}}P(t)x$ with $P(t)=P(t{)}^{\mathsf{T}}$; then ${\nabla }_{x}V=2P(t)x$ and the partial in time is $x^{\mathsf{T}}\dot{P}(t)x$. The bracketed expression to minimise becomes

$$x^{\mathsf{T}}Qx+u^{\mathsf{T}}Ru+2x^{\mathsf{T}}P(Ax+Bu).$$

This is strictly convex in $u$ because $R\succ 0$; setting its $u$-gradient to zero, $2Ru+2B^{\mathsf{T}}Px=0$, yields the minimiser

$$u^{\star }=-R^{-1}{B}^{\mathsf{T}}Px.$$

Substitute $u^{\star }$ back. The control terms collapse: $u^{\star \mathsf{T}}R{u}^{\star }+2x^{\mathsf{T}}PB{u}^{\star }=x^{\mathsf{T}}PB{R}^{-1}{B}^{\mathsf{T}}Px-2x^{\mathsf{T}}PB{R}^{-1}{B}^{\mathsf{T}}Px=-x^{\mathsf{T}}PB{R}^{-1}{B}^{\mathsf{T}}Px$. The HJB equation reads

$$-x^{\mathsf{T}}\dot{P}x=x^{\mathsf{T}}(Q+A^{\mathsf{T}}P+PA-PB{R}^{-1}{B}^{\mathsf{T}}P)x\text{for all }x,$$

where $2x^{\mathsf{T}}PAx=x^{\mathsf{T}}(A^{\mathsf{T}}P+PA)x$ symmetrises the cross term. Matching the symmetric quadratic forms gives exactly the Riccati differential equation $-\dot{P}=A^{\mathsf{T}}P+PA-PB{R}^{-1}{B}^{\mathsf{T}}P+Q$ with terminal data $P(t_f)=S$.

Existence on the whole interval: the right-hand side is a quadratic (hence locally Lipschitz) field in $P$, so a local solution exists backward from $t_f$. Finite escape is excluded because $P(t)$ is sandwiched: $0⪯P(t)⪯P_{\mathrm{lin}}(t)$, where $P_{\mathrm{lin}}$ solves the linear Lyapunov equation obtained by fixing any stabilising feedback and evaluating its (finite) cost-to-go. The lower bound $P(t)⪰0$ holds because $V(t,x)\ge 0$ as the minimum of a non-negative cost; the upper bound holds because the optimal cost cannot exceed the cost of a fixed admissible policy, which is the finite quadratic form $P_{\mathrm{lin}}$. A solution bounded in the Loewner order on a compact interval cannot blow up, so $P(t)$ extends to all of $[t_0,t_f]$.

Optimality and uniqueness: for any admissible $u$ with resulting trajectory $x(\cdot )$, the function $t\mapsto x(t{)}^{\mathsf{T}}P(t)x(t)$ has derivative, using the dynamics and the RDE,

$$\frac{d}{dt}\left(x^{\mathsf{T}}Px\right)=x^{\mathsf{T}}\dot{P}x+2x^{\mathsf{T}}P(Ax+Bu)=-\ell (x,u)+\Vert R^{1/2}(u-u^{\star }){\Vert }^2,$$

after completing the square (the algebra is the same collapse as above, with the leftover $\Vert R^{1/2}(u-u^{\star }){\Vert }^2\ge 0$). Integrating from $t_0$ to $t_f$ and using $P(t_f)=S$,

$$J(u)=x_0^{\mathsf{T}}P(t_0)x_0+{\int }_{t_0}^{t_f}\Vert R^{1/2}(u-u^{\star }){\Vert }^{2}dt\ge x_0^{\mathsf{T}}P(t_0)x_0=V(t_0,x_0),$$

with equality iff $u\equiv u^{\star }=-R^{-1}{B}^{\mathsf{T}}Px$ almost everywhere. So $u^{\star }$ is the unique optimal control and $V(t_0,x_0)=x_0^{\mathsf{T}}P(t_0)x_0$. $\square$

Bridge. This completion-of-squares identity builds toward the infinite-horizon theory of the Advanced results, where letting $t_f\to \infty$ sends the Riccati matrix $P(t)$ to a constant $P$ solving the algebraic Riccati equation, and the same square-completion certifies that the resulting static feedback is optimal over the infinite horizon. The Riccati differential equation is exactly the Bellman recursion of dynamic programming written in continuous time, so this is the foundational reason a single backward sweep yields the optimal policy at every state at once rather than one trajectory at a time. The argument generalises the scalar $p^2=1$ of the Beginner worked example to a matrix quadratic, and the structure of that matrix quadratic is dual to the Lyapunov equation of 43.06.12: dropping the quadratic term $PB{R}^{-1}{B}^{\mathsf{T}}P$ leaves the closed-loop Lyapunov equation $A_{\mathrm{cl}}^{\mathsf{T}}P+P{A}_{\mathrm{cl}}=-(Q+K^{\mathsf{T}}RK)$ that the optimal $P$ also satisfies. The central insight is that quadratic cost plus linear dynamics forces a quadratic value function, collapsing an infinite-dimensional minimisation over control functions into a finite-dimensional matrix equation; putting these together, the bridge is the recognition that the cost-to-go Hessian $P$ is simultaneously the object that prices every future and the object whose feedback $K=R^{-1}{B}^{\mathsf{T}}P$ acts now. The Riccati matrix appears again in the Kalman filter, the estimation dual treated downstream.

Exercises Intermediate+

Advanced results Master

Theorem (existence and uniqueness of the stabilising CARE solution; Anderson-Moore Ch. 3-5 ^{[source pending]}; Lancaster-Rodman Ch. 7-9 ^{[source pending]}). Let $(A,B)$ be stabilisable and $(A,Q^{1/2})$ detectable, with $R\succ 0$, $Q⪰0$. Then the continuous-time algebraic Riccati equation $A^{\mathsf{T}}P+PA-PB{R}^{-1}{B}^{\mathsf{T}}P+Q=0$ has a unique symmetric positive-semidefinite solution $P$ that is stabilising, meaning $A_{\mathrm{cl}}=A-B{R}^{-1}{B}^{\mathsf{T}}P$ is Hurwitz. This $P$ is the maximal symmetric solution in the Loewner order, equals the infinite-horizon value Hessian, and the static feedback $u=-R^{-1}{B}^{\mathsf{T}}Px$ minimises ${\int }_0^{\infty }(x^{\mathsf{T}}Qx+u^{\mathsf{T}}Ru)dt$. If detectability is strengthened to observability of $(A,Q^{1/2})$, then $P\succ 0$.

Stabilisability is what makes the infinite-horizon cost finite for at least one policy, so the value function is well defined; detectability is what forces the optimal closed loop to be stable rather than merely cost-finite, by ruling out an undetected unstable mode that contributes zero running cost. Under both, the finite-horizon Riccati solution $P(t_0)$ converges monotonically as $t_f\to \infty$ to the stabilising $P$, independent of the terminal weight $S⪰0$, which is how the algebraic equation inherits its solution from the differential one.

Theorem (Hamiltonian / invariant-subspace construction; Lancaster-Rodman Ch. 7-9 ^{[source pending]}). Form the Hamiltonian matrix

$$\mathcal{H}=\left(\begin{array}{cc}A& -B{R}^{-1}{B}^{\mathsf{T}}\\ -Q& -A^{\mathsf{T}}\end{array}\right)\in {\mathbb{R}}^{2n\times 2n}.$$

Under stabilisability and detectability $\mathcal{H}$ has no eigenvalues on the imaginary axis, so it has an $n$-dimensional stable invariant subspace; choosing a basis \begin{psmallmatrix} X \\ Y\end{psmallmatrix} for it with $X$ invertible, the stabilising solution is $P=Y{X}^{-1}$, which is symmetric, and the eigenvalues of the optimal closed loop $A_{\mathrm{cl}}$ are exactly the $n$ stable eigenvalues of $\mathcal{H}$. The eigenvalues of $\mathcal{H}$ come in pairs $\{\lambda ,-\lambda \}$ by the Hamiltonian symmetry ${\mathcal{J}}^{-1}{\mathcal{H}}^{\mathsf{T}}\mathcal{J}=-\mathcal{H}$, with \mathcal J = \begin{psmallmatrix} 0 & I \\ -I & 0\end{psmallmatrix}, so the imaginary axis splits the $2n$ eigenvalues into $n$ stable and $n$ unstable.

This reduces solving a quadratic matrix equation to a linear-algebra computation: a real Schur reduction of $\mathcal{H}$ that orders the stable eigenvalues first, after which the leading $n$ Schur columns furnish \begin{psmallmatrix} X \\ Y\end{psmallmatrix} and $P=Y{X}^{-1}$ follows by one triangular solve. This is the Schur method of Laub (1979), built on the same real Schur form of [43.06] that the Bartels-Stewart solver of 43.06.12 uses; it is the production route in LAPACK/SLICOT, numerically preferable to iterating the Riccati equation. The ordered Schur form is what selects the stabilising root from the lattice of symmetric solutions: each Lagrangian $\mathcal{H}$-invariant subspace gives a different symmetric solution, and the stable one alone yields the Hurwitz closed loop.

Theorem (Newton-Kleinman iteration; Anderson-Moore Ch. 3 ^{[source pending]}). Start from any stabilising feedback $K_0$ (so $A-B{K}_0$ Hurwitz). Iterate: given $K_i$, solve the Lyapunov equation $(A-B{K}_i{)}^{\mathsf{T}}{P}_i+P_i(A-B{K}_i)=-(Q+K_i^{\mathsf{T}}R{K}_i)$ for $P_i⪰0$, then set $K_{i+1}=R^{-1}{B}^{\mathsf{T}}{P}_i$. The iterates $P_i$ decrease monotonically, $P_0⪰P_1⪰\cdots ⪰P$, and converge quadratically to the stabilising solution $P$, with every $A-B{K}_i$ Hurwitz along the way.

Each Newton-Kleinman step is exactly one Lyapunov solve of 43.06.12 — a Bartels-Stewart / Hammarling computation — so the entire CARE solver is a short loop over the matrix-equation solver of the numerical-analysis spine. The iteration is Newton's method applied to the Riccati operator $\mathcal{R}(P)=A^{\mathsf{T}}P+PA-PB{R}^{-1}{B}^{\mathsf{T}}P+Q$, whose Fréchet derivative at $P$ in direction $H$ is the Lyapunov operator $H\mapsto A_{\mathrm{cl}}^{\mathsf{T}}H+H{A}_{\mathrm{cl}}$; the Hurwitz closed loop makes that derivative invertible (its spectrum avoids zero by 43.06.12), which is the source of the quadratic local rate.

Theorem (discrete-time regulator and the DARE; Kwakernaak-Sivan Ch. 3 ^{[source pending]}). For $x_{k+1}=A{x}_k+B{u}_k$ with $(A,B)$ stabilisable and $(A,Q^{1/2})$ detectable, the discrete-time algebraic Riccati equation $P=A^{\mathsf{T}}PA-A^{\mathsf{T}}PB(R+B^{\mathsf{T}}PB{)}^{-1}{B}^{\mathsf{T}}PA+Q$ has a unique stabilising symmetric positive-semidefinite solution, the optimal gain is $K=(R+B^{\mathsf{T}}PB{)}^{-1}{B}^{\mathsf{T}}PA$, and the closed loop $A-BK$ has spectral radius below one. The discrete solution arises from the stable deflating subspace of the symplectic matrix pencil \lambda \begin{psmallmatrix} I & BR^{-1}B^{\mathsf T} \\ 0 & A^{\mathsf T}\end{psmallmatrix} - \begin{psmallmatrix} A & 0 \\ -Q & I\end{psmallmatrix}, whose eigenvalues are symmetric under $\lambda \leftrightarrow 1/\bar{\lambda }$ about the unit circle rather than the imaginary axis; the unit circle plays the role the imaginary axis plays in continuous time.

Synthesis. The linear-quadratic regulator is one object — the value function $V=x^{\mathsf{T}}Px$ of a quadratic cost over linear dynamics — and the Riccati matrix $P$ is the central insight that organises every face of it: the differential equation $-\dot{P}=A^{\mathsf{T}}P+PA-PB{R}^{-1}{B}^{\mathsf{T}}P+Q$ is the continuous-time Bellman recursion, its steady state is the algebraic equation, its gain $K=R^{-1}{B}^{\mathsf{T}}P$ is the controller, and its Hamiltonian $\mathcal{H}$ converts the quadratic solve into a linear invariant-subspace computation. This is exactly why dropping the quadratic term recovers the Lyapunov equation of 43.06.12: the optimal $P$ satisfies the closed-loop Lyapunov identity $A_{\mathrm{cl}}^{\mathsf{T}}P+P{A}_{\mathrm{cl}}=-(Q+K^{\mathsf{T}}RK)$, so $V$ is simultaneously the optimal cost-to-go and a Lyapunov function certifying that $A_{\mathrm{cl}}$ is Hurwitz, and the foundational reason LQR guarantees stability is that optimality and the Lyapunov stability test of 02.12.08 are the same equation viewed from two sides.

The construction generalises the scalar quadratic $p^2=2ap+Q$ to a matrix quadratic whose stabilising root is selected by the stable invariant subspace of $\mathcal{H}$, and that selection is dual to the eigenvalue split: the $n$ left-half-plane eigenvalues of $\mathcal{H}$ become the closed-loop poles. Putting these together, the entire solver pipeline rides on the numerical-linear-algebra spine — the Schur form of [43.06] orders the Hamiltonian spectrum, and each Newton-Kleinman step is a Bartels-Stewart Lyapunov solve of 43.06.12 — so the bridge from optimal control to computation is the recognition that a quadratic matrix equation is an eigenvalue problem in disguise, and its conditioning is governed by the same separation that conditions the Lyapunov solve underneath it.

Full proof set Master

Proposition 1 (the HJB ansatz reduces LQR to the Riccati equation). The value function of the finite-horizon LQR is $V(t,x)=x^{\mathsf{T}}P(t)x$ with $P$ solving the RDE $-\dot{P}=A^{\mathsf{T}}P+PA-PB{R}^{-1}{B}^{\mathsf{T}}P+Q$, $P(t_f)=S$, and the optimal control is $u^{\star }=-R^{-1}{B}^{\mathsf{T}}Px$.

Proof. This is the Key theorem. The HJB equation $-{\partial }_{t}V={\min }_u[\ell (x,u)+{\nabla }_{x}V\cdot (Ax+Bu)]$ with $V=x^{\mathsf{T}}Px$ has inner minimiser $u^{\star }=-R^{-1}{B}^{\mathsf{T}}Px$ from the first-order condition $2Ru+2B^{\mathsf{T}}Px=0$ and the Hessian $2R\succ 0$; substituting collapses the control terms to $-x^{\mathsf{T}}PB{R}^{-1}{B}^{\mathsf{T}}Px$ and matching symmetric quadratic forms yields the RDE. The completion-of-squares identity $\frac{d}{dt}(x^{\mathsf{T}}Px)=-\ell (x,u)+\Vert R^{1/2}(u-u^{\star }){\Vert }^2$ integrates to $J(u)=x_0^{\mathsf{T}}P(t_0)x_0+\int \Vert R^{1/2}(u-u^{\star }){\Vert }^{2}dt$, proving optimality and uniqueness of $u^{\star }$. $\square$

Proposition 2 (the optimal $P$ solves a closed-loop Lyapunov equation). If $P$ solves the CARE and $K=R^{-1}{B}^{\mathsf{T}}P$, $A_{\mathrm{cl}}=A-BK$, then $A_{\mathrm{cl}}^{\mathsf{T}}P+P{A}_{\mathrm{cl}}=-(Q+K^{\mathsf{T}}RK)$.

Proof. Expand $A_{\mathrm{cl}}^{\mathsf{T}}P+P{A}_{\mathrm{cl}}=A^{\mathsf{T}}P+PA-K^{\mathsf{T}}{B}^{\mathsf{T}}P-PBK$. With $K=R^{-1}{B}^{\mathsf{T}}P$ and $P=P^{\mathsf{T}}$, $R=R^{\mathsf{T}}$, both $PBK$ and $K^{\mathsf{T}}{B}^{\mathsf{T}}P$ equal $PB{R}^{-1}{B}^{\mathsf{T}}P$, so the sum is $A^{\mathsf{T}}P+PA-2PB{R}^{-1}{B}^{\mathsf{T}}P$. The CARE gives $A^{\mathsf{T}}P+PA=PB{R}^{-1}{B}^{\mathsf{T}}P-Q$, hence the sum equals $-Q-PB{R}^{-1}{B}^{\mathsf{T}}P$. Since $K^{\mathsf{T}}RK=PB{R}^{-1}{B}^{\mathsf{T}}P$, this is $-(Q+K^{\mathsf{T}}RK)$. $\square$

Proposition 3 (stabilising CARE solution implies a Hurwitz closed loop, given detectability). Let $P⪰0$ solve the CARE with $(A,Q^{1/2})$ detectable, and suppose $V(x)=x^{\mathsf{T}}Px$ is non-increasing along $\dot{x}=A_{\mathrm{cl}}x$. Then $A_{\mathrm{cl}}$ is Hurwitz.

Proof. By Proposition 2, $\dot{V}=-x^{\mathsf{T}}(Q+K^{\mathsf{T}}RK)x\le 0$, so $V$ is a Lyapunov function and the closed loop is at least marginally stable. Suppose $A_{\mathrm{cl}}$ had an eigenvalue $\lambda$ with $\mathrm{Re}\lambda \ge 0$ and eigenvector $v$. Along the mode $x(t)=e^{\lambda t}v$ the LaSalle argument forces $\dot{V}\equiv 0$, hence $Q^{1/2}x\equiv 0$ and $Kx\equiv 0$ on this trajectory; with $Kx\equiv 0$ the closed-loop dynamics reduce to $\dot{x}=Ax$, so $(A,Q^{1/2})$ has an unobservable mode at $\lambda$ with $\mathrm{Re}\lambda \ge 0$, contradicting detectability. Therefore every eigenvalue of $A_{\mathrm{cl}}$ has strictly negative real part. $\square$

Proposition 4 (Hamiltonian invariant subspace yields the symmetric solution). If \begin{psmallmatrix} X \\ Y\end{psmallmatrix} spans an $\mathcal{H}$-invariant subspace with $X$ invertible and the subspace is Lagrangian (i.e. $X^{\mathsf{T}}Y=Y^{\mathsf{T}}X$), then $P=Y{X}^{-1}$ is symmetric and solves the CARE.

Proof. Invariance gives \mathcal H\begin{psmallmatrix} X \\ Y\end{psmallmatrix} = \begin{psmallmatrix} X \\ Y\end{psmallmatrix}M for some $M\in {\mathbb{R}}^{n\times n}$. The top block reads $AX-B{R}^{-1}{B}^{\mathsf{T}}Y=XM$, so $M=X^{-1}(AX-B{R}^{-1}{B}^{\mathsf{T}}Y)$. The bottom block reads $-QX-A^{\mathsf{T}}Y=YM=Y{X}^{-1}(AX-B{R}^{-1}{B}^{\mathsf{T}}Y)$. Right-multiply by $X^{-1}$ and write $P=Y{X}^{-1}$: $-Q-A^{\mathsf{T}}P=P(A-B{R}^{-1}{B}^{\mathsf{T}}P)$, which rearranges to the CARE $A^{\mathsf{T}}P+PA-PB{R}^{-1}{B}^{\mathsf{T}}P+Q=0$. Symmetry $P=P^{\mathsf{T}}$ is the Lagrangian condition $X^{\mathsf{T}}Y=Y^{\mathsf{T}}X$: $X^{\mathsf{T}}(P-P^{\mathsf{T}})X=X^{\mathsf{T}}Y{X}^{-1}X-X^{\mathsf{T}}{X}^{-\mathsf{T}}{Y}^{\mathsf{T}}X=X^{\mathsf{T}}Y-Y^{\mathsf{T}}X=0$, and $X$ invertible forces $P=P^{\mathsf{T}}$. Selecting the stable invariant subspace makes $M$ similar to $A_{\mathrm{cl}}$ Hurwitz, so $P$ is the stabilising solution. $\square$

Proposition 5 (Hamiltonian eigenvalue symmetry). The spectrum of $\mathcal{H}$ is symmetric about the imaginary axis: $\lambda \in \sigma (\mathcal{H})⟺-\lambda \in \sigma (\mathcal{H})$.

Proof. With \mathcal J = \begin{psmallmatrix} 0 & I \\ -I & 0\end{psmallmatrix}, a direct block computation gives ${\mathcal{J}}^{-1}{\mathcal{H}}^{\mathsf{T}}\mathcal{J}=-\mathcal{H}$ (the defining property of a Hamiltonian matrix; here ${\mathcal{J}}^{-1}=-\mathcal{J}={\mathcal{J}}^{\mathsf{T}}$). Thus ${\mathcal{H}}^{\mathsf{T}}$ and $-\mathcal{H}$ are similar, so they share a spectrum. Since ${\mathcal{H}}^{\mathsf{T}}$ and $\mathcal{H}$ always share a spectrum, $\sigma (\mathcal{H})=\sigma (-\mathcal{H})=-\sigma (\mathcal{H})$, the claimed symmetry. With no eigenvalues on the imaginary axis (guaranteed under stabilisability and detectability), the $2n$ eigenvalues split into $n$ in the open left half-plane and $n$ in the open right half-plane, and the stable invariant subspace of Proposition 4 is well defined. $\square$

Connections Master

The Riccati equation is the quadratic enrichment of the Lyapunov equation of 43.06.12: deleting its quadratic term $PB{R}^{-1}{B}^{\mathsf{T}}P$ leaves the closed-loop Lyapunov equation $A_{\mathrm{cl}}^{\mathsf{T}}P+P{A}_{\mathrm{cl}}=-(Q+K^{\mathsf{T}}RK)$ that the optimal cost-to-go also satisfies, and every Newton-Kleinman step solves exactly one such Lyapunov equation by the Bartels-Stewart / Hammarling algorithm of that unit, so the conditioning of the regulator solve is governed by the separation $\mathrm{sep}(A_{\mathrm{cl}},A_{\mathrm{cl}}^{\mathsf{T}})$ that conditions the underlying matrix-equation solve.

The LQR closed-loop guarantee is the algebraic face of the Lyapunov direct method of 02.12.08: the value function $V(x)=x^{\mathsf{T}}Px$ built from the stabilising CARE solution is a quadratic Lyapunov function whose derivative along $\dot{x}=A_{\mathrm{cl}}x$ is $-x^{\mathsf{T}}(Q+K^{\mathsf{T}}RK)x\le 0$, so optimality and asymptotic stability are the same statement, and detectability of $(A,Q^{1/2})$ is exactly the observability condition that upgrades the non-strict decrease to strict asymptotic stability.

The flow $\dot{x}=A_{\mathrm{cl}}x$ and the cost integrals rest on the matrix-exponential theory of 02.06.03: the infinite-horizon cost of any stabilising feedback is the convergent Gramian ${\int }_0^{\infty }{e}^{A_{\mathrm{cl}}^{\mathsf{T}}t}(Q+K^{\mathsf{T}}RK)e^{A_{\mathrm{cl}}t}dt$, the same exponential representation that solves the Lyapunov equation, and the Riccati matrix is the minimum of these costs over all stabilising feedbacks.

The positive-semidefinite structure of $Q$, $R$, and $P$ is the spectral/PSD theory of 01.01.13: the weights live in the positive cone, the CARE solution is selected to be the maximal positive-semidefinite root, and the principal-axes decomposition of $P$ is what diagonalises the value function $x^{\mathsf{T}}Px$ into independent quadratic costs along the eigendirections, the geometric content of the cost-to-go.

The Hamiltonian construction forward-links to the Pontryagin maximum principle of this chapter (44.07.01) and to the dynamic-programming / HJB framework of [44.08]: the Hamiltonian matrix $\mathcal{H}$ is the linear-quadratic specialisation of the Pontryagin costate dynamics, and the value function $V=x^{\mathsf{T}}Px$ is the linear-quadratic instance of the HJB value function, so LQR is simultaneously the worked example that makes both general theories concrete and the case where they are solvable in closed form.

The estimation dual forward-links to the Kalman filter (a downstream statistics/estimation unit): the filter's error-covariance recursion is a Riccati equation with $(A,B,Q,R)$ replaced by the transposed/observation data, so the CARE solver and its Hamiltonian-subspace machinery transfer verbatim to optimal linear estimation by duality.

Historical & philosophical context Master

The Riccati equation descends from the scalar nonlinear differential equation $y^{\prime }=q_0(x)+q_1(x)y+q_2(x)y^2$ studied by Jacopo Riccati in the 1720s, whose defining feature — a quadratic nonlinearity linearisable by a change to a second-order linear equation — is the same structure that lets the matrix Riccati equation be solved through the linear Hamiltonian system ^{[Riccati 1724]}. The control-theoretic matrix version emerged from Rudolf Kalman's 1960 reformulation of optimal control and filtering in state-space terms; his work on the linear-quadratic regulator and, with Richard Bucy, on the dual filtering problem established the Riccati equation as the central computational object of modern control ^{[Kalman 1960]}.

The dynamic-programming derivation used here is Richard Bellman's, whose principle of optimality and the resulting Hamilton-Jacobi-Bellman equation reduce the search over control functions to a backward recursion on the value function ^{[Bellman 1957]}. The existence-uniqueness theory under stabilisability and detectability, the maximal-solution characterisation, and the closed-loop stability guarantee were systematised by Wonham, Kucera, and Willems in the 1960s and 1970s, and consolidated in the monographs of Anderson and Moore and of Kwakernaak and Sivan ^{[Willems 1971]}. The numerically stable Hamiltonian-Schur solver was given by Alan Laub in 1979, replacing direct iteration of the Riccati equation with an ordered real Schur reduction; it and the Newton-Kleinman iteration are the algorithms behind the Riccati routines of LAPACK, SLICOT, and the control toolboxes.

Bibliography Master

@book{Bertsekas2017DP,
  author    = {Bertsekas, Dimitri P.},
  title     = {Dynamic Programming and Optimal Control, Volume I},
  edition   = {4th},
  publisher = {Athena Scientific},
  address   = {Belmont, MA},
  year      = {2017}
}

@book{AndersonMoore1990,
  author    = {Anderson, Brian D. O. and Moore, John B.},
  title     = {Optimal Control: Linear Quadratic Methods},
  publisher = {Prentice-Hall},
  address   = {Englewood Cliffs, NJ},
  year      = {1990}
}

@book{KwakernaakSivan1972,
  author    = {Kwakernaak, Huibert and Sivan, Raphael},
  title     = {Linear Optimal Control Systems},
  publisher = {Wiley-Interscience},
  address   = {New York},
  year      = {1972}
}

@book{LancasterRodman1995,
  author    = {Lancaster, Peter and Rodman, Leiba},
  title     = {Algebraic Riccati Equations},
  publisher = {Oxford University Press},
  address   = {Oxford},
  year      = {1995}
}

@article{Kalman1960,
  author  = {Kalman, Rudolf E.},
  title   = {Contributions to the Theory of Optimal Control},
  journal = {Bolet\'in de la Sociedad Matem\'atica Mexicana},
  volume  = {5},
  year    = {1960},
  pages   = {102--119}
}

@book{Bellman1957,
  author    = {Bellman, Richard E.},
  title     = {Dynamic Programming},
  publisher = {Princeton University Press},
  address   = {Princeton, NJ},
  year      = {1957}
}

@article{Willems1971,
  author  = {Willems, Jan C.},
  title   = {Least Squares Stationary Optimal Control and the Algebraic Riccati Equation},
  journal = {IEEE Transactions on Automatic Control},
  volume  = {16},
  number  = {6},
  year    = {1971},
  pages   = {621--634}
}

@article{Laub1979,
  author  = {Laub, Alan J.},
  title   = {A Schur Method for Solving Algebraic Riccati Equations},
  journal = {IEEE Transactions on Automatic Control},
  volume  = {24},
  number  = {6},
  year    = {1979},
  pages   = {913--921}
}

@article{Kleinman1968,
  author  = {Kleinman, David L.},
  title   = {On an Iterative Technique for Riccati Equation Computations},
  journal = {IEEE Transactions on Automatic Control},
  volume  = {13},
  number  = {1},
  year    = {1968},
  pages   = {114--115}
}