In previous state-feedback lessons, the controller was usually selected by assigning desired poles or by augmenting the system for tracking. In this chapter, the emphasis shifted from “place the eigenvalues here” to “measure performance through a quadratic form.” The natural question is: can the gain \( K \) be chosen by minimizing a precise state-and-input performance index?

For a continuous-time LTI plant \( \dot{\mathbf{x} }=A\mathbf{x}+B\mathbf{u} \), the simplest optimal state-feedback problem penalizes state deviation and control effort:

\[ J(\mathbf{x}_0,\mathbf{u})=\int_0^\infty \left(\mathbf{x}(t)^TQ\mathbf{x}(t)+ \mathbf{u}(t)^TR\mathbf{u}(t)\right)dt,\quad Q=Q^T\succeq 0, \quad R=R^T\succ 0. \]

The preview result is that, under the standard stabilizability and detectability assumptions, the minimizing controller has the linear state-feedback form \( \mathbf{u}^*(t)=-K\mathbf{x}(t) \). Thus the state-feedback structure introduced in Chapter 22 is not only a design template; it is also the solution of a mathematically well-defined optimization problem.

flowchart TD
  A["Choose performance: state error and input effort"] --> B["Select Q and R"]
  B --> C["Solve Riccati equation for P"]
  C --> D["Compute K = inv(R) B^T P"]
  D --> E["Close loop: Acl = A - B K"]
  E --> F["Check stability and cost"]
  F --> G["Interpret trade-off: speed vs effort"]
        

The finite-horizon regulator includes a terminal penalty and is the cleanest bridge to the optimal control course. Given \( t_0=0 \) and terminal time \( T \), define

\[ J_T=\mathbf{x}(T)^TS\mathbf{x}(T)+ \int_0^T\left(\mathbf{x}(t)^TQ\mathbf{x}(t)+ \mathbf{u}(t)^TR\mathbf{u}(t)\right)dt,\quad S=S^T\succeq 0. \]

In an optimal control course, this problem is usually derived by the Hamilton-Jacobi-Bellman equation or by the calculus of variations. For a linear plant and quadratic cost, the value function is assumed to be quadratic:

\[ V(t,\mathbf{x})=\mathbf{x}^TP(t)\mathbf{x},\quad P(t)=P(t)^T. \]

Substitution into the optimality condition gives the time-varying optimal feedback

\[ \mathbf{u}^*(t)=-R^{-1}B^TP(t)\mathbf{x}(t), \]

where \( P(t) \) satisfies the differential Riccati equation integrated backward from \( P(T)=S \):

\[ -\dot{P}(t)=A^TP(t)+P(t)A-P(t)BR^{-1}B^TP(t)+Q,\quad P(T)=S. \]

The finite-horizon gain is therefore time-varying: \( K(t)=R^{-1}B^TP(t) \). This is an important conceptual difference from ordinary pole placement, where a constant gain is usually selected directly.

A compact proof of the LQR feedback law can be obtained by differentiating \( \mathbf{x}^TP(t)\mathbf{x} \) along the trajectory:

\[ \frac{d}{dt}\left(\mathbf{x}^TP\mathbf{x}\right)= \mathbf{x}^T\dot{P}\mathbf{x}+ \mathbf{x}^T(A^TP+PA)\mathbf{x}+ 2\mathbf{x}^TPB\mathbf{u}. \]

Add and subtract the quadratic term that completes the square:

\[ \mathbf{u}^TR\mathbf{u}+ 2\mathbf{x}^TPB\mathbf{u}= \left(\mathbf{u}+R^{-1}B^TP\mathbf{x}\right)^TR \left(\mathbf{u}+R^{-1}B^TP\mathbf{x}\right) -\mathbf{x}^TPBR^{-1}B^TP\mathbf{x}. \]

If \( P(t) \) satisfies the differential Riccati equation, then the finite-horizon cost becomes

\[ J_T=\mathbf{x}(0)^TP(0)\mathbf{x}(0)+ \int_0^T \left(\mathbf{u}+R^{-1}B^TP\mathbf{x}\right)^TR \left(\mathbf{u}+R^{-1}B^TP\mathbf{x}\right)dt. \]

