Riccati Equation, Optimal Feedback Law, and Numerical Implementation

Matrix Riccati Differential Equation¶

From the quadratic-form identity derived in Part I,

Solving for $\dot{\boldsymbol{S}}$ gives the differential Riccati equation:

The terminal condition is

This equation is integrated backward in time.

Equivalent Backward-Time Form¶

Define the backward-time variable

Then

Therefore,

with

This form can be convenient for numerical integration because the independent variable increases forward.

Optimal Feedback Law¶

From Part I,

Using

we obtain

Define the time-varying feedback gain

Then

Closed-Loop Dynamics¶

Substituting the feedback law into the state equation gives

Therefore,

Once $\boldsymbol{S}(t)$ has been computed backward in time, the closed-loop state is propagated forward from $\boldsymbol{x}(t_0)=\boldsymbol{x}_0$.

Backward Planning and Forward Execution¶

The LQR solution consists of two distinct integrations:

1. Integrate the Riccati equation backward:

2. Compute $\boldsymbol{K}(t)$ from $\boldsymbol{S}(t)$.

3. Integrate the closed-loop state equation forward:

4. Recover the optimal control:

This is the exact mathematical realization of Bellman’s idea: plan backward and execute forward.

Interpretation of the Matrix $\boldsymbol{S}(t)$¶

The value function is

Thus, $\boldsymbol{S}(t)$ describes the local curvature of the optimal cost-to-go with respect to the state.

Large eigenvalues of $\boldsymbol{S}(t)$ indicate directions in state space for which deviations are expensive. Small eigenvalues indicate directions that are relatively inexpensive.

Relationship to the Costate¶

Along the optimal trajectory,

Hence,

Substituting this into the Pontryagin control law

recovers the HJB feedback law.

Consistency with Pontryagin’s Minimum Principle¶

For the LQR problem, the PMP Hamiltonian is

Stationarity gives

The costate equation is

Assume

Then

Substituting the closed-loop dynamics and equating terms yields the same differential Riccati equation.

Therefore, HJB and PMP are fully consistent for the LQR problem.

Why the LQR Problem Is Special¶

The LQR problem is unusually tractable because:

• the dynamics are linear;

• the cost is quadratic;

• the optimal value function remains quadratic;

• the HJB PDE reduces to a matrix ODE;

• the optimal control becomes a linear feedback law.

A small nonlinear perturbation may destroy this structure and return the problem to a high-dimensional nonlinear PDE.

Curse of Dimensionality¶

For a general nonlinear system with $n_x$ states, the HJB equation is defined over an $(n_x+1)$-dimensional state-time domain.

If each state dimension is discretized with $M$ grid points, the number of state-grid points is approximately

This exponential growth is the curse of dimensionality.

The LQR problem avoids this difficulty because the unknown value function is represented by the matrix $\boldsymbol{S}(t)$ rather than by values on a full state-space grid.

Time-Varying Versus Time-Invariant LQR¶

For time-varying matrices, the Riccati equation is

For an infinite-horizon time-invariant problem, one seeks a constant solution $\boldsymbol{S}_\infty$ satisfying

This is the continuous algebraic Riccati equation.

Finite-Horizon MATLAB Implementation¶

The following template integrates the Riccati equation backward and the state equation forward.

