Chapter 28: Advanced Directions and Links to Other Courses

Lesson 1: Adaptive Critics and Approximate Dynamic Programming (Pointer to Neural/Intelligent Control)

This lesson connects adaptive control to approximate dynamic programming (ADP), reinforcement learning, and neural control. The central idea is to adapt an approximation of the long-horizon value function, called the critic, and to use that approximation to improve a control policy, called the actor. The treatment remains control-theoretic: Bellman equations, Hamilton-Jacobi-Bellman equations, policy iteration, temporal- difference residuals, approximation error, admissibility, and stability limitations are emphasized before neural-network implementation details.

1. Why Adaptive Critics Belong Near the End of an Adaptive-Control Course

Earlier chapters adapted plant parameters, controller gains, or both in order to obtain tracking and regulation despite uncertainty. Adaptive critics adapt a different mathematical object: an approximation to the cost-to-go. For a state \(x_k\) and policy \(\mu\), the cost-to-go measures the cumulative future penalty generated when the policy is followed from that state.

\[ V^{\mu}(x_k) = \mathbb{E}_{\mu}\!\left[ \sum_{i=0}^{\infty}\gamma^i \ell(x_{k+i},u_{k+i}) \;\middle|\;x_k \right], \qquad 0 < \gamma < 1. \]

Here \(\ell(x,u)\ge 0\) is a stage cost and \(\gamma\) is a discount factor. The critic learns or approximates \(V^{\mu}\); the actor changes \(\mu\) to reduce this value. Thus, the adaptation target is not merely a plant coefficient or tracking-error gain, but a solution of a dynamic-programming equation.

This distinction is important. A model-reference adaptive controller is usually designed around a chosen error model and a Lyapunov proof. An ADP controller is designed around an optimization objective and a Bellman or HJB equation. Its weights may be updated online, but that fact alone does not provide Lyapunov stability, robustness margins, safe exploration, or parameter convergence.

Question Classical adaptive control Adaptive critic / ADP
Main learned object Plant parameters or controller parameters Value, action-value, value gradient, and/or policy
Primary equation Error dynamics and Lyapunov derivative Bellman equation or HJB equation
Primary goal Tracking/regulation with uncertainty Minimize cumulative closed-loop cost
Typical excitation issue Persistent excitation for parameter convergence State-action coverage for value/policy learning
Typical proof burden Boundedness, tracking convergence Bellman approximation, policy improvement, stability, safety

2. Discrete-Time Optimal-Control Foundation

Consider a deterministic nonlinear system

\[ x_{k+1}=f(x_k,u_k), \qquad x_k\in\mathbb{R}^{n},\quad u_k\in\mathcal{U}. \]

For a stationary policy \(u_k=\mu(x_k)\), the value function satisfies the policy Bellman equation

\[ V^{\mu}(x) = \ell(x,\mu(x)) + \gamma V^{\mu}\!\left(f(x,\mu(x))\right). \]

The optimal value is the pointwise minimum over all admissible policies:

\[ V^{\star}(x)=\inf_{\mu}V^{\mu}(x). \]

Bellman's principle of optimality yields the Bellman optimality equation

\[ V^{\star}(x) = \min_{u\in\mathcal{U}} \left\{ \ell(x,u)+\gamma V^{\star}(f(x,u)) \right\}. \]

Define the Bellman optimality operator \(T\) by

\[ (TV)(x) = \min_{u\in\mathcal{U}} \left\{\ell(x,u)+\gamma V(f(x,u))\right\}. \]

Under bounded costs and \(0<\gamma<1\), this operator is a contraction in the sup norm. For arbitrary bounded functions \(V\) and \(W\), choose an action that minimizes the expression for one function and use the elementary inequality between minima:

\[ \begin{aligned} |(TV)(x)-(TW)(x)| &\le \gamma\max_{u\in\mathcal{U}} |V(f(x,u))-W(f(x,u))| \\ &\le \gamma\|V-W\|_{\infty}. \end{aligned} \]

Taking the supremum over \(x\) gives

\[ \|TV-TW\|_{\infty} \le \gamma\|V-W\|_{\infty}. \]

Therefore, the Banach fixed-point theorem gives a unique fixed point \(V^{\star}=TV^{\star}\), and exact value iteration \(V_{j+1}=TV_j\) converges geometrically. ADP begins when exact storage or exact application of \(T\) is impossible and a parameterized approximation replaces the unrestricted function.

3. Critic, Actor, and Optional Model

The critic estimates a value-related quantity, while the actor generates the control input. Some architectures also contain an identifier or learned model. The most common mathematical choices are:

  • State-value critic: \(\widehat V(x;w_c)\) approximates \(V^{\mu}(x)\).
  • Action-value critic: \(\widehat Q(x,u;w_c)\) approximates the cost of taking \(u\) and then following the policy.
  • Derivative critic: approximates \(\nabla_x V(x)\), which is closely related to dual-heuristic programming.
  • Actor: \(\widehat\mu(x;w_a)\) approximates the minimizing control law.
flowchart TD
  R["Reference, objective, and constraints"] --> A["Actor: control u = policy(x)"]
  A --> P["Plant or environment"]
  P --> X["Measured transition: x, u, cost, x_next"]
  X --> C["Critic: estimate long-horizon cost"]
  C --> TD["Bellman or TD residual"]
  TD --> C
  C --> PI["Policy improvement signal"]
  PI --> A
  X --> ID["Optional identifier or learned model"]
  ID --> C
  ID --> PI
        

Historically, adaptive-critic terminology includes heuristic dynamic programming (value critic), dual heuristic programming (value-gradient critic), and globalized dual heuristic programming (both value and gradient information). Modern reinforcement-learning terminology more often uses critic, actor-critic, Q-learning, and policy gradient. The underlying control equations remain Bellman or HJB equations.

4. Function Approximation and the Projected Bellman Equation

Let a linear-in-the-weights critic be

\[ \widehat V(x;w_c)=w_c^{\mathsf T}\phi(x), \qquad \phi(x)\in\mathbb{R}^{N_c}. \]

The features may be monomials, radial basis functions, Fourier features, or the penultimate-layer outputs of a neural network. A neural critic can be written more generally as

\[ \widehat V(x;w_c)=\mathcal{N}_c(x;w_c). \]

For a fixed policy, the Bellman operator is

\[ (T_{\mu}V)(x) =\ell(x,\mu(x)) +\gamma V(f(x,\mu(x))). \]

Because the exact Bellman image need not lie in the approximation space, the critic commonly solves a projected Bellman equation

\[ \Phi w_c = \Pi T_{\mu}(\Phi w_c), \]

where \(\Pi\) is a projection induced by a sampling distribution and inner product. For sampled transitions \((x_k,u_k,c_k,x_{k+1})\), define the temporal-difference residual

\[ \delta_k =c_k+\gamma\widehat V(x_{k+1};w_c) -\widehat V(x_k;w_c). \]

Classical semi-gradient TD(0) treats the bootstrap target as locally constant:

\[ w_{c,k+1} =w_{c,k}+\alpha_c\delta_k\phi(x_k). \]