Because \( R\succ 0 \), the integral is minimized exactly when the squared term vanishes. Hence \( \mathbf{u}^*=-R^{-1}B^TP\mathbf{x} \), and the optimal cost equals \( \mathbf{x}(0)^TP(0)\mathbf{x}(0) \).

For a time-invariant plant and an infinite horizon, the steady-state solution of the Riccati differential equation is a constant symmetric matrix \( P=P^T \). Setting \( \dot{P}=0 \) gives the continuous-time algebraic Riccati equation:

\[ A^TP+PA-PBR^{-1}B^TP+Q=0. \]

The corresponding gain and closed-loop matrix are

\[ K=R^{-1}B^TP,\quad A_{cl}=A-BK. \]

The stabilizing Riccati solution is the solution for which \( A-BR^{-1}B^TP \) is Hurwitz. When \( (A,B) \) is stabilizable and \( (Q^{1/2},A) \) is detectable, there exists a unique stabilizing positive semidefinite solution \( P\succeq 0 \), and the feedback above minimizes the infinite-horizon cost.

\[ J_\infty^*(\mathbf{x}_0)=\mathbf{x}_0^TP\mathbf{x}_0. \]

This expression gives an interpretation of \( P \): it is not merely an algebraic artifact. It is the matrix representing the optimal cost-to-go as a quadratic form in the initial state.

Let \( \mathbf{u}=-K\mathbf{x} \), where \( K=R^{-1}B^TP \). Substitute the Riccati equation into the closed-loop Lyapunov derivative:

\[ A_{cl}^TP+PA_{cl} =A^TP+PA-PBR^{-1}B^TP-PBR^{-1}B^TP. \]

Since the Riccati equation implies \( A^TP+PA-PBR^{-1}B^TP=-Q \), we obtain

\[ A_{cl}^TP+PA_{cl}=-(Q+K^TRK). \]

Therefore, for \( V(\mathbf{x})=\mathbf{x}^TP\mathbf{x} \),

\[ \dot{V}(\mathbf{x})= -\mathbf{x}^T(Q+K^TRK)\mathbf{x}\leq 0. \]

This proves nonincreasing optimal cost along closed-loop trajectories. Under the usual detectability condition, the only trajectory that can keep \( \dot{V}=0 \) indefinitely is the zero trajectory, which gives asymptotic stability of \( A_{cl} \). Thus the same matrix \( P \) simultaneously encodes cost and proves stability.

flowchart TD
  X["state x"] --> U["u = -K x"]
  U --> PLANT["plant: x_dot = A x + B u"]
  PLANT --> X
  QW["state penalty Q"] --> RIC["Riccati equation"]
  RW["input penalty R"] --> RIC
  RIC --> P["cost matrix P"]
  P --> K["gain K"]
  K --> U
        

Pole placement chooses eigenvalues first and then constructs a feedback gain. LQR chooses a quadratic performance index first and then lets the optimality conditions determine both the gain and the closed-loop poles. These approaches answer different design questions:

\[ \text{pole placement: desired poles } \Longrightarrow K, \qquad \text{LQR: } (Q,R) \Longrightarrow P \Longrightarrow K \Longrightarrow \lambda(A-BK). \]

Increasing \( Q \) relative to \( R \) usually makes state deviation more expensive, leading to faster regulation and larger control effort. Increasing \( R \) makes input more expensive, often producing slower but less aggressive motion. Unlike arbitrary pole placement, LQR tends to avoid unnecessarily large gains because input effort is part of the objective.

A useful scalar scaling experiment is \( Q_\alpha=\alpha Q \) with \( R \) fixed. For a controllable system, large \( \alpha \) typically moves dominant closed-loop poles farther left, but the precise movement is constrained by the system matrices and the Riccati equation rather than by direct pole assignment.

The algebraic Riccati equation is nonlinear in \( P \), but it can be connected to a linear invariant-subspace problem. Define the Hamiltonian matrix

\[ \mathcal{H}= \begin{bmatrix} A & -BR^{-1}B^T \\ -Q & -A^T \end{bmatrix}. \]