function demo_finite_horizon_lqr

    t0 = 0;
    tf = 10;

    A = [0 1; -2 -0.5];
    B = [0; 1];
    Q = eye(2);
    R = 1;
    Sf = eye(2);

    n = size(A,1);

    % Backward Riccati equation
    riccati = @(t,svec) riccati_rhs( ...
        t, svec, A, B, Q, R, n);

    [tS, Svec] = ode45( ...
        riccati, [tf t0], Sf(:));

    % Reverse for increasing-time interpolation
    tS = flipud(tS);
    Svec = flipud(Svec);

    Sfun = @(t) reshape( ...
        interp1(tS, Svec, t, 'pchip'), n, n);

    x0 = [1; 0];

    % Closed-loop state dynamics
    closed_loop = @(t,x) ...
        (A - B*(R\(B.'*Sfun(t))))*x;

    [tx, x] = ode45(closed_loop, [t0 tf], x0);

    u = zeros(length(tx),1);
    for k = 1:length(tx)
        S = Sfun(tx(k));
        K = R\(B.'*S);
        u(k) = -K*x(k,:).';
    end

    figure;
    plot(tx,x,'LineWidth',1.5);
    xlabel('Time');
    ylabel('State');
    legend('x_1','x_2');
    grid on;

    figure;
    plot(tx,u,'LineWidth',1.5);
    xlabel('Time');
    ylabel('Control');
    grid on;
end

function ds = riccati_rhs(~,svec,A,B,Q,R,n)

    S = reshape(svec,n,n);

    dS = -A.'*S - S*A ...
         + S*B*(R\B.')*S - Q;

    ds = dS(:);
end

Numerical Symmetry Preservation¶

In exact arithmetic, $\boldsymbol{S}(t)$ remains symmetric if $\boldsymbol{S}_f$ is symmetric.

Numerical integration may introduce small asymmetry. A common correction is

This should be used only to remove numerical roundoff, not to hide a significant implementation error.

Verification of the Riccati Solution¶

A computed solution should be checked using the residual

A valid numerical solution should satisfy

Also verify

Verification of the Value Function¶

For the optimal state and control,

The same quantity should be recovered from numerical integration:

The difference

should be small.

Feedback Interpretation¶

The HJB derivation generates

Thus, the control depends on the current state, not only on the nominal state trajectory.

If a disturbance changes the state, the same gain matrix maps the new state to a corrective control.

Open-Loop and Closed-Loop Views¶

The LQR solution can be interpreted in two ways:

1. An open-loop optimal trajectory is obtained by propagating the optimal closed-loop system from a specified initial state.

2. The feedback law itself is valid for every state in the domain.

The HJB solution therefore contains more information than a single open-loop trajectory.

Common Derivation Errors¶

Dropping the Factor of One-Half¶

The factors of one-half simplify derivatives. Removing them inconsistently produces factor-of-two errors.

Incorrect Gradient Orientation¶

Use a consistent column-gradient convention:

Forgetting Matrix Symmetry¶

The derivative of

is $\boldsymbol{S}\boldsymbol{x}$ only when $\boldsymbol{S}$ is symmetric.

Incorrect Riccati Sign Convention¶

The finite-horizon Riccati equation is commonly written either in forward-time derivative form with a terminal condition or in backward-time variable form. Mixing the two conventions causes sign errors.

Integrating from the Wrong Boundary¶

The Riccati equation must satisfy

Assuming Numerical Optimality Without Verification¶

The Riccati residual, terminal condition, closed-loop dynamics, and cost identity must all be checked.

Recommended Computational Workflow¶

1. Define $\boldsymbol{A}(t)$, $\boldsymbol{B}(t)$, $\boldsymbol{Q}(t)$, $\boldsymbol{R}(t)$, and $\boldsymbol{S}_f$.

2. Confirm symmetry and definiteness assumptions.

3. Integrate the Riccati equation backward.

4. Verify the terminal condition.

5. Verify symmetry of $\boldsymbol{S}(t)$.

6. Compute $\boldsymbol{K}(t)$.

7. Integrate the closed-loop state equation forward.

8. Compute $\boldsymbol{u}^*(t)$.

9. Evaluate the cost numerically.

10. Compare it with the value-function prediction.

11. Check Riccati and state-equation residuals.

Summary¶

1. The HJB equation reduces to the differential Riccati equation for the LQR problem.

2. The matrix $\boldsymbol{S}(t)$ is integrated backward from $\boldsymbol{S}_f$.

3. The optimal feedback gain is

4. The optimal control is

5. The costate is

6. HJB and PMP produce the same Riccati equation.

7. The LQR problem avoids the curse of dimensionality because the value function is represented by a matrix.

8. Numerical solutions require residual, symmetry, boundary-condition, and cost checks.

Connection. A feedback law is one component of a larger system, so the discussion now places optimization alongside guidance, navigation, estimation, and inner-loop control.