A residual-gradient update differentiates both terms:

\[ w_{c,k+1} =w_{c,k} -\alpha_c\delta_k \left[ \gamma\nabla_{w_c}\widehat V(x_{k+1};w_c) -\nabla_{w_c}\widehat V(x_k;w_c) \right]. \]

The implementations in this lesson use a normalized residual-gradient step because it is transparent for a one-parameter quadratic critic. In larger problems, TD, least-squares TD, gradient-TD, replay buffers, target networks, and regularized nonlinear optimization address different bias, variance, conditioning, and off-policy stability issues.

For linear features under a stationary policy, the expected TD normal equation has the form

\[ A w_c=b, \]

\[ A=\mathbb{E}\!\left[ \phi(x_k) \left(\phi(x_k)-\gamma\phi(x_{k+1})\right)^{\mathsf T} \right], \qquad b=\mathbb{E}\!\left[\phi(x_k)c_k\right]. \]

This equation reveals the analogue of excitation in adaptive control: the data must make \(A\) sufficiently informative and well-conditioned. Poor state visitation or redundant features can make critic identification ill-conditioned even when the closed loop remains bounded.

5. Policy Iteration and the Policy-Improvement Argument

Exact policy iteration alternates between:

  1. Policy evaluation: solve \(V^{\mu_j}=T_{\mu_j}V^{\mu_j}\).
  2. Policy improvement: choose \(\mu_{j+1}\) greedily with respect to \(V^{\mu_j}\).

\[ \mu_{j+1}(x) \in\arg\min_{u\in\mathcal{U}} \left\{ \ell(x,u)+\gamma V^{\mu_j}(f(x,u)) \right\}. \]

By the greedy definition,

\[ T_{\mu_{j+1}}V^{\mu_j} \le T_{\mu_j}V^{\mu_j} =V^{\mu_j}. \]

Repeated application of the monotone operator \(T_{\mu_{j+1}}\) then gives

\[ V^{\mu_{j+1}} =\lim_{m\to\infty}T_{\mu_{j+1}}^{m}V^{\mu_j} \le V^{\mu_j}. \]

Thus exact policy improvement cannot increase the cost. Approximate policy iteration replaces exact evaluation and exact minimization by approximations. If the critic error is large, the actor optimization is inaccurate, or the policy leaves the region represented by the data, monotonic improvement is no longer automatic.

An action-value critic can remove explicit model dependence. Define

\[ Q^{\mu}(x,u) =\ell(x,u)+\gamma V^{\mu}(f(x,u)). \]

Then the greedy actor is simply \(\mu^{+}(x)\in\arg\min_u Q^{\mu}(x,u)\). This is the route from model-assisted ADP to model-free Q-learning and actor-critic methods.

6. Continuous-Time Link: HJB and Value Gradients

For the input-affine continuous-time system

\[ \dot x=f(x)+g(x)u \]

and infinite-horizon cost

\[ J(x_0,u)=\int_{0}^{\infty} \left[x^{\mathsf T}Qx+u^{\mathsf T}Ru\right]dt, \qquad Q\succeq 0,\quad R\succ 0, \]

the stationary HJB equation is

\[ 0=\min_u\left\{ x^{\mathsf T}Qx+u^{\mathsf T}Ru +\nabla V^{\star}(x)^{\mathsf T}[f(x)+g(x)u] \right\}. \]

Differentiating the Hamiltonian with respect to \(u\) gives

\[ 2Ru+g(x)^{\mathsf T}\nabla V^{\star}(x)=0, \]

\[ u^{\star}(x) =-\frac{1}{2}R^{-1}g(x)^{\mathsf T}\nabla V^{\star}(x). \]

This formula explains why value-gradient critics are attractive in nonlinear control: the optimal feedback is directly related to \(\nabla V^{\star}\). It also exposes a difficulty: an inaccurate critic gradient can produce an unsafe actor even when the value approximation itself appears numerically reasonable. Continuous- time integral reinforcement learning and actor-critic methods avoid direct differentiation of unknown drift dynamics by integrating Bellman identities over finite intervals; a full treatment belongs in a neural or learning-based control course.

7. Worked Scalar LQR: Exact Solution and Quadratic Critic

The multilingual laboratory uses the scalar plant

\[ x_{k+1}=a x_k+b u_k \]

with discounted quadratic cost

\[ J=\sum_{k=0}^{\infty}\gamma^k \left(qx_k^2+ru_k^2\right), \qquad q>0,\quad r>0. \]

Assume the value function is quadratic: \(V(x)=p x^2\). For a current coefficient \(p\), the one-step greedy problem is

\[ \min_u\left\{ qx^2+ru^2+\gamma p(ax+bu)^2 \right\}. \]

The first-order condition is

\[ 2ru+2\gamma p b(ax+bu)=0. \]

Therefore the greedy actor has the form \(u=-k(p)x\), where

\[ k(p)=\frac{\gamma bpa}{r+\gamma b^2p}. \]

Substitution into the Bellman equation gives the scalar discounted Riccati mapping

\[ \mathcal{F}(p) =q+\gamma a^2p -\frac{(\gamma abp)^2}{r+\gamma b^2p}. \]

The optimal coefficient satisfies \(p^{\star}=\mathcal{F}(p^{\star})\), and the optimal gain is \(k^{\star}=k(p^{\star})\).

The critic used in the code is \(\widehat V(x;p_c)=p_cx^2\). For the observed transition,

\[ \delta_k =qx_k^2+ru_k^2 +\gamma p_c x_{k+1}^2-p_cx_k^2. \]

Its normalized residual-gradient update is

\[ p_{c,k+1} =\operatorname{Proj}_{[p_{\min},p_{\max}]} \left[ p_{c,k} -\alpha_c\delta_k \frac{\gamma x_{k+1}^2-x_k^2} {1+(\gamma x_{k+1}^2-x_k^2)^2} \right]. \]

The actor is moved gradually toward the greedy gain implied by the critic:

\[ k_{a,k+1} =\operatorname{Proj}_{[k_{\min},k_{\max}]} \left[ k_{a,k}+\alpha_a(k(p_{c,k+1})-k_{a,k}) \right]. \]

This is a model-assisted actor update because the formula uses \(a\) and \(b\). The critic itself is updated from sampled transitions. To remove model dependence, one can learn a quadratic action-value function

\[ \widehat Q(x,u) =\begin{bmatrix}x\\u\end{bmatrix}^{\mathsf T} \begin{bmatrix}h_{xx}&h_{xu}\\h_{xu}&h_{uu}\end{bmatrix} \begin{bmatrix}x\\u\end{bmatrix}, \]

for which the greedy input is \(u=-h_{uu}^{-1}h_{ux}x\), provided \(h_{uu}>0\).

8. Admissibility, Exploration, and Stability Limitations

ADP proofs commonly begin from an admissible policy: a policy that stabilizes the plant in the region of interest and produces finite cost. For the scalar linear policy \(u=-kx\), the closed-loop dynamics are

\[ x_{k+1}=(a-bk)x_k. \]

A sufficient and necessary condition for asymptotic stability of this scalar discrete-time closed loop is