If the stable invariant subspace of \( \mathcal{H} \) is spanned by columns of \( \begin{bmatrix}X_1^T & X_2^T\end{bmatrix}^T \) and \( X_1 \) is nonsingular, then

\[ P=X_2X_1^{-1}. \]

This viewpoint explains why numerically reliable Riccati solvers often use Schur decompositions or ordered invariant subspaces rather than directly solving the quadratic matrix equation coefficient by coefficient.

For a sampled-data system \( \mathbf{x}_{k+1}=A_d\mathbf{x}_k+B_d\mathbf{u}_k \), the infinite-horizon quadratic cost is

\[ J_d=\sum_{k=0}^\infty \left(\mathbf{x}_k^TQ\mathbf{x}_k+ \mathbf{u}_k^TR\mathbf{u}_k\right). \]

The discrete-time algebraic Riccati equation is

\[ P=A_d^TPA_d-A_d^TPB_d(R+B_d^TPB_d)^{-1}B_d^TPA_d+Q. \]

The gain is

\[ K_d=(R+B_d^TPB_d)^{-1}B_d^TPA_d,\quad \mathbf{u}_k^*=-K_d\mathbf{x}_k. \]

This preview is important because many modern controllers are deployed digitally. The conceptual structure is the same: quadratic value function, Riccati equation, stabilizing feedback, and cost-to-go matrix.

Common libraries for this topic include Python packages scipy.linalg, python-control, and slycot; C++ libraries such as Eigen plus custom Riccati or Lyapunov routines; Java linear algebra packages such as EJML, ojAlgo, or Apache Commons Math; MATLAB Control System Toolbox functions lqr and care; and Wolfram Mathematica symbolic or numerical matrix solvers.

Chapter28_Lesson5.py

"""
Chapter28_Lesson5.py

Preview of optimal state-feedback for the continuous-time LQR problem.

System:
    x_dot = A x + B u

Cost:
    J = integral_0^infty (x.T Q x + u.T R u) dt

The script solves the continuous-time algebraic Riccati equation (CARE),
computes the feedback gain K = R^{-1} B.T P, simulates the closed-loop
system, and estimates the finite-horizon approximation of the cost.
"""

import numpy as np
from numpy.linalg import inv, eigvals
from scipy.linalg import solve_continuous_are
from scipy.integrate import solve_ivp


def lqr_gain(A: np.ndarray, B: np.ndarray, Q: np.ndarray, R: np.ndarray):
    """Return P and K for the continuous-time LQR problem."""
    P = solve_continuous_are(A, B, Q, R)
    K = inv(R) @ B.T @ P
    return P, K


def simulate_closed_loop(A, B, K, Q, R, x0, tf=8.0, npts=1200):
    """Simulate x_dot = (A - B K) x and estimate the quadratic cost."""
    Acl = A - B @ K

    def rhs(t, x):
        return Acl @ x

    t_eval = np.linspace(0.0, tf, npts)
    sol = solve_ivp(rhs, (0.0, tf), x0, t_eval=t_eval, rtol=1e-9, atol=1e-11)

    X = sol.y.T
    U = -(K @ X.T).T
    running_cost = np.einsum("bi,ij,bj->b", X, Q, X) + np.einsum(
        "bi,ij,bj->b", U, R, U
    )
    J_est = np.trapezoid(running_cost, sol.t)
    return sol.t, X, U, J_est


def main():
    # A standard mass-normalized double integrator.
    A = np.array([[0.0, 1.0],
                  [0.0, 0.0]])
    B = np.array([[0.0],
                  [1.0]])

    # Position is penalized more than velocity; control is relatively cheap.
    Q = np.diag([10.0, 1.0])
    R = np.array([[0.2]])

    P, K = lqr_gain(A, B, Q, R)
    Acl = A - B @ K

    x0 = np.array([1.0, 0.0])
    t, X, U, J_est = simulate_closed_loop(A, B, K, Q, R, x0)

    print("P =")
    print(P)
    print("\nK =")
    print(K)
    print("\nclosed-loop eigenvalues =")
    print(eigvals(Acl))
    print("\nfinite-horizon estimated cost over [0, 8] =", J_est)
    print("infinite-horizon value x0.T P x0 =", x0 @ P @ x0)

    try:
        import matplotlib.pyplot as plt

        plt.figure()
        plt.plot(t, X[:, 0], label="position")
        plt.plot(t, X[:, 1], label="velocity")
        plt.plot(t, U[:, 0], label="control")
        plt.xlabel("time")
        plt.grid(True)
        plt.legend()
        plt.title("Chapter 28 Lesson 5: LQR closed-loop response")
        plt.show()
    except Exception as exc:
        print("Plot skipped:", exc)


if __name__ == "__main__":
    main()
      

The Python implementation uses scipy.linalg.solve_continuous_are. The output gain is computed as \( K=R^{-1}B^TP \), and the script verifies the value identity \( J^*(\mathbf{x}_0)=\mathbf{x}_0^TP\mathbf{x}_0 \) numerically.

Chapter28_Lesson5.cpp

/*
Chapter28_Lesson5.cpp

Continuous-time LQR preview using Eigen.

This implementation demonstrates Kleinman's policy-iteration method for a
stabilizing continuous-time algebraic Riccati equation (CARE) solution.

Compile example:
    g++ -std=c++17 Chapter28_Lesson5.cpp -I /path/to/eigen -O2 -o Chapter28_Lesson5
*/

#include <Eigen/Dense>
#include <iostream>
#include <iomanip>

using Eigen::MatrixXd;
using Eigen::VectorXd;

MatrixXd kron(const MatrixXd& A, const MatrixXd& B) {
    MatrixXd K(A.rows() * B.rows(), A.cols() * B.cols());
    for (int i = 0; i < A.rows(); ++i) {
        for (int j = 0; j < A.cols(); ++j) {
            K.block(i * B.rows(), j * B.cols(), B.rows(), B.cols()) = A(i, j) * B;
        }
    }
    return K;
}

VectorXd vec(const MatrixXd& M) {
    VectorXd v(M.rows() * M.cols());
    int k = 0;
    for (int j = 0; j < M.cols(); ++j) {
        for (int i = 0; i < M.rows(); ++i) {
            v(k++) = M(i, j);
        }
    }
    return v;
}

MatrixXd unvec(const VectorXd& v, int rows, int cols) {
    MatrixXd M(rows, cols);
    int k = 0;
    for (int j = 0; j < cols; ++j) {
        for (int i = 0; i < rows; ++i) {
            M(i, j) = v(k++);
        }
    }
    return M;
}

MatrixXd solveContinuousLyapunov(const MatrixXd& A, const MatrixXd& C) {
    // Solves A^T P + P A + C = 0.
    const int n = A.rows();
    MatrixXd I = MatrixXd::Identity(n, n);
    MatrixXd L = kron(I, A.transpose()) + kron(A.transpose(), I);
    VectorXd rhs = -vec(C);
    VectorXd sol = L.fullPivLu().solve(rhs);
    MatrixXd P = unvec(sol, n, n);
    return 0.5 * (P + P.transpose());
}

MatrixXd rk4Step(const MatrixXd& Acl, const MatrixXd& x, double h) {
    MatrixXd k1 = Acl * x;
    MatrixXd k2 = Acl * (x + 0.5 * h * k1);
    MatrixXd k3 = Acl * (x + 0.5 * h * k2);
    MatrixXd k4 = Acl * (x + h * k3);
    return x + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
}