\[ |a-bk|<1. \]

The code starts from such a gain. Exploration is added as a decaying perturbation to the input so that the critic observes informative transitions. This creates a fundamental conflict:

  • Too little exploration can make critic parameters unidentifiable.
  • Too much exploration can violate constraints or destabilize the plant.
  • Rapid actor changes can invalidate the critic's policy-evaluation target.
  • Function-approximation error can destroy exact policy improvement.

Projection, normalization, gain-rate separation, state constraints, backup controllers, and Lyapunov or barrier certificates are therefore not optional details in safety-critical control. The simple script uses projection and a slower actor tracking step, but it is not a general safety certificate.

\[ \alpha_a \ll \alpha_c \]

is a common two-time-scale design intuition: let the critic approximately evaluate the current policy before the actor moves substantially. The code uses equal nominal rates only because the scalar critic converges quickly; larger systems generally require more careful time-scale design.

9. Algorithm Used in the Multilingual Laboratory

flowchart TD
  S["Choose plant, cost, discount, and admissible initial gain"] --> E["Reset state for a training episode"]
  E --> U["Apply actor control plus bounded exploration"]
  U --> T["Measure next state and stage cost"]
  T --> D["Compute Bellman or TD residual"]
  D --> C["Update and project critic parameter"]
  C --> G["Compute greedy gain from current critic"]
  G --> A["Move and project actor gain"]
  A --> Q["More samples in this episode?"]
  Q -->|yes| U
  Q -->|no| N["More episodes?"]
  N -->|yes| E
  N -->|no| V["Compare learned policy with Riccati reference"]
        

The numerical parameters are \(a=1.05\), \(b=0.5\), \(q=1\), \(r=0.1\), and \(\gamma=0.98\). The open-loop plant is unstable because \(|a|>1\). The initial gain \(k_0=0.8\) gives the stable pole \(a-bk_0=0.65\). With the fixed random seed, the learned values are expected to be close to \(p^{\star}\approx1.337912\) and \(k^{\star}\approx1.609103\).

10. Python Implementation — Chapter28_Lesson1.py

This implementation uses NumPy only for reproducible random sampling and clipping. For larger ADP systems, relevant Python libraries include NumPy/SciPy for numerical linear algebra, PyTorch or JAX for differentiable neural critics, and Gymnasium-style environments for standardized simulation interfaces.

"""Chapter 28, Lesson 1: adaptive critic for a scalar discounted LQR.

The critic estimates V(x) = p*x^2 from temporal-difference data. The actor
uses the critic in the scalar greedy policy-improvement formula u = -k*x.
This is a model-assisted pedagogical ADP example: the critic is data-driven,
while the actor update uses the known input coefficient b.
"""

from __future__ import annotations

import math
from dataclasses import dataclass

import numpy as np


@dataclass(frozen=True)
class Problem:
    a: float = 1.05
    b: float = 0.50
    q: float = 1.00
    r: float = 0.10
    gamma: float = 0.98


def discounted_riccati_reference(problem: Problem) -> tuple[float, float]:
    """Return the scalar discounted-Riccati value coefficient and feedback gain."""
    p = problem.q
    for _ in range(100_000):
        denominator = problem.r + problem.gamma * problem.b**2 * p
        p_next = (
            problem.q
            + problem.gamma * problem.a**2 * p
            - (problem.gamma * problem.a * problem.b * p) ** 2 / denominator
        )
        if abs(p_next - p) < 1.0e-14:
            p = p_next
            break
        p = p_next

    k = (
        problem.gamma
        * problem.b
        * p
        * problem.a
        / (problem.r + problem.gamma * problem.b**2 * p)
    )
    return p, k


def train_adaptive_critic(
    problem: Problem,
    episodes: int = 500,
    steps_per_episode: int = 80,
    seed: int = 7,
) -> tuple[float, float, list[tuple[int, float, float]]]:
    """Train a quadratic critic and scalar actor with normalized TD updates."""
    rng = np.random.default_rng(seed)

    critic_p = 0.50
    actor_k = 0.80  # admissible initial controller: |a - b*k| < 1
    critic_rate = 0.08
    actor_rate = 0.08
    history: list[tuple[int, float, float]] = []

    for episode in range(episodes):
        x = float(rng.uniform(-2.0, 2.0))
        exploration_scale = 0.25 * math.exp(-episode / 180.0)

        for _ in range(steps_per_episode):
            exploration = exploration_scale * float(rng.normal())
            u = -actor_k * x + exploration
            x_next = problem.a * x + problem.b * u
            stage_cost = problem.q * x * x + problem.r * u * u

            # Bellman/TD residual for V_hat(x) = critic_p*x^2.
            delta = (
                stage_cost
                + problem.gamma * critic_p * x_next * x_next
                - critic_p * x * x
            )
            critic_gradient = problem.gamma * x_next * x_next - x * x

            # Normalization bounds a large one-sample update.
            critic_p -= (
                critic_rate
                * delta
                * critic_gradient
                / (1.0 + critic_gradient * critic_gradient)
            )
            critic_p = float(np.clip(critic_p, 1.0e-6, 50.0))

            # Greedy policy improvement induced by the current quadratic critic.
            target_k = (
                problem.gamma
                * problem.b
                * critic_p
                * problem.a
                / (problem.r + problem.gamma * problem.b**2 * critic_p)
            )
            actor_k += actor_rate * (target_k - actor_k)
            actor_k = float(np.clip(actor_k, 0.0, 3.0))
            x = x_next

        if episode % 50 == 0 or episode == episodes - 1:
            history.append((episode, critic_p, actor_k))

    return critic_p, actor_k, history


def rollout_cost(
    problem: Problem,
    gain: float,
    x0: float = 2.0,
    horizon: int = 60,
) -> tuple[float, list[float]]:
    """Evaluate a deterministic discounted closed-loop rollout."""
    x = x0
    total = 0.0
    trajectory = [x]
    discount = 1.0

    for _ in range(horizon):
        u = -gain * x
        total += discount * (problem.q * x * x + problem.r * u * u)
        x = problem.a * x + problem.b * u
        trajectory.append(x)
        discount *= problem.gamma

    return total, trajectory


def main() -> None:
    problem = Problem()
    p_star, k_star = discounted_riccati_reference(problem)
    learned_p, learned_k, history = train_adaptive_critic(problem)

    initial_cost, _ = rollout_cost(problem, gain=0.80)
    learned_cost, learned_trajectory = rollout_cost(problem, gain=learned_k)

    print("episode, critic_p, actor_k")
    for episode, critic_p, actor_k in history:
        print(f"{episode:4d}, {critic_p:10.6f}, {actor_k:10.6f}")

    print("\nReference discounted-Riccati solution")
    print(f"p* = {p_star:.9f}")
    print(f"k* = {k_star:.9f}")
    print("\nLearned adaptive-critic solution")
    print(f"p_hat = {learned_p:.9f}")
    print(f"k_hat = {learned_k:.9f}")
    print(f"closed-loop pole = {problem.a - problem.b * learned_k:.9f}")
    print(f"initial-policy rollout cost = {initial_cost:.9f}")
    print(f"learned-policy rollout cost = {learned_cost:.9f}")
    print("first ten learned-policy states:")
    print(" ".join(f"{x:.6f}" for x in learned_trajectory[:10]))