int main() {
    MatrixXd A(2, 2);
    A << 0.0, 1.0,
         0.0, 0.0;

    MatrixXd B(2, 1);
    B << 0.0,
         1.0;

    MatrixXd Q = MatrixXd::Zero(2, 2);
    Q(0, 0) = 10.0;
    Q(1, 1) = 1.0;

    MatrixXd R(1, 1);
    R << 0.2;

    // Initial stabilizing gain for the double integrator.
    MatrixXd K(1, 2);
    K << 2.0, 2.0;

    MatrixXd Rinv = R.inverse();
    MatrixXd P = MatrixXd::Zero(2, 2);

    for (int iter = 0; iter < 50; ++iter) {
        MatrixXd Acl = A - B * K;
        MatrixXd C = Q + K.transpose() * R * K;
        MatrixXd Pnew = solveContinuousLyapunov(Acl, C);
        MatrixXd Knew = Rinv * B.transpose() * Pnew;

        double err = (Knew - K).norm();
        K = Knew;
        P = Pnew;

        if (err < 1e-11) {
            break;
        }
    }

    MatrixXd Acl = A - B * K;
    MatrixXd x(2, 1);
    x << 1.0,
         0.0;

    double h = 0.002;
    int steps = static_cast<int>(8.0 / h);
    double J = 0.0;

    for (int i = 0; i < steps; ++i) {
        MatrixXd u = -K * x;
        double integrand = (x.transpose() * Q * x)(0, 0) + (u.transpose() * R * u)(0, 0);
        J += integrand * h;
        x = rk4Step(Acl, x, h);
    }

    std::cout << std::setprecision(8);
    std::cout << "P =\n" << P << "\n\n";
    std::cout << "K =\n" << K << "\n\n";
    std::cout << "finite-horizon estimated cost over [0, 8] = " << J << "\n";
    std::cout << "infinite-horizon value x0^T P x0 = " << P(0, 0) << "\n";

    return 0;
}
      

The C++ implementation uses Eigen and illustrates a computational principle: Riccati equations can be solved by repeated Lyapunov equations when an initial stabilizing gain is available. This is the basic idea behind Kleinman iteration.

Chapter28_Lesson5.java

/*
Chapter28_Lesson5.java

Analytic continuous-time LQR preview for the double integrator.

No external library is required. The program uses the special closed-form
CARE solution for:
    A = [[0, 1], [0, 0]]
    B = [[0], [1]]
    Q = diag(q1, q2)
    R = r
*/

public class Chapter28_Lesson5 {
    static double[] dx(double[][] Acl, double[] x) {
        return new double[] {
            Acl[0][0] * x[0] + Acl[0][1] * x[1],
            Acl[1][0] * x[0] + Acl[1][1] * x[1]
        };
    }

    static double[] rk4(double[][] Acl, double[] x, double h) {
        double[] k1 = dx(Acl, x);
        double[] x2 = {x[0] + 0.5 * h * k1[0], x[1] + 0.5 * h * k1[1]};
        double[] k2 = dx(Acl, x2);
        double[] x3 = {x[0] + 0.5 * h * k2[0], x[1] + 0.5 * h * k2[1]};
        double[] k3 = dx(Acl, x3);
        double[] x4 = {x[0] + h * k3[0], x[1] + h * k3[1]};
        double[] k4 = dx(Acl, x4);

        return new double[] {
            x[0] + h * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
            x[1] + h * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0
        };
    }

    public static void main(String[] args) {
        double q1 = 10.0;
        double q2 = 1.0;
        double r = 0.2;

        // CARE for double integrator:
        // p12 = sqrt(q1*r)
        // p22 = sqrt(r*(2*p12 + q2))
        // p11 = p12*p22/r
        double p12 = Math.sqrt(q1 * r);
        double p22 = Math.sqrt(r * (2.0 * p12 + q2));
        double p11 = p12 * p22 / r;

        double k1 = p12 / r;
        double k2 = p22 / r;

        double[][] Acl = {
            {0.0, 1.0},
            {-k1, -k2}
        };

        double[] x = {1.0, 0.0};
        double h = 0.002;
        int steps = (int)(8.0 / h);
        double J = 0.0;

        for (int i = 0; i < steps; i++) {
            double u = -(k1 * x[0] + k2 * x[1]);
            double integrand = q1 * x[0] * x[0] + q2 * x[1] * x[1] + r * u * u;
            J += integrand * h;
            x = rk4(Acl, x, h);
        }

        System.out.printf("P = [[%.8f, %.8f], [%.8f, %.8f]]%n", p11, p12, p12, p22);
        System.out.printf("K = [%.8f, %.8f]%n", k1, k2);
        System.out.printf("finite-horizon estimated cost over [0, 8] = %.8f%n", J);
        System.out.printf("infinite-horizon value x0^T P x0 = %.8f%n", p11);
    }
}
      

The Java implementation intentionally avoids external dependencies by exploiting the analytic Riccati solution for the double integrator. In larger Java projects, use a tested numerical linear algebra library instead of hand-writing general eigenvalue or Schur routines.

Chapter28_Lesson5.m

% Chapter28_Lesson5.m
%
% Preview of optimal state-feedback using continuous-time LQR.
% Requires Control System Toolbox for lqr/care. The optional Simulink part
% creates a closed-loop State-Space block if Simulink is available.

clear; clc; close all;

A = [0 1; 0 0];
B = [0; 1];
Q = diag([10 1]);
R = 0.2;

[K, P, closedLoopPoles] = lqr(A, B, Q, R);
Acl = A - B*K;

fprintf('P =\n');
disp(P);
fprintf('K =\n');
disp(K);
fprintf('closed-loop poles =\n');
disp(closedLoopPoles);

x0 = [1; 0];
tf = 8;
[t, x] = ode45(@(t, x) Acl*x, [0 tf], x0);
u = -(K*x')';

runningCost = sum((x*Q).*x, 2) + (u.^2)*R;
J_est = trapz(t, runningCost);

fprintf('finite-horizon estimated cost over [0, 8] = %.8f\n', J_est);
fprintf('infinite-horizon value x0''*P*x0 = %.8f\n', x0'*P*x0);

figure;
plot(t, x(:,1), 'LineWidth', 1.5); hold on;
plot(t, x(:,2), 'LineWidth', 1.5);
plot(t, u, 'LineWidth', 1.5);
grid on;
xlabel('time');
legend('position', 'velocity', 'control');
title('Chapter 28 Lesson 5: LQR closed-loop response');

% Optional Simulink generation: closed-loop model x_dot = (A-BK)x.
if exist('new_system', 'file') == 2
    mdl = 'Chapter28_Lesson5_Simulink';
    if bdIsLoaded(mdl)
        close_system(mdl, 0);
    end
    new_system(mdl);
    open_system(mdl);

    add_block('simulink/Continuous/State-Space', [mdl '/ClosedLoopStateSpace'], ...
        'A', mat2str(Acl), ...
        'B', mat2str(zeros(2,1)), ...
        'C', mat2str(eye(2)), ...
        'D', mat2str(zeros(2,1)), ...
        'X0', mat2str(x0), ...
        'Position', [150 100 330 170]);

    add_block('simulink/Sinks/Scope', [mdl '/Scope'], ...
        'Position', [430 105 500 165]);

    add_line(mdl, 'ClosedLoopStateSpace/1', 'Scope/1');
    set_param(mdl, 'StopTime', num2str(tf));
    save_system(mdl);
    fprintf('Created optional Simulink model: %s.slx\n', mdl);
end
      

The MATLAB file uses lqr and also creates a simple optional Simulink model for the closed-loop system \( \dot{\mathbf{x} }=(A-BK)\mathbf{x} \) when Simulink is available.

Chapter28_Lesson5.nb

Notebook[{
Cell["Chapter28_Lesson5.nb", "Title"],
Cell["Continuous-time LQR preview for a double integrator.", "Text"],
Cell[BoxData[
 RowBox[{
  RowBox[{"A", "=", RowBox[{"{ {", RowBox[{
      RowBox[{"0", ",", "1"}], ",", RowBox[{"0", ",", "0"}]}], "} }"}]}], ";",
  RowBox[{"B", "=", RowBox[{"{ {", "0", "}", ",", "{", "1", "} }"}]}], ";",
  RowBox[{"Q", "=", RowBox[{"DiagonalMatrix", "[", RowBox[{"{", RowBox[{"10", ",", "1"}], "}"}], "]"}]}], ";",
  RowBox[{"R", "=", RowBox[{"{ {", "0.2", "} }"}]}], ";"}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"p", "=", RowBox[{"{ {", RowBox[{
      RowBox[{"p11", ",", "p12"}], ",", RowBox[{"p12", ",", "p22"}]}], "} }"}]}], ";",
  RowBox[{"care", "=", RowBox[{
      RowBox[{"Transpose", "[", "A", "]"}], ".", "p", "+", RowBox[{"p", ".", "A"}], "-",
      RowBox[{"p", ".", "B", ".", RowBox[{"Inverse", "[", "R", "]"}], ".", RowBox[{"Transpose", "[", "B", "]"}], ".", "p"}], "+", "Q"}]}], ";",
  RowBox[{"solutions", "=", RowBox[{"NSolve", "[", RowBox[{
      RowBox[{"Thread", "[", RowBox[{
          RowBox[{"Flatten", "[", "care", "]"}], "==", RowBox[{"ConstantArray", "[", RowBox[{"0", ",", "4"}], "]"}]}], "]"}],
      ",", RowBox[{"{", RowBox[{"p11", ",", "p12", ",", "p22"}], "}"}], ",", "Reals"}], "]"}]}], ";"}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"stableSolution", "=", RowBox[{"SelectFirst", "[", RowBox[{"solutions", ",",
      RowBox[{RowBox[{"Max", "[", RowBox[{"Re", "[",
          RowBox[{"Eigenvalues", "[", RowBox[{"A", "-", RowBox[{"B", ".", RowBox[{"Inverse", "[", "R", "]"}], ".", RowBox[{"Transpose", "[", "B", "]"}], ".", RowBox[{"(", RowBox[{"p", "/.", "#"}], ")"}]}]}], "]"}], "]"}], "]"}], "<", "0"}], "&"}]}], "]"}]}], ";",
  RowBox[{"P", "=", RowBox[{"p", "/.", "stableSolution"}]}], ";",
  RowBox[{"K", "=", RowBox[{RowBox[{"Inverse", "[", "R", "]"}], ".", RowBox[{"Transpose", "[", "B", "]"}], ".", "P"}]}], ";",
  RowBox[{"Acl", "=", RowBox[{"A", "-", RowBox[{"B", ".", "K"}]}]}], ";"}]], "Input"],