if __name__ == "__main__":
    main()

11. C++ Implementation — Chapter28_Lesson1.cpp

The C++17 implementation uses only the standard library. In larger systems, Eigen provides matrix algebra, mlpack provides machine-learning components, and LibTorch provides the C++ interface to PyTorch neural networks.

// Chapter 28, Lesson 1: adaptive critic for a scalar discounted LQR.
// Compile: g++ -std=c++17 -O2 Chapter28_Lesson1.cpp -o Chapter28_Lesson1

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <random>
#include <stdexcept>
#include <utility>
#include <vector>

struct Problem {
    double a = 1.05;
    double b = 0.50;
    double q = 1.00;
    double r = 0.10;
    double gamma = 0.98;
};

struct Snapshot {
    int episode;
    double criticP;
    double actorK;
};

std::pair<double, double> discountedRiccatiReference(const Problem& problem) {
    double p = problem.q;
    for (int iteration = 0; iteration < 100000; ++iteration) {
        const double denominator =
            problem.r + problem.gamma * problem.b * problem.b * p;
        const double cross = problem.gamma * problem.a * problem.b * p;
        const double pNext = problem.q + problem.gamma * problem.a * problem.a * p
                           - cross * cross / denominator;
        if (std::abs(pNext - p) < 1.0e-14) {
            p = pNext;
            break;
        }
        p = pNext;
    }

    const double k = problem.gamma * problem.b * p * problem.a
                   / (problem.r + problem.gamma * problem.b * problem.b * p);
    return {p, k};
}

std::tuple<double, double, std::vector<Snapshot>> trainAdaptiveCritic(
    const Problem& problem,
    int episodes = 500,
    int stepsPerEpisode = 80,
    unsigned int seed = 7U) {

    if (episodes <= 0 || stepsPerEpisode <= 0) {
        throw std::invalid_argument("episodes and stepsPerEpisode must be positive");
    }

    std::mt19937 generator(seed);
    std::uniform_real_distribution<double> initialState(-2.0, 2.0);
    std::normal_distribution<double> standardNormal(0.0, 1.0);

    double criticP = 0.50;
    double actorK = 0.80;
    const double criticRate = 0.08;
    const double actorRate = 0.08;
    std::vector<Snapshot> history;

    for (int episode = 0; episode < episodes; ++episode) {
        double x = initialState(generator);
        const double explorationScale = 0.25 * std::exp(-episode / 180.0);

        for (int step = 0; step < stepsPerEpisode; ++step) {
            const double exploration = explorationScale * standardNormal(generator);
            const double u = -actorK * x + exploration;
            const double xNext = problem.a * x + problem.b * u;
            const double stageCost = problem.q * x * x + problem.r * u * u;

            const double delta = stageCost
                               + problem.gamma * criticP * xNext * xNext
                               - criticP * x * x;
            const double criticGradient = problem.gamma * xNext * xNext - x * x;

            criticP -= criticRate * delta * criticGradient
                     / (1.0 + criticGradient * criticGradient);
            criticP = std::clamp(criticP, 1.0e-6, 50.0);

            const double targetK = problem.gamma * problem.b * criticP * problem.a
                                 / (problem.r
                                    + problem.gamma * problem.b * problem.b * criticP);
            actorK += actorRate * (targetK - actorK);
            actorK = std::clamp(actorK, 0.0, 3.0);
            x = xNext;
        }

        if (episode % 50 == 0 || episode == episodes - 1) {
            history.push_back({episode, criticP, actorK});
        }
    }

    return {criticP, actorK, history};
}

std::pair<double, std::vector<double>> rolloutCost(
    const Problem& problem,
    double gain,
    double x0 = 2.0,
    int horizon = 60) {

    double x = x0;
    double total = 0.0;
    double discount = 1.0;
    std::vector<double> trajectory{x};

    for (int step = 0; step < horizon; ++step) {
        const double u = -gain * x;
        total += discount * (problem.q * x * x + problem.r * u * u);
        x = problem.a * x + problem.b * u;
        trajectory.push_back(x);
        discount *= problem.gamma;
    }

    return {total, trajectory};
}