Cell[BoxData[
 RowBox[{"{", RowBox[{
  RowBox[{"MatrixForm", "[", "P", "]"}], ",",
  RowBox[{"MatrixForm", "[", "K", "]"}], ",",
  RowBox[{"Eigenvalues", "[", "Acl", "]"}]}], "}"}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"x0", "=", RowBox[{"{", RowBox[{"1", ",", "0"}], "}"}]}], ";",
  RowBox[{"x", "=", RowBox[{"NDSolveValue", "[", RowBox[{
      RowBox[{"{", RowBox[{
          RowBox[{RowBox[{"state", "'"}], "[", "t", "]"}], "==", RowBox[{"Acl", ".", RowBox[{"state", "[", "t", "]"}]}]}], ",",
          RowBox[{RowBox[{"state", "[", "0", "]"}], "==", "x0"}]}], "}"}], ",", "state", ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "8"}], "}"}]}], "]"}]}], ";",
  RowBox[{"Plot", "[", RowBox[{
      RowBox[{"Evaluate", "[", RowBox[{"x", "[", "t", "]"}], "]"}], ",",
      RowBox[{"{", RowBox[{"t", ",", "0", ",", "8"}], "}"}], ",",
      RowBox[{"PlotLegends", "->", RowBox[{"{", RowBox[{"\"position\"", ",", "\"velocity\""}], "}"}]}]}], "]"}]}]], "Input"]
}]
      

The Mathematica notebook expression solves the CARE symbolically or numerically as algebraic equations, selects the stabilizing solution, and simulates the closed-loop trajectory.

Problem 1 (Completing the Square): For the finite-horizon LQR problem, show that if \( P(t) \) satisfies the differential Riccati equation with terminal condition \( P(T)=S \), then the optimal input is \( \mathbf{u}^*=-R^{-1}B^TP\mathbf{x} \).

Solution: Differentiate \( \mathbf{x}^TP(t)\mathbf{x} \), substitute \( \dot{\mathbf{x} }=A\mathbf{x}+B\mathbf{u} \), and add the running cost. The cross term \( 2\mathbf{x}^TPB\mathbf{u} \) combines with \( \mathbf{u}^TR\mathbf{u} \) as

\[ \mathbf{u}^TR\mathbf{u}+2\mathbf{x}^TPB\mathbf{u}= \left(\mathbf{u}+R^{-1}B^TP\mathbf{x}\right)^TR \left(\mathbf{u}+R^{-1}B^TP\mathbf{x}\right) -\mathbf{x}^TPBR^{-1}B^TP\mathbf{x}. \]

The Riccati equation cancels the remaining state-only terms. Since \( R\succ 0 \), the square is minimized by setting it to zero, giving the stated feedback law.

Problem 2 (Double Integrator CARE): Let \( A=\begin{bmatrix}0&1\\0&0\end{bmatrix} \), \( B=\begin{bmatrix}0\\1\end{bmatrix} \), \( Q=\operatorname{diag}(q_1,q_2) \), and \( R=r>0 \). Derive the positive stabilizing solution \( P \).

Solution: Write \( P=\begin{bmatrix}p_{11}&p_{12}\\p_{12}&p_{22}\end{bmatrix} \). The CARE gives

\[ -\frac{p_{12}^2}{r}+q_1=0,\quad p_{11}-\frac{p_{12}p_{22} }{r}=0,\quad 2p_{12}-\frac{p_{22}^2}{r}+q_2=0. \]

Choose the positive branch for the stabilizing solution:

\[ p_{12}=\sqrt{q_1r},\quad p_{22}=\sqrt{r(2\sqrt{q_1r}+q_2)},\quad p_{11}=\frac{p_{12}p_{22} }{r}. \]

Hence \( K=R^{-1}B^TP=\begin{bmatrix}p_{12}/r&p_{22}/r\end{bmatrix} \).

Problem 3 (Lyapunov Stability from the CARE): Suppose \( P=P^T\succeq 0 \) solves the CARE and \( K=R^{-1}B^TP \). Show that the closed-loop derivative of \( V=\mathbf{x}^TP\mathbf{x} \) is negative semidefinite.

Solution: With \( A_{cl}=A-BK \),

\[ \dot{V}=\mathbf{x}^T(A_{cl}^TP+PA_{cl})\mathbf{x} =-\mathbf{x}^T(Q+K^TRK)\mathbf{x}\leq 0. \]

If the detectability condition holds, trajectories that produce zero penalty cannot hide unstable unobservable modes. Therefore the stabilizing Riccati solution makes \( A_{cl} \) Hurwitz.

Problem 4 (Cost-to-Go Identity): For the infinite-horizon LQR, prove that the optimal cost from \( \mathbf{x}_0 \) equals \( \mathbf{x}_0^TP\mathbf{x}_0 \).

Solution: Under the optimal feedback, \( \dot{V}=-(\mathbf{x}^TQ\mathbf{x}+ \mathbf{u}^{*T}R\mathbf{u}^*) \). Integrating from \( 0 \) to \( T \) gives

\[ \int_0^T(\mathbf{x}^TQ\mathbf{x}+ \mathbf{u}^{*T}R\mathbf{u}^*)dt =V(\mathbf{x}_0)-V(\mathbf{x}(T)). \]

Since the stabilizing closed loop gives \( \mathbf{x}(T)→ 0 \), the limit yields \( J_\infty^*(\mathbf{x}_0)=V(\mathbf{x}_0)= \mathbf{x}_0^TP\mathbf{x}_0 \).

Problem 5 (Effect of Input Weight): Explain why increasing \( R \) usually reduces gain magnitude.

Solution: The feedback gain is \( K=R^{-1}B^TP \). Although \( P \) also changes with \( R \), the direct effect of a larger positive definite \( R \) is to make input energy more expensive. The optimizer therefore accepts slower state regulation if that reduces \( \mathbf{u}^TR\mathbf{u} \). This is the formal version of the speed-versus-effort trade-off.

Optimal state feedback converts quadratic performance measures into a constructive feedback law. The finite-horizon problem produces a time-varying Riccati differential equation; the infinite-horizon LTI problem produces the algebraic Riccati equation. The stabilizing solution \( P \) is both a value-function matrix and a Lyapunov matrix, while \( K=R^{-1}B^TP \) gives the optimal state-feedback gain. This lesson is therefore the bridge between modern state-space design and the next course in optimal control.