int main() {
    try {
        const Problem problem;
        const auto [pStar, kStar] = discountedRiccatiReference(problem);
        const auto [learnedP, learnedK, history] = trainAdaptiveCritic(problem);
        const auto [initialCost, ignored] = rolloutCost(problem, 0.80);
        const auto [learnedCost, learnedTrajectory] = rolloutCost(problem, learnedK);

        std::cout << std::fixed << std::setprecision(9);
        std::cout << "episode, critic_p, actor_k\n";
        for (const Snapshot& snapshot : history) {
            std::cout << std::setw(4) << snapshot.episode << ", "
                      << snapshot.criticP << ", " << snapshot.actorK << '\n';
        }

        std::cout << "\nReference discounted-Riccati solution\n";
        std::cout << "p* = " << pStar << '\n';
        std::cout << "k* = " << kStar << '\n';
        std::cout << "\nLearned adaptive-critic solution\n";
        std::cout << "p_hat = " << learnedP << '\n';
        std::cout << "k_hat = " << learnedK << '\n';
        std::cout << "closed-loop pole = " << problem.a - problem.b * learnedK << '\n';
        std::cout << "initial-policy rollout cost = " << initialCost << '\n';
        std::cout << "learned-policy rollout cost = " << learnedCost << '\n';
        std::cout << "first ten learned-policy states:\n";
        for (std::size_t index = 0;
             index < std::min<std::size_t>(10, learnedTrajectory.size());
             ++index) {
            std::cout << learnedTrajectory[index]
                      << (index + 1 == 10 ? '\n' : ' ');
        }
        return 0;
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
}

12. Java Implementation — Chapter28_Lesson1.java

The Java implementation also uses only the standard library. EJML is a lightweight choice for dense and sparse matrix computations, while ND4J and Deeplearning4j support tensor operations and neural critics.

// Chapter 28, Lesson 1: adaptive critic for a scalar discounted LQR.
// Compile: javac Chapter28_Lesson1.java
// Run:     java Chapter28_Lesson1

import java.util.ArrayList;
import java.util.List;
import java.util.Locale;
import java.util.Random;

public final class Chapter28_Lesson1 {
    private Chapter28_Lesson1() {}

    private static final class Problem {
        final double a = 1.05;
        final double b = 0.50;
        final double q = 1.00;
        final double r = 0.10;
        final double gamma = 0.98;
    }

    private static final class Snapshot {
        final int episode;
        final double criticP;
        final double actorK;

        Snapshot(int episode, double criticP, double actorK) {
            this.episode = episode;
            this.criticP = criticP;
            this.actorK = actorK;
        }
    }

    private static final class TrainingResult {
        final double criticP;
        final double actorK;
        final List<Snapshot> history;

        TrainingResult(double criticP, double actorK, List<Snapshot> history) {
            this.criticP = criticP;
            this.actorK = actorK;
            this.history = history;
        }
    }

    private static final class RolloutResult {
        final double cost;
        final List<Double> trajectory;

        RolloutResult(double cost, List<Double> trajectory) {
            this.cost = cost;
            this.trajectory = trajectory;
        }
    }

    private static double clamp(double value, double lower, double upper) {
        return Math.max(lower, Math.min(upper, value));
    }

    private static double[] discountedRiccatiReference(Problem problem) {
        double p = problem.q;
        for (int iteration = 0; iteration < 100_000; iteration++) {
            double denominator = problem.r
                + problem.gamma * problem.b * problem.b * p;
            double cross = problem.gamma * problem.a * problem.b * p;
            double pNext = problem.q
                + problem.gamma * problem.a * problem.a * p
                - cross * cross / denominator;
            if (Math.abs(pNext - p) < 1.0e-14) {
                p = pNext;
                break;
            }
            p = pNext;
        }

        double k = problem.gamma * problem.b * p * problem.a
            / (problem.r + problem.gamma * problem.b * problem.b * p);
        return new double[] {p, k};
    }

    private static TrainingResult trainAdaptiveCritic(
        Problem problem,
        int episodes,
        int stepsPerEpisode,
        long seed) {

        if (episodes <= 0 || stepsPerEpisode <= 0) {
            throw new IllegalArgumentException(
                "episodes and stepsPerEpisode must be positive");
        }

        Random random = new Random(seed);
        double criticP = 0.50;
        double actorK = 0.80;
        double criticRate = 0.08;
        double actorRate = 0.08;
        List<Snapshot> history = new ArrayList<>();

        for (int episode = 0; episode < episodes; episode++) {
            double x = -2.0 + 4.0 * random.nextDouble();
            double explorationScale = 0.25 * Math.exp(-episode / 180.0);

            for (int step = 0; step < stepsPerEpisode; step++) {
                double exploration = explorationScale * random.nextGaussian();
                double u = -actorK * x + exploration;
                double xNext = problem.a * x + problem.b * u;
                double stageCost = problem.q * x * x + problem.r * u * u;

                double delta = stageCost
                    + problem.gamma * criticP * xNext * xNext
                    - criticP * x * x;
                double criticGradient = problem.gamma * xNext * xNext - x * x;

                criticP -= criticRate * delta * criticGradient
                    / (1.0 + criticGradient * criticGradient);
                criticP = clamp(criticP, 1.0e-6, 50.0);

                double targetK = problem.gamma * problem.b * criticP * problem.a
                    / (problem.r
                       + problem.gamma * problem.b * problem.b * criticP);
                actorK += actorRate * (targetK - actorK);
                actorK = clamp(actorK, 0.0, 3.0);
                x = xNext;
            }

            if (episode % 50 == 0 || episode == episodes - 1) {
                history.add(new Snapshot(episode, criticP, actorK));
            }
        }

        return new TrainingResult(criticP, actorK, history);
    }

    private static RolloutResult rolloutCost(
        Problem problem,
        double gain,
        double x0,
        int horizon) {

        double x = x0;
        double total = 0.0;
        double discount = 1.0;
        List<Double> trajectory = new ArrayList<>();
        trajectory.add(x);

        for (int step = 0; step < horizon; step++) {
            double u = -gain * x;
            total += discount * (problem.q * x * x + problem.r * u * u);
            x = problem.a * x + problem.b * u;
            trajectory.add(x);
            discount *= problem.gamma;
        }

        return new RolloutResult(total, trajectory);
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        Problem problem = new Problem();
        double[] reference = discountedRiccatiReference(problem);
        TrainingResult learned = trainAdaptiveCritic(problem, 500, 80, 7L);
        RolloutResult initialRollout = rolloutCost(problem, 0.80, 2.0, 60);
        RolloutResult learnedRollout = rolloutCost(problem, learned.actorK, 2.0, 60);

        System.out.println("episode, critic_p, actor_k");
        for (Snapshot snapshot : learned.history) {
            System.out.printf(
                "%4d, %10.6f, %10.6f%n",
                snapshot.episode,
                snapshot.criticP,
                snapshot.actorK);
        }

        System.out.println("\nReference discounted-Riccati solution");
        System.out.printf("p* = %.9f%n", reference[0]);
        System.out.printf("k* = %.9f%n", reference[1]);
        System.out.println("\nLearned adaptive-critic solution");
        System.out.printf("p_hat = %.9f%n", learned.criticP);
        System.out.printf("k_hat = %.9f%n", learned.actorK);
        System.out.printf(
            "closed-loop pole = %.9f%n",
            problem.a - problem.b * learned.actorK);
        System.out.printf(
            "initial-policy rollout cost = %.9f%n",
            initialRollout.cost);
        System.out.printf(
            "learned-policy rollout cost = %.9f%n",
            learnedRollout.cost);
        System.out.println("first ten learned-policy states:");
        for (int index = 0; index < Math.min(10, learnedRollout.trajectory.size()); index++) {
            System.out.printf(
                "%.6f%s",
                learnedRollout.trajectory.get(index),
                index == 9 ? System.lineSeparator() : " ");
        }
    }
}

13. MATLAB Implementation — Chapter28_Lesson1.m

The MATLAB script implements the critic, actor, Riccati reference, and rollout directly. Reinforcement Learning Toolbox supports actor and critic representations, Deep Learning Toolbox supports neural approximators, and Control System Toolbox provides linear-system analysis used to check the learned closed-loop pole.

% Chapter 28, Lesson 1: adaptive critic for a scalar discounted LQR.
% The critic estimates V(x) = p*x^2 from TD data. The actor uses the
% scalar greedy policy-improvement formula u = -k*x.

clear; clc;
rng(7, 'twister');

% Plant and discounted cost.
a = 1.05;
b = 0.50;
q = 1.00;
r = 0.10;
gamma = 0.98;

% Reference discounted Riccati iteration.
pStar = q;
for iteration = 1:100000
    denominator = r + gamma*b^2*pStar;
    pNext = q + gamma*a^2*pStar - (gamma*a*b*pStar)^2/denominator;
    if abs(pNext - pStar) < 1.0e-14
        pStar = pNext;
        break;
    end
    pStar = pNext;
end
kStar = gamma*b*pStar*a/(r + gamma*b^2*pStar);

% Adaptive-critic initialization.
criticP = 0.50;
actorK = 0.80;  % admissible: abs(a - b*actorK) < 1
criticRate = 0.08;
actorRate = 0.08;
episodes = 500;
stepsPerEpisode = 80;
history = zeros(11, 3);
historyIndex = 1;

for episode = 0:(episodes - 1)
    x = -2.0 + 4.0*rand();
    explorationScale = 0.25*exp(-episode/180.0);

    for step = 1:stepsPerEpisode
        exploration = explorationScale*randn();
        u = -actorK*x + exploration;
        xNext = a*x + b*u;
        stageCost = q*x^2 + r*u^2;

        delta = stageCost + gamma*criticP*xNext^2 - criticP*x^2;
        criticGradient = gamma*xNext^2 - x^2;
        criticP = criticP - criticRate*delta*criticGradient ...
            /(1.0 + criticGradient^2);
        criticP = min(max(criticP, 1.0e-6), 50.0);

        targetK = gamma*b*criticP*a/(r + gamma*b^2*criticP);
        actorK = actorK + actorRate*(targetK - actorK);
        actorK = min(max(actorK, 0.0), 3.0);
        x = xNext;
    end

    if mod(episode, 50) == 0 || episode == episodes - 1
        history(historyIndex, :) = [episode, criticP, actorK];
        historyIndex = historyIndex + 1;
    end
end
history = history(1:(historyIndex - 1), :);

% Deterministic rollout comparison.
x0 = 2.0;
horizon = 60;
[initialCost, ~] = rolloutCost(a, b, q, r, gamma, 0.80, x0, horizon);
[learnedCost, learnedTrajectory] = rolloutCost( ...
    a, b, q, r, gamma, actorK, x0, horizon);

fprintf('episode, critic_p, actor_k\n');
for row = 1:size(history, 1)
    fprintf('%4d, %10.6f, %10.6f\n', ...
        history(row, 1), history(row, 2), history(row, 3));
end

fprintf('\nReference discounted-Riccati solution\n');
fprintf('p* = %.9f\n', pStar);
fprintf('k* = %.9f\n', kStar);
fprintf('\nLearned adaptive-critic solution\n');
fprintf('p_hat = %.9f\n', criticP);
fprintf('k_hat = %.9f\n', actorK);
fprintf('closed-loop pole = %.9f\n', a - b*actorK);
fprintf('initial-policy rollout cost = %.9f\n', initialCost);
fprintf('learned-policy rollout cost = %.9f\n', learnedCost);
fprintf('first ten learned-policy states:\n');
fprintf('%.6f ', learnedTrajectory(1:10));
fprintf('\n');

figure('Name', 'Chapter28 Lesson1 Adaptive Critic');
plot(0:horizon, learnedTrajectory, 'LineWidth', 1.5);
grid on;
xlabel('Step');
ylabel('State x_k');
title('Learned-policy closed-loop state');

function [total, trajectory] = rolloutCost(a, b, q, r, gamma, gain, x0, horizon)
    x = x0;
    total = 0.0;
    discount = 1.0;
    trajectory = zeros(horizon + 1, 1);
    trajectory(1) = x;

    for step = 1:horizon
        u = -gain*x;
        total = total + discount*(q*x^2 + r*u^2);
        x = a*x + b*u;
        trajectory(step + 1) = x;
        discount = discount*gamma;
    end
end

14. Simulink Deployment Model — Chapter28_Lesson1_Simulink.m

This script programmatically creates a discrete-time Simulink model of the learned actor and plant. The online learning loop remains in the MATLAB script so students can inspect every update. In a production design, the actor, critic, replay/data buffer, projection logic, safety supervisor, and plant interface would be separate subsystems with explicit rates and saturation handling.

% Chapter 28, Lesson 1: programmatically create a Simulink rollout model.
% This model uses the learned gain generated by Chapter28_Lesson1.m.
% It isolates the deployed actor and plant; online critic learning remains in
% the MATLAB training script because variable-step algebraic learning loops
% are clearer there for a first university treatment.

clear; clc;
modelName = 'Chapter28_Lesson1_ActorDeployment';

% Learned gain from the deterministic seeded training experiment.
a = 1.05;
b = 0.50;
kLearned = 1.60914;
sampleTime = 1.0;
stopTime = 20;

if bdIsLoaded(modelName)
    close_system(modelName, 0);
end
if exist([modelName '.slx'], 'file')
    delete([modelName '.slx']);
end

new_system(modelName);
open_system(modelName);
set_param(modelName, 'Solver', 'FixedStepDiscrete');
set_param(modelName, 'FixedStep', num2str(sampleTime));
set_param(modelName, 'StopTime', num2str(stopTime));

add_block('simulink/Sources/Constant', [modelName '/Zero'], ...
    'Value', '0', 'Position', [40 45 80 75]);
add_block('simulink/Math Operations/Gain', [modelName '/Actor_-K'], ...
    'Gain', num2str(-kLearned), 'Position', [180 110 260 150]);
add_block('simulink/Math Operations/Gain', [modelName '/Plant_B'], ...
    'Gain', num2str(b), 'Position', [320 110 390 150]);
add_block('simulink/Math Operations/Gain', [modelName '/Plant_A'], ...
    'Gain', num2str(a), 'Position', [320 210 390 250]);
add_block('simulink/Math Operations/Sum', [modelName '/State_Update'], ...
    'Inputs', '++', 'Position', [450 145 480 215]);
add_block('simulink/Discrete/Unit Delay', [modelName '/State_x'], ...
    'InitialCondition', '2', 'SampleTime', num2str(sampleTime), ...
    'Position', [540 160 610 200]);
add_block('simulink/Sinks/Scope', [modelName '/State_Scope'], ...
    'Position', [700 150 760 210]);
add_block('simulink/Sinks/To Workspace', [modelName '/State_Log'], ...
    'VariableName', 'xSim', 'SaveFormat', 'Structure With Time', ...
    'Position', [690 245 790 275]);

add_line(modelName, 'State_x/1', 'Actor_-K/1', 'autorouting', 'on');
add_line(modelName, 'Actor_-K/1', 'Plant_B/1', 'autorouting', 'on');
add_line(modelName, 'Plant_B/1', 'State_Update/1', 'autorouting', 'on');
add_line(modelName, 'State_x/1', 'Plant_A/1', 'autorouting', 'on');
add_line(modelName, 'Plant_A/1', 'State_Update/2', 'autorouting', 'on');
add_line(modelName, 'State_Update/1', 'State_x/1', 'autorouting', 'on');
add_line(modelName, 'State_x/1', 'State_Scope/1', 'autorouting', 'on');
add_line(modelName, 'State_x/1', 'State_Log/1', 'autorouting', 'on');

save_system(modelName);
simulationOutput = sim(modelName);
assignin('base', 'Chapter28_Lesson1_SimulationOutput', simulationOutput);

fprintf('Created and simulated %s.slx\n', modelName);
fprintf('Closed-loop pole: %.6f\n', a - b*kLearned);

15. Wolfram Mathematica Implementation — Chapter28_Lesson1.nb

The Wolfram Language implementation uses associations, modules, and symbolic-numeric iteration. For larger studies, FindMinimum, NDSolve, LinearSolve, and the neural-network functions NetChain and NetTrain can support actor- critic approximation and HJB-residual minimization.

Notebook[{
  Cell["Chapter 28, Lesson 1: Adaptive Critics and Approximate Dynamic Programming", "Title"],
  Cell["Quadratic critic and actor implementation for a scalar discounted LQR.", "Text"],
  Cell[BoxData[
    ToBoxes[
      Defer[
        ClearAll["Global`*"];
        SeedRandom[7];
        
        problem = <|"a" -> 1.05, "b" -> 0.50, "q" -> 1.00,
           "r" -> 0.10, "gamma" -> 0.98|>;
        
        DiscountedRiccatiReference[pr_Association] := Module[
          {p = pr["q"], pNext, denominator, k, iteration},
          For[iteration = 1, iteration <= 100000, iteration++,
           denominator = pr["r"] + pr["gamma"] pr["b"]^2 p;
           pNext = pr["q"] + pr["gamma"] pr["a"]^2 p -
             (pr["gamma"] pr["a"] pr["b"] p)^2/denominator;
           If[Abs[pNext - p] < 10^-14, p = pNext; Break[]];
           p = pNext;
          ];
          k = pr["gamma"] pr["b"] p pr["a"]/
            (pr["r"] + pr["gamma"] pr["b"]^2 p);
          {p, k}
        ];
        
        TrainAdaptiveCritic[pr_Association, episodes_Integer : 500,
           stepsPerEpisode_Integer : 80] := Module[
          {criticP = 0.50, actorK = 0.80, criticRate = 0.08,
           actorRate = 0.08, x, explorationScale, exploration, u, xNext,
           stageCost, delta, criticGradient, targetK, history = {}, episode,
           step},
        
          For[episode = 0, episode < episodes, episode++,
           x = RandomReal[{-2.0, 2.0}];
           explorationScale = 0.25 Exp[-episode/180.0];
        
           For[step = 1, step <= stepsPerEpisode, step++,
            exploration = explorationScale RandomVariate[NormalDistribution[]];
            u = -actorK x + exploration;
            xNext = pr["a"] x + pr["b"] u;
            stageCost = pr["q"] x^2 + pr["r"] u^2;
        
            delta = stageCost + pr["gamma"] criticP xNext^2 - criticP x^2;
            criticGradient = pr["gamma"] xNext^2 - x^2;
            criticP = criticP - criticRate delta criticGradient/
               (1.0 + criticGradient^2);
            criticP = Clip[criticP, {10^-6, 50.0}];
        
            targetK = pr["gamma"] pr["b"] criticP pr["a"]/
              (pr["r"] + pr["gamma"] pr["b"]^2 criticP);
            actorK = actorK + actorRate (targetK - actorK);
            actorK = Clip[actorK, {0.0, 3.0}];
            x = xNext;
           ];
        
           If[Mod[episode, 50] == 0 || episode == episodes - 1,
            AppendTo[history, {episode, criticP, actorK}]
           ];
          ];
          <|"CriticP" -> criticP, "ActorK" -> actorK, "History" -> history|>
        ];
        
        RolloutCost[pr_Association, gain_?NumericQ, x0_ : 2.0,
           horizon_Integer : 60] := Module[
          {x = x0, total = 0.0, discount = 1.0, trajectory = {x0}, u,
           step},
          For[step = 1, step <= horizon, step++,
           u = -gain x;
           total = total + discount (pr["q"] x^2 + pr["r"] u^2);
           x = pr["a"] x + pr["b"] u;
           AppendTo[trajectory, x];
           discount = discount pr["gamma"];
          ];
          {total, trajectory}
        ];
        
        {pStar, kStar} = DiscountedRiccatiReference[problem];
        learned = TrainAdaptiveCritic[problem];
        {initialCost, initialTrajectory} = RolloutCost[problem, 0.80];
        {learnedCost, learnedTrajectory} = RolloutCost[problem, learned["ActorK"]];
        
        Print["episode, critic_p, actor_k"];
        Scan[Print[NumberForm[#, {12, 6}]] &, learned["History"]];
        Print["Reference p* = ", NumberForm[pStar, {12, 9}]];
        Print["Reference k* = ", NumberForm[kStar, {12, 9}]];
        Print["Learned p_hat = ", NumberForm[learned["CriticP"], {12, 9}]];
        Print["Learned k_hat = ", NumberForm[learned["ActorK"], {12, 9}]];
        Print["Closed-loop pole = ",
         NumberForm[problem["a"] - problem["b"] learned["ActorK"], {12, 9}]];
        Print["Initial-policy rollout cost = ", NumberForm[initialCost, {12, 9}]];
        Print["Learned-policy rollout cost = ", NumberForm[learnedCost, {12, 9}]];
        
        ListLinePlot[learnedTrajectory, PlotRange -> All,
         AxesLabel -> {"step", "x_k"},
         PlotLabel -> "Learned-policy closed-loop state"]
      ]
    ]
  ], "Input"]
},
WindowTitle -> "Chapter28_Lesson1",
StyleDefinitions -> "Default.nb"]

16. Interpreting the Numerical Result

For the selected problem, discounted Riccati iteration produces

\[ p^{\star}\approx1.337911639, \qquad k^{\star}\approx1.609103043. \]

A representative seeded Python run produces a critic coefficient and actor gain within approximately \(10^{-4}\) of these reference values. The learned closed-loop pole is approximately \(0.2454\), so the deterministic state decays much faster than under the initial pole \(0.65\). The initial gain is stabilizing but not optimal; the critic supplies the long-horizon information needed to trade state decay against control effort.

Agreement in this example does not establish a general convergence theorem for nonlinear neural actor-critic control. It occurs because the value class contains the exact scalar LQR value, the actor class contains the exact linear policy, the plant derivative used in improvement is correct, the initial policy is admissible, and exploration visits a useful range of states.

17. When a Neural Critic Is Necessary

A quadratic basis is exact for LQR but generally inadequate for nonlinear plants, constraints, hybrid effects, and nonquadratic objectives. A neural critic replaces fixed features by learned nonlinear features:

\[ \widehat V(x;w_c) =W_L\sigma_{L-1}\!\left( W_{L-1}\sigma_{L-2}(\cdots\sigma_1(W_1x+b_1)\cdots)+b_{L-1} \right)+b_L. \]

Universal approximation results justify representational capacity on compact sets, but they do not guarantee that online optimization finds suitable weights, that the Bellman residual is small everywhere, or that the induced actor is stable. Neural adaptive control therefore needs more than a large network: data coverage, regularization, bounded updates, robust supervision, and a stability argument remain central.

This lesson is intentionally a pointer. Deep-network architectures, replay buffers, target critics, stochastic policy gradients, entropy regularization, distributional critics, and offline reinforcement learning belong in Neural Control, Intelligent Control, or Learning-Based Control courses. The transferable foundation is the Bellman/HJB structure developed here.

18. Problems and Solutions

Problem 1 — Bellman contraction. Let \((TV)(x)=\min_u\{\ell(x,u)+\gamma V(f(x,u))\}\) with \(0<\gamma<1\). Prove that \(T\) is a contraction under the sup norm.

Solution. For every state, use the inequality

\[ |\min_u a_u-\min_u b_u|\le\max_u|a_u-b_u|. \]

Set

\[ a_u=\ell(x,u)+\gamma V(f(x,u)), \qquad b_u=\ell(x,u)+\gamma W(f(x,u)). \]

Then

\[ \begin{aligned} |(TV)(x)-(TW)(x)| &\le\gamma\max_u|V(f(x,u))-W(f(x,u))|\\ &\le\gamma\|V-W\|_{\infty}. \end{aligned} \]

Taking the supremum gives \(\|TV-TW\|_{\infty}\le\gamma\|V-W\|_{\infty}\). Since \(\gamma<1\), the operator is a contraction and has a unique fixed point.

Problem 2 — Derive the scalar greedy actor. For \(x^{+}=ax+bu\) and \(\widehat V(x)=px^2\), derive the minimizing control for \(qx^2+ru^2+\gamma p(ax+bu)^2\).

Solution. Differentiate with respect to the scalar input:

\[ \frac{\partial}{\partial u} \left[qx^2+ru^2+\gamma p(ax+bu)^2\right] =2ru+2\gamma pb(ax+bu). \]

Set the derivative to zero and collect terms in \(u\):

\[ (r+\gamma pb^2)u=-\gamma pba x. \]

\[ u=-\frac{\gamma pba}{r+\gamma pb^2}x=-k(p)x. \]

Since \(r>0\) and \(p\ge0\), the second derivative \(2(r+\gamma pb^2)\) is positive, so this stationary point is the unique minimizer.

Problem 3 — One critic/actor update. Use \(a=1.05\), \(b=0.5\), \(q=1\), \(r=0.1\), \(\gamma=0.98\), \(x=1\), \(p_c=0.5\), \(k_a=0.8\), no exploration, and \(\alpha_c=\alpha_a=0.08\). Compute one update.

Solution. The input and next state are

\[ u=-0.8, \qquad x^{+}=1.05+0.5(-0.8)=0.65. \]

The stage cost and TD residual are

\[ c=1^2+0.1(0.8)^2=1.064, \]

\[ \delta =1.064+0.98(0.5)(0.65)^2-0.5 =0.771025. \]

The critic residual gradient is

\[ g_c=0.98(0.65)^2-1=-0.58595. \]

Hence

\[ p_c^{+} =0.5-0.08(0.771025) \frac{-0.58595}{1+(-0.58595)^2} \approx0.526905. \]

The greedy target and actor update are

\[ k_{\mathrm{target}} =\frac{0.98(0.5)(0.526905)(1.05)} {0.1+0.98(0.5)^2(0.526905)} \approx1.183337, \]

\[ k_a^{+}=0.8+0.08(1.183337-0.8)\approx0.830667. \]

Problem 4 — TD normal equations. For a fixed policy and linear critic \(\widehat V(x)=w^{\mathsf T}\phi(x)\), derive the expected equation obtained by requiring the TD residual to be orthogonal to \(\phi(x_k)\).

Solution. The orthogonality condition is

\[ \mathbb{E}[\phi(x_k)\delta_k]=0. \]

Substitute

\[ \delta_k=c_k+\gamma w^{\mathsf T}\phi(x_{k+1}) -w^{\mathsf T}\phi(x_k). \]

Rearrangement gives

\[ \mathbb{E}\!\left[ \phi(x_k)(\phi(x_k)-\gamma\phi(x_{k+1}))^{\mathsf T} \right]w =\mathbb{E}[\phi(x_k)c_k]. \]

Thus \(Aw=b\). Nonsingularity and conditioning of \(A\) depend on the features and the transition distribution generated by the behavior policy.

Problem 5 — Why approximate policy improvement can fail. Suppose \(\widehat V=V^{\mu}+e\). Explain why the greedy policy based on \(\widehat V\) need not improve the true value, even if \(\|e\|_{\infty}\) is small.

Solution. The actor compares actions through

\[ \ell(x,u)+\gamma\widehat V(f(x,u)) =Q^{\mu}(x,u)+\gamma e(f(x,u)). \]

If two candidate actions have nearly equal true action values, even a small approximation error can reverse their ordering. The selected policy can then move the state toward regions where the critic has poorer data, magnifying subsequent error. A uniform error bound can produce an approximate performance bound, but exact monotonic improvement requires stronger conditions than a small average training loss.

Problem 6 — Model-free quadratic Q critic. Let \(Q(x,u)=h_{xx}x^2+2h_{xu}xu+h_{uu}u^2\) with \(h_{uu}>0\). Derive the greedy control and explain why the plant coefficients need not appear explicitly.

Solution. Differentiate with respect to the action:

\[ \frac{\partial Q}{\partial u}=2h_{xu}x+2h_{uu}u. \]

Setting the derivative to zero gives

\[ u^{+}=-h_{uu}^{-1}h_{ux}x. \]

The coefficients \(h_{xu}\) and \(h_{uu}\) can be estimated from transition data using a Q-Bellman residual. Once they are estimated, greedy improvement uses only the critic coefficients; explicit values of \(a\) and \(b\) are unnecessary. Adequate state-action exploration is still required to identify the quadratic terms.

19. Summary

Adaptive critics approximate value information rather than only plant or controller parameters. Bellman and HJB equations define the target; temporal-difference residuals provide data-driven learning signals; policy iteration connects critic evaluation to actor improvement. Exact dynamic programming has contraction and policy-improvement guarantees, but approximation, limited data, nonlinear optimization, exploration, and policy-induced distribution shift complicate those guarantees. The scalar LQR laboratory showed a transparent actor-critic mechanism whose learned quadratic value and feedback gain approach the discounted Riccati solution. Neural networks extend representational power, but stability, safety, and data coverage remain control-design obligations.

20. References

  1. Barto, A.G., Sutton, R.S., & Anderson, C.W. (1983). Neuronlike adaptive elements that can solve difficult learning control problems. IEEE Transactions on Systems, Man, and Cybernetics, SMC-13(5), 834–846. doi:10.1109/TSMC.1983.6313077.
  2. Sutton, R.S. (1988). Learning to predict by the methods of temporal differences. Machine Learning, 3, 9–44. doi:10.1007/BF00115009.
  3. Watkins, C.J.C.H., & Dayan, P. (1992). Q-learning. Machine Learning, 8, 279–292. doi:10.1007/BF00992698.
  4. Bradtke, S.J., & Barto, A.G. (1996). Linear least-squares algorithms for temporal difference learning. Machine Learning, 22, 33–57. doi:10.1007/BF00114723.
  5. Tsitsiklis, J.N., & Van Roy, B. (1997). An analysis of temporal-difference learning with function approximation. IEEE Transactions on Automatic Control, 42(5), 674–690. doi:10.1109/9.580874.
  6. Doya, K. (2000). Reinforcement learning in continuous time and space. Neural Computation, 12(1), 219–245. doi:10.1162/089976600300015961.
  7. Abu-Khalaf, M., & Lewis, F.L. (2005). Nearly optimal control laws for nonlinear systems with saturating actuators using a neural network HJB approach. Automatica, 41(5), 779–791. doi:10.1016/j.automatica.2004.11.034.
  8. Vamvoudakis, K.G., & Lewis, F.L. (2010). Online actor–critic algorithm to solve the continuous-time infinite-horizon optimal control problem. Automatica, 46(5), 878–888. doi:10.1016/j.automatica.2010.02.018.
Support CaaT Academy

Help keep these engineering tutorials free and growing

If these lessons, examples, and project pages help you, a small donation supports the continued creation and improvement of free control, robotics, software, and engineering education resources.

Created and maintained by Abolfazl Mohammadijoo.