Consider the outlet-temperature loop of a steam-heated process unit. The manipulated variable \( u(k) \) is the steam-valve command in percent, the controlled variable \( y(k) \) is outlet temperature, and \( r(k) \) is the operator or supervisory-layer setpoint. Heat-transfer gain decreases as fouling develops, while effective thermal inertia increases. A fixed PI controller can remain stable but become slow, oscillatory, or unnecessarily aggressive across the full operating campaign.

The adaptive deployment must satisfy five non-negotiable requirements:

1. Never remove the validated PI controller; retain it as an immediate fallback.
2. Estimate in shadow mode before the adaptive controller receives authority.
3. Constrain estimates and control signals to physically credible regions.
4. Record enough evidence to reproduce every mode change and parameter update.
5. Demonstrate performance improvement without violating safety or production limits.

flowchart TD
  R["Setpoint from operator or supervisory layer"] --> SUP["Supervisory state machine"]
  Y["Validated temperature \nmeasurement"] --> PI["Backup PI controller"]
  Y --> ID["Projected RLS estimator"]
  ID --> AD["Certainty-equivalent \nadaptive controller"]
  PI --> BLEND["Bumpless authority blend"]
  AD --> BLEND
  SUP --> ID
  SUP --> BLEND
  BLEND --> LIM["Rate limit and saturation"]
  LIM --> P["Steam-heated process"]
  P --> Y
  SUP --> HIST["Historian, alarms, \nand audit record"]
  ID --> HIST
  LIM --> HIST
        

Let \( T_a \) denote a measured or slowly varying ambient/feed reference temperature, and define the deviation state \( x(k)=y(k)-T_a \). A sampled first-order model is

\[ x(k+1)=a(k)x(k)+b(k)u(k)+d(k)+v(k), \]

where \( 0<a(k)<1 \) is the discrete thermal-memory coefficient, \( b(k)>0 \) is the heater gain, \( d(k) \) is a bounded load disturbance, and \( v(k) \) collects process and discretization errors. If a continuous first-order model has time constant \( \tau \), static gain \( K \), and sample period \( T_s \), then

\[ a=e^{-T_s/\tau}, \qquad b=K\left(1-e^{-T_s/\tau}\right). \]

Fouling typically produces \( \tau\uparrow \) and \( K\downarrow \); therefore \( a \) moves closer to one while \( b \) decreases. The simulation changes the true parameters from \( (a,b)=(0.940,0.035) \) to \( (0.965,0.025) \) during operation.

The model is linear in unknown parameters. With \( \boldsymbol{\theta}=[a\;b]^T \) and \( \boldsymbol{\phi}(k)=[x(k)\;u(k)]^T \),

\[ x(k+1)=\boldsymbol{\theta}^T(k)\boldsymbol{\phi}(k)+d(k)+v(k). \]

This low-order model is intentionally a controller model, not a claim that all exchanger physics are first order. Unmodeled valve dynamics, transport delay, nonlinear steam enthalpy, and sensor lag are handled through conservative bandwidth, robust gating, and validation rather than by pretending they do not exist.

The one-step prediction and innovation are

\[ \widehat{x}(k+1|k)=\widehat{\boldsymbol{\theta} }^T(k)\boldsymbol{\phi}(k), \qquad \varepsilon(k+1)=x(k+1)-\widehat{x}(k+1|k). \]

Exponentially weighted RLS uses

\[ \mathbf{K}(k)=\frac{\mathbf{P}(k)\boldsymbol{\phi}(k)} {\lambda+\boldsymbol{\phi}^T(k)\mathbf{P}(k)\boldsymbol{\phi}(k)}, \]

\[ \widehat{\boldsymbol{\theta} }_u(k+1)= \widehat{\boldsymbol{\theta} }(k)+\mathbf{K}(k)\varepsilon(k+1), \]

\[ \mathbf{P}(k+1)=\frac{1}{\lambda} \left[\mathbf{P}(k)-\mathbf{K}(k)\boldsymbol{\phi}^T(k)\mathbf{P}(k)\right], \qquad 0<\lambda\le 1. \]

The unconstrained estimate \( \widehat{\boldsymbol{\theta} }_u \) is projected into the industrially approved set

\[ \Omega = \left\{ (a,b) : a_{\min} \le a \le a_{\max}, \; b_{\min} \le b \le b_{\max} \right\}, \qquad \widehat{\boldsymbol{\theta}}(k+1) = \operatorname{Proj}_{\Omega}\!\left( \widehat{\boldsymbol{\theta}}_{u}(k+1) \right). \]

Projection prevents a temporary data defect from creating an unstable pole estimate or a near-zero heater-gain estimate. In the case study, \( 0.80\le\widehat{a}\le0.995 \) and \( 0.005\le\widehat{b}\le0.10 \).

3.1 Update gating

Estimation is frozen when the measurement-quality bit is bad, the command is at a hard limit, the innovation is physically implausible, the loop is in maintenance/manual mode, or a safety interlock is active. Freezing the estimator does not mean setting \( \mathbf{P} \) to zero. Mild covariance inflation can preserve the ability to relearn after a gated interval.

3.2 Excitation and identifiability

Good tracking alone does not imply parameter convergence. In steady closed-loop operation, \( x(k) \) and \( u(k) \) can become nearly collinear. The information matrix

\[ \mathbf{G}_N=\sum_{k=0}^{N-1} \boldsymbol{\phi}(k)\boldsymbol{\phi}^T(k) \]

must remain sufficiently well conditioned for both parameters to be identified. The example injects a small bounded square-wave probe only when tracking error is small and the actuator is far from its limits. In a real facility, such excitation requires an approved operating window and a production-impact assessment.

Select a desired first-order reference model

\[ x_m(k+1)=p_m x_m(k)+(1-p_m)x_r(k), \qquad x_r(k)=r(k)-T_a, \qquad 0<p_m<1. \]

Enforcing the same one-step relation on the estimated plant yields the raw adaptive command

\[ u_a(k)=\frac{\left(p_m-\widehat{a}(k)\right)x(k)+ (1-p_m)x_r(k)}{\widehat{b}(k)}. \]

To avoid division by an unsafe gain estimate, the denominator is protected by projection and by a software lower bound. The deployed command is not \( u_a \) directly:

\[ u_c(k)=\operatorname{sat}_{[u_{\min},u_{\max}]} \left(\operatorname{rate}_{\Delta u_{\max} } \left((1-\alpha(k))u_{PI}(k)+\alpha(k)u_a(k)\right)\right), \]

where \( \alpha(k) \) ramps from zero to one during bumpless transfer. The rate limiter protects the valve, steam header, and process from abrupt commands.

Assume temporarily that the reference is constant, saturation is inactive, and the true parameters are known. Substituting the ideal command into the plant gives

\[ x(k+1)=p_m x(k)+(1-p_m)x_r. \]

For the tracking error \( e(k)=x(k)-x_r \),

\[ e(k+1)=p_m e(k), \qquad e(k)=p_m^k e(0). \]

Hence the ideal loop is exponentially stable for \( |p_m|<1 \). With parameter error, disturbance, model mismatch, and command modification, write

\[ e(k+1)=p_m e(k)+w_e(k), \]

\[ w_e(k)=\widetilde{\boldsymbol{\theta} }^T(k)\boldsymbol{\phi}(k) +d(k)+v(k)+b(k)\Delta u(k)+\Delta_r(k), \]

where \( \Delta u(k)=u_c(k)-u_a(k) \) captures rate limiting, saturation, and blending, while \( \Delta_r(k) \) captures reference variation. Iteration gives

\[ |e(k)|\le |p_m|^k|e(0)|+ \sum_{i=0}^{k-1}|p_m|^{k-1-i}|w_e(i)|. \]

If \( |w_e(k)|\le\bar{w} \), then

\[ \limsup_{k\to\infty}|e(k)|\le \frac{\bar{w} }{1-|p_m|}. \]

This input-to-state bound explains an important deployment trade-off. Choosing \( p_m \) very close to one makes the nominal response slow and can enlarge the disturbance-to-error bound; choosing it too small demands aggressive valve motion and excites unmodeled dynamics. Industrial tuning therefore begins from a conservative closed-loop time constant and is tightened only after staged tests.

A deployable adaptive loop should have explicit states rather than a single enable bit:

1. Baseline: PI controls; estimator may be disabled.
2. Shadow: PI controls; estimator runs and is evaluated offline.
3. Blend: authority transfers gradually from PI to adaptive control.
4. Adaptive: adaptive command has authority under all limiters.
5. Hold: estimates freeze; last valid or backup control is used.
6. Fallback: authority returns to PI after a trip or failed monitor.

Typical invariants checked every sample are

\[ u_{\min}\le u_c(k)\le u_{\max}, \qquad |u_c(k)-u_c(k-1)|\le\Delta u_{\max}, \]

\[ \widehat{\boldsymbol{\theta} }(k)\in\Omega, \qquad |\varepsilon(k)|\le\varepsilon_{\max} \quad \text{when adaptation is enabled}. \]

A fallback trigger may combine hard conditions and persistence timers. Examples include bad sensor quality, repeated projection-bound contact, covariance loss of positive definiteness, innovation variance above its baseline envelope, excessive valve travel, tracking error above a limit for a specified duration, or disagreement with a parallel safety model.

flowchart TD
  A["Freeze requirements and hazard review"] --> B["Validate sensors, actuator limits, and fallback PI"]
  B --> C["Replay historian data \nand digital-twin tests"]
  C --> D["Run estimator in shadow mode"]
  D --> E{"Prediction and parameter \nchecks pass?"}
  E -->|"no"| F["Retune model, filters, \ngates, or bounds"]
  F --> C
  E -->|"yes"| G["Limited-authority plant trial"]
  G --> H["Bumpless blend under \noperator supervision"]
  H --> I{"Acceptance metrics and \nsafety monitors pass?"}
  I -->|"no"| J["Automatic fallback \nand root-cause review"]
  J --> C
  I -->|"yes"| K["Controlled production release"]
  K --> L["Periodic audit, drift review, \nand rollback test"]
        

7.1 Offline evidence

Historian replay should include normal campaigns, grade or load changes, shutdown/restart segments, sensor bias, frozen measurements, valve stiction, saturation, communication delay, and parameter drift. The purpose is not only to obtain low average error; it is to expose failure modes and verify deterministic fallback.

7.2 Shadow-mode acceptance

During shadow mode, compare one-step prediction error, parameter trajectories, projection events, covariance conditioning, and the hypothetical adaptive command against the actual PI command. Shadow operation should cover the intended operating envelope before authority is granted.

7.3 Limited-authority trial

Initial trials should constrain adaptive authority, setpoint range, operating grade, rate limits, and duration. An operator-visible mode indication and one-action rollback are required. Every transition should be time-stamped with the software version, tuning set, model bounds, and reason code.

Performance metrics should be computed over declared windows:

\[ \operatorname{IAE}=\sum_{k=k_0}^{k_1}|e(k)|T_s, \qquad \operatorname{RMSE}=\sqrt{\frac{1}{N} \sum_{k=k_0}^{k_1}e^2(k)}, \]

\[ \operatorname{TV}(u)=\sum_{k=k_0+1}^{k_1}|u(k)-u(k-1)|, \qquad \rho_{sat}=\frac{1}{N}\sum_{k=k_0}^{k_1} \mathbf{1}_{\{u(k)\text{ modified by a hard limit}\} }. \]

A release decision should compare adaptive and baseline operation using comparable production conditions. At minimum, archive setpoint, measurement, command before and after limiting, active controller mode, parameter estimates, covariance diagnostics, innovation, gate status, projection flags, fallback causes, operator actions, and software/configuration hashes.

Parameter convergence is not itself the primary business objective. A controller can achieve safe and accurate tracking while estimates remain biased because of weak excitation or model mismatch. Conversely, apparently plausible estimates do not justify release if valve wear, overshoot, or fallback frequency is unacceptable.

The Python version uses NumPy for matrix operations and Matplotlib for the acceptance plots. It writes a CSV audit trace and a PNG.

Chapter27_Lesson5.py

"""Chapter27_Lesson5.py
Industrial adaptive-control case study for a steam-heated thermal process.

Dependencies:
    pip install numpy matplotlib

The simulation demonstrates:
  * baseline PI operation and shadow-mode identification,
  * projected RLS with exponential forgetting,
  * certainty-equivalent adaptive control,
  * bumpless blending, rate limits, saturation, probing, and adaptation gating,
  * a plant change representing heat-transfer fouling.
"""

from __future__ import annotations

import csv
import math
from dataclasses import dataclass

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class Limits:
    u_min: float = 0.0
    u_max: float = 100.0
    du_max: float = 3.0
    a_min: float = 0.80
    a_max: float = 0.995
    b_min: float = 0.005
    b_max: float = 0.10


def clamp(value: float, lower: float, upper: float) -> float:
    return max(lower, min(upper, value))


def setpoint(k: int) -> float:
    if k < 300:
        return 55.0
    if k < 600:
        return 65.0
    return 58.0


def true_parameters(k: int) -> tuple[float, float]:
    # At k = 450, fouling increases thermal inertia and reduces heater gain.
    if k < 450:
        return 0.940, 0.035
    return 0.965, 0.025


def projected_rls_update(
    theta: np.ndarray,
    covariance: np.ndarray,
    phi: np.ndarray,
    target: float,
    forgetting: float,
    limits: Limits,
) -> tuple[np.ndarray, np.ndarray, float]:
    prediction = float(theta @ phi)
    innovation = target - prediction
    denominator = forgetting + float(phi @ covariance @ phi)
    gain = covariance @ phi / denominator

    theta_new = theta + gain * innovation
    theta_new[0] = clamp(float(theta_new[0]), limits.a_min, limits.a_max)
    theta_new[1] = clamp(float(theta_new[1]), limits.b_min, limits.b_max)

    covariance_new = (
        covariance - np.outer(gain, phi) @ covariance
    ) / forgetting
    covariance_new = 0.5 * (covariance_new + covariance_new.T)

    # Numerical safeguard: keep the covariance symmetric positive definite
    # and prevent extreme eigenvalues after long periods of weak excitation.
    eigenvalues, eigenvectors = np.linalg.eigh(covariance_new)
    eigenvalues = np.clip(eigenvalues, 1.0e-5, 1.0e4)
    covariance_new = eigenvectors @ np.diag(eigenvalues) @ eigenvectors.T
    return theta_new, covariance_new, innovation


def run_simulation(seed: int = 27) -> dict[str, np.ndarray]:
    rng = np.random.default_rng(seed)
    limits = Limits()

    sample_count = 900
    sample_time = 1.0
    ambient_temperature = 25.0
    desired_time_constant = 20.0
    desired_pole = math.exp(-sample_time / desired_time_constant)
    forgetting = 0.995

    y = np.zeros(sample_count)
    measured_y = np.zeros(sample_count)
    u = np.zeros(sample_count)
    r = np.zeros(sample_count)
    a_true = np.zeros(sample_count)
    b_true = np.zeros(sample_count)
    a_hat = np.zeros(sample_count)
    b_hat = np.zeros(sample_count)
    innovation = np.zeros(sample_count)
    saturation = np.zeros(sample_count, dtype=bool)
    adaptation_enabled = np.zeros(sample_count, dtype=bool)
    adaptive_blend = np.zeros(sample_count)

    y[0] = ambient_temperature
    theta = np.array([0.900, 0.030], dtype=float)
    covariance = 500.0 * np.eye(2)
    a_hat[0], b_hat[0] = theta

    # Baseline PI controller used during commissioning and as fallback.
    kp = 2.0
    ki = 0.04
    pi_integral = 20.0
    previous_u = 40.0

    for k in range(sample_count - 1):
        r[k] = setpoint(k)
        a, b = true_parameters(k)
        a_true[k], b_true[k] = a, b

        measured_y[k] = y[k] + rng.normal(0.0, 0.08)
        tracking_error = r[k] - measured_y[k]

        # Backup PI with conditional integration anti-windup.
        integral_candidate = pi_integral + ki * tracking_error
        pi_raw = 20.0 + kp * tracking_error + integral_candidate
        pi_output = clamp(pi_raw, limits.u_min, limits.u_max)
        integration_is_safe = (
            abs(pi_output - pi_raw) < 1.0e-12
            or (pi_output >= limits.u_max and tracking_error < 0.0)
            or (pi_output <= limits.u_min and tracking_error > 0.0)
        )
        if integration_is_safe:
            pi_integral = integral_candidate

        # Certainty-equivalent adaptive command for the deviation model
        # x(k+1) = a*x(k) + b*u(k), x = y - ambient.
        x = measured_y[k] - ambient_temperature
        x_reference = r[k] - ambient_temperature
        effective_b = max(float(theta[1]), limits.b_min)
        adaptive_raw = (
            (desired_pole - float(theta[0])) * x
            + (1.0 - desired_pole) * x_reference
        ) / effective_b

        # Small bounded excitation is injected only near the setpoint and away
        # from actuator limits. In a plant this requires an approved test window.
        probe = 0.0
        if k >= 80 and abs(tracking_error) < 2.0 and 3.0 < adaptive_raw < 97.0:
            probe = 0.8 if (k // 15) % 2 == 0 else -0.8
        adaptive_raw += probe

        # Shadow mode for 80 samples, followed by a 40-sample bumpless blend.
        if k < 80:
            blend = 0.0
            desired_u = pi_output
        else:
            blend = min(1.0, (k - 80) / 40.0)
            desired_u = (1.0 - blend) * pi_output + blend * adaptive_raw
        adaptive_blend[k] = blend

        rate_limited_u = clamp(
            desired_u,
            previous_u - limits.du_max,
            previous_u + limits.du_max,
        )
        commanded_u = clamp(rate_limited_u, limits.u_min, limits.u_max)
        u[k] = commanded_u
        saturation[k] = (
            abs(commanded_u - desired_u) > 1.0e-9
            or commanded_u < 0.5
            or commanded_u > 99.5
        )

        # First-order thermal process. The negative pulse models a colder feed.
        feed_disturbance = -0.45 if 620 <= k < 645 else 0.0
        y[k + 1] = (
            a * y[k]
            + b * commanded_u
            + (1.0 - a) * ambient_temperature
            + feed_disturbance
            + rng.normal(0.0, 0.03)
        )

        # RLS predicts x(k+1) from [x(k), u(k)]. Update only with plausible
        # innovations and while the manipulated variable is not hard saturated.
        phi = np.array([measured_y[k] - ambient_temperature, commanded_u])
        target = y[k + 1] - ambient_temperature
        predicted_innovation = target - float(theta @ phi)
        estimator_is_valid = (
            k >= 20
            and 1.0 < commanded_u < 99.0
            and abs(predicted_innovation) < 5.0
        )

        if estimator_is_valid:
            theta, covariance, innovation[k] = projected_rls_update(
                theta, covariance, phi, target, forgetting, limits
            )
            adaptation_enabled[k] = True
        else:
            # Covariance inflation preserves sensitivity after gated intervals.
            covariance /= forgetting
            eigenvalues, eigenvectors = np.linalg.eigh(covariance)
            eigenvalues = np.clip(eigenvalues, 1.0e-5, 1.0e4)
            covariance = eigenvectors @ np.diag(eigenvalues) @ eigenvectors.T
            innovation[k] = predicted_innovation

        a_hat[k + 1], b_hat[k + 1] = theta
        previous_u = commanded_u

    r[-1] = r[-2]
    u[-1] = u[-2]
    measured_y[-1] = y[-1]
    a_true[-1], b_true[-1] = true_parameters(sample_count - 1)
    adaptive_blend[-1] = adaptive_blend[-2]

    return {
        "time": np.arange(sample_count) * sample_time,
        "setpoint": r,
        "temperature": y,
        "measured_temperature": measured_y,
        "heater_command": u,
        "a_true": a_true,
        "b_true": b_true,
        "a_hat": a_hat,
        "b_hat": b_hat,
        "innovation": innovation,
        "saturation": saturation.astype(float),
        "adaptation_enabled": adaptation_enabled.astype(float),
        "adaptive_blend": adaptive_blend,
    }


def save_results(results: dict[str, np.ndarray]) -> None:
    csv_name = "Chapter27_Lesson5_results.csv"
    keys = list(results.keys())
    with open(csv_name, "w", newline="", encoding="utf-8") as stream:
        writer = csv.writer(stream)
        writer.writerow(keys)
        writer.writerows(zip(*(results[key] for key in keys)))

    time = results["time"]
    fig, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=True)
    axes[0].plot(time, results["setpoint"], "--", label="setpoint")
    axes[0].plot(time, results["temperature"], label="temperature")
    axes[0].set_ylabel("Temperature (deg C)")
    axes[0].legend()
    axes[0].grid(True)

    axes[1].plot(time, results["heater_command"])
    axes[1].set_ylabel("Heater command (%)")
    axes[1].grid(True)

    axes[2].plot(time, results["a_true"], "--", label="true a")
    axes[2].plot(time, results["a_hat"], label="estimated a")
    axes[2].plot(time, results["b_true"], "--", label="true b")
    axes[2].plot(time, results["b_hat"], label="estimated b")
    axes[2].set_xlabel("Time (s)")
    axes[2].set_ylabel("Parameters")
    axes[2].legend(ncol=2)
    axes[2].grid(True)

    fig.tight_layout()
    fig.savefig("Chapter27_Lesson5_results.png", dpi=160)
    plt.close(fig)


def print_acceptance_metrics(results: dict[str, np.ndarray]) -> None:
    start = 120  # evaluate after completion of the bumpless transfer
    error = results["setpoint"][start:] - results["temperature"][start:]
    rmse = float(np.sqrt(np.mean(error**2)))
    mae = float(np.mean(np.abs(error)))
    max_error = float(np.max(np.abs(error)))
    saturation_fraction = float(np.mean(results["saturation"][start:]))
    adaptation_fraction = float(np.mean(results["adaptation_enabled"][start:]))

    print(f"RMSE after commissioning: {rmse:.3f} deg C")
    print(f"Mean absolute error:       {mae:.3f} deg C")
    print(f"Maximum absolute error:    {max_error:.3f} deg C")
    print(f"Saturation fraction:       {100.0*saturation_fraction:.2f} %")
    print(f"Estimator update fraction: {100.0*adaptation_fraction:.2f} %")
    print(
        "Final parameter estimates: "
        f"a_hat={results['a_hat'][-1]:.5f}, "
        f"b_hat={results['b_hat'][-1]:.5f}"
    )


if __name__ == "__main__":
    simulation_results = run_simulation()
    save_results(simulation_results)
    print_acceptance_metrics(simulation_results)

The C++17 implementation uses only the standard library and illustrates a deterministic, allocation-light structure suitable for translation into a real-time task after platform timing, I/O, and safety requirements are defined.

Chapter27_Lesson5.cpp

// Chapter27_Lesson5.cpp
// Industrial adaptive-control case study for a steam-heated thermal process.
// Build: g++ -std=c++17 -O2 Chapter27_Lesson5.cpp -o Chapter27_Lesson5

#include <algorithm>
#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <random>
#include <vector>

struct Sample {
    double time{};
    double reference{};
    double temperature{};
    double heater{};
    double aTrue{};
    double bTrue{};
    double aHat{};
    double bHat{};
    double innovation{};
    int saturation{};
    int adaptationEnabled{};
    double blend{};
};

double clampValue(double value, double lower, double upper) {
    return std::max(lower, std::min(upper, value));
}

double setpoint(int k) {
    if (k < 300) return 55.0;
    if (k < 600) return 65.0;
    return 58.0;
}

std::array<double, 2> trueParameters(int k) {
    if (k < 450) return {0.940, 0.035};
    return {0.965, 0.025};
}

int main() {
    constexpr int N = 900;
    constexpr double Ts = 1.0;
    constexpr double ambient = 25.0;
    constexpr double lambda = 0.995;
    constexpr double uMin = 0.0;
    constexpr double uMax = 100.0;
    constexpr double duMax = 3.0;
    const double desiredPole = std::exp(-Ts / 20.0);

    std::mt19937 generator(27);
    std::normal_distribution<double> measurementNoise(0.0, 0.08);
    std::normal_distribution<double> processNoise(0.0, 0.03);

    std::vector<double> y(N, ambient);
    std::vector<Sample> log(N);
    std::array<double, 2> theta{0.900, 0.030};
    double P[2][2]{ {500.0, 0.0}, {0.0, 500.0} };

    double previousU = 40.0;
    double integral = 20.0;
    constexpr double kp = 2.0;
    constexpr double ki = 0.04;

    for (int k = 0; k < N - 1; ++k) {
        const double reference = setpoint(k);
        const auto plant = trueParameters(k);
        const double a = plant[0];
        const double b = plant[1];
        const double measuredY = y[k] + measurementNoise(generator);
        const double error = reference - measuredY;

        const double integralCandidate = integral + ki * error;
        const double piRaw = 20.0 + kp * error + integralCandidate;
        const double piOutput = clampValue(piRaw, uMin, uMax);
        const bool safeIntegration =
            std::abs(piOutput - piRaw) < 1.0e-12 ||
            (piOutput >= uMax && error < 0.0) ||
            (piOutput <= uMin && error > 0.0);
        if (safeIntegration) integral = integralCandidate;

        const double x = measuredY - ambient;
        const double xr = reference - ambient;
        const double effectiveB = std::max(theta[1], 0.005);
        double adaptiveRaw =
            ((desiredPole - theta[0]) * x + (1.0 - desiredPole) * xr) /
            effectiveB;

        if (k >= 80 && std::abs(error) < 2.0 && adaptiveRaw > 3.0 && adaptiveRaw < 97.0) {
            adaptiveRaw += ((k / 15) % 2 == 0) ? 0.8 : -0.8;
        }

        double blend = 0.0;
        double desiredU = piOutput;
        if (k >= 80) {
            blend = std::min(1.0, (k - 80) / 40.0);
            desiredU = (1.0 - blend) * piOutput + blend * adaptiveRaw;
        }

        const double rateLimited = clampValue(desiredU, previousU - duMax, previousU + duMax);
        const double commandedU = clampValue(rateLimited, uMin, uMax);
        const bool saturated =
            std::abs(commandedU - desiredU) > 1.0e-9 ||
            commandedU < 0.5 || commandedU > 99.5;

        const double feedDisturbance = (k >= 620 && k < 645) ? -0.45 : 0.0;
        y[k + 1] = a * y[k] + b * commandedU + (1.0 - a) * ambient +
                   feedDisturbance + processNoise(generator);

        const std::array<double, 2> phi{measuredY - ambient, commandedU};
        const double target = y[k + 1] - ambient;
        const double prediction = theta[0] * phi[0] + theta[1] * phi[1];
        const double innovation = target - prediction;
        const bool estimatorValid =
            k >= 20 && commandedU > 1.0 && commandedU < 99.0 &&
            std::abs(innovation) < 5.0;

        if (estimatorValid) {
            const double Pphi0 = P[0][0] * phi[0] + P[0][1] * phi[1];
            const double Pphi1 = P[1][0] * phi[0] + P[1][1] * phi[1];
            const double denominator = lambda + phi[0] * Pphi0 + phi[1] * Pphi1;
            const double K0 = Pphi0 / denominator;
            const double K1 = Pphi1 / denominator;

            theta[0] = clampValue(theta[0] + K0 * innovation, 0.80, 0.995);
            theta[1] = clampValue(theta[1] + K1 * innovation, 0.005, 0.10);

            const double row0P0 = phi[0] * P[0][0] + phi[1] * P[1][0];
            const double row0P1 = phi[0] * P[0][1] + phi[1] * P[1][1];
            double newP00 = (P[0][0] - K0 * row0P0) / lambda;
            double newP01 = (P[0][1] - K0 * row0P1) / lambda;
            double newP10 = (P[1][0] - K1 * row0P0) / lambda;
            double newP11 = (P[1][1] - K1 * row0P1) / lambda;
            const double offDiagonal = 0.5 * (newP01 + newP10);
            P[0][0] = clampValue(newP00, 1.0e-5, 1.0e4);
            P[0][1] = offDiagonal;
            P[1][0] = offDiagonal;
            P[1][1] = clampValue(newP11, 1.0e-5, 1.0e4);
        } else {
            P[0][0] = clampValue(P[0][0] / lambda, 1.0e-5, 1.0e4);
            P[1][1] = clampValue(P[1][1] / lambda, 1.0e-5, 1.0e4);
            P[0][1] /= lambda;
            P[1][0] = P[0][1];
        }

        log[k] = {k * Ts, reference, y[k], commandedU, a, b, theta[0], theta[1],
                  innovation, saturated ? 1 : 0, estimatorValid ? 1 : 0, blend};
        previousU = commandedU;
    }

    log[N - 1] = log[N - 2];
    log[N - 1].time = (N - 1) * Ts;
    log[N - 1].temperature = y[N - 1];
    log[N - 1].reference = setpoint(N - 1);

    std::ofstream output("Chapter27_Lesson5_results_cpp.csv");
    output << "time,setpoint,temperature,heater,a_true,b_true,a_hat,b_hat,innovation,"
              "saturation,adaptation_enabled,blend\n";
    output << std::setprecision(10);
    for (const auto& sample : log) {
        output << sample.time << ',' << sample.reference << ',' << sample.temperature << ','
               << sample.heater << ',' << sample.aTrue << ',' << sample.bTrue << ','
               << sample.aHat << ',' << sample.bHat << ',' << sample.innovation << ','
               << sample.saturation << ',' << sample.adaptationEnabled << ','
               << sample.blend << '\n';
    }

    double squaredError = 0.0;
    double absoluteError = 0.0;
    double maximumError = 0.0;
    int saturationCount = 0;
    const int start = 120;
    for (int k = start; k < N; ++k) {
        const double error = log[k].reference - log[k].temperature;
        squaredError += error * error;
        absoluteError += std::abs(error);
        maximumError = std::max(maximumError, std::abs(error));
        saturationCount += log[k].saturation;
    }
    const double count = static_cast<double>(N - start);
    std::cout << std::fixed << std::setprecision(3)
              << "RMSE after commissioning: " << std::sqrt(squaredError / count) << " deg C\n"
              << "Mean absolute error:       " << absoluteError / count << " deg C\n"
              << "Maximum absolute error:    " << maximumError << " deg C\n"
              << "Saturation fraction:       " << 100.0 * saturationCount / count << " %\n"
              << std::setprecision(5)
              << "Final estimates: a_hat=" << theta[0] << ", b_hat=" << theta[1] << '\n';
    return 0;
}

The Java implementation uses standard classes only. In an industrial gateway, the algorithm would normally run inside a scheduled service with explicit monotonic timing, fieldbus quality flags, bounded queues, and a separate safety path.

Chapter27_Lesson5.java

// Chapter27_Lesson5.java
// Industrial adaptive-control case study for a steam-heated thermal process.
// Build and run: javac Chapter27_Lesson5.java && java Chapter27_Lesson5

import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.util.Locale;
import java.util.Random;

public final class Chapter27_Lesson5 {
    private static final int N = 900;
    private static final double TS = 1.0;
    private static final double AMBIENT = 25.0;
    private static final double LAMBDA = 0.995;
    private static final double U_MIN = 0.0;
    private static final double U_MAX = 100.0;
    private static final double DU_MAX = 3.0;

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

    private static double setpoint(int k) {
        if (k < 300) return 55.0;
        if (k < 600) return 65.0;
        return 58.0;
    }

    private static double[] trueParameters(int k) {
        if (k < 450) return new double[] {0.940, 0.035};
        return new double[] {0.965, 0.025};
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);
        final double desiredPole = Math.exp(-TS / 20.0);
        final Random random = new Random(27L);

        final double[] y = new double[N];
        final double[] reference = new double[N];
        final double[] heater = new double[N];
        final double[] aTrue = new double[N];
        final double[] bTrue = new double[N];
        final double[] aHat = new double[N];
        final double[] bHat = new double[N];
        final double[] innovationLog = new double[N];
        final int[] saturation = new int[N];
        final int[] adaptationEnabled = new int[N];
        final double[] blendLog = new double[N];

        y[0] = AMBIENT;
        final double[] theta = {0.900, 0.030};
        final double[][] p = { {500.0, 0.0}, {0.0, 500.0} };
        aHat[0] = theta[0];
        bHat[0] = theta[1];

        double previousU = 40.0;
        double integral = 20.0;
        final double kp = 2.0;
        final double ki = 0.04;

        for (int k = 0; k < N - 1; k++) {
            reference[k] = setpoint(k);
            final double[] plant = trueParameters(k);
            final double a = plant[0];
            final double b = plant[1];
            aTrue[k] = a;
            bTrue[k] = b;

            final double measuredY = y[k] + 0.08 * random.nextGaussian();
            final double error = reference[k] - measuredY;

            final double integralCandidate = integral + ki * error;
            final double piRaw = 20.0 + kp * error + integralCandidate;
            final double piOutput = clamp(piRaw, U_MIN, U_MAX);
            final boolean safeIntegration =
                Math.abs(piOutput - piRaw) < 1.0e-12
                || (piOutput >= U_MAX && error < 0.0)
                || (piOutput <= U_MIN && error > 0.0);
            if (safeIntegration) integral = integralCandidate;

            final double x = measuredY - AMBIENT;
            final double xr = reference[k] - AMBIENT;
            final double effectiveB = Math.max(theta[1], 0.005);
            double adaptiveRaw =
                ((desiredPole - theta[0]) * x + (1.0 - desiredPole) * xr)
                / effectiveB;

            if (k >= 80 && Math.abs(error) < 2.0 && adaptiveRaw > 3.0 && adaptiveRaw < 97.0) {
                adaptiveRaw += ((k / 15) % 2 == 0) ? 0.8 : -0.8;
            }

            double blend = 0.0;
            double desiredU = piOutput;
            if (k >= 80) {
                blend = Math.min(1.0, (k - 80) / 40.0);
                desiredU = (1.0 - blend) * piOutput + blend * adaptiveRaw;
            }
            blendLog[k] = blend;

            final double rateLimited = clamp(desiredU, previousU - DU_MAX, previousU + DU_MAX);
            final double commandedU = clamp(rateLimited, U_MIN, U_MAX);
            heater[k] = commandedU;
            saturation[k] =
                (Math.abs(commandedU - desiredU) > 1.0e-9 || commandedU < 0.5 || commandedU > 99.5)
                ? 1 : 0;

            final double feedDisturbance = (k >= 620 && k < 645) ? -0.45 : 0.0;
            y[k + 1] = a * y[k] + b * commandedU + (1.0 - a) * AMBIENT
                + feedDisturbance + 0.03 * random.nextGaussian();

            final double[] phi = {measuredY - AMBIENT, commandedU};
            final double target = y[k + 1] - AMBIENT;
            final double prediction = theta[0] * phi[0] + theta[1] * phi[1];
            final double innovation = target - prediction;
            innovationLog[k] = innovation;
            final boolean estimatorValid =
                k >= 20 && commandedU > 1.0 && commandedU < 99.0
                && Math.abs(innovation) < 5.0;

            if (estimatorValid) {
                final double pPhi0 = p[0][0] * phi[0] + p[0][1] * phi[1];
                final double pPhi1 = p[1][0] * phi[0] + p[1][1] * phi[1];
                final double denominator = LAMBDA + phi[0] * pPhi0 + phi[1] * pPhi1;
                final double k0 = pPhi0 / denominator;
                final double k1 = pPhi1 / denominator;

                theta[0] = clamp(theta[0] + k0 * innovation, 0.80, 0.995);
                theta[1] = clamp(theta[1] + k1 * innovation, 0.005, 0.10);

                final double phiTP0 = phi[0] * p[0][0] + phi[1] * p[1][0];
                final double phiTP1 = phi[0] * p[0][1] + phi[1] * p[1][1];
                final double newP00 = (p[0][0] - k0 * phiTP0) / LAMBDA;
                final double newP01 = (p[0][1] - k0 * phiTP1) / LAMBDA;
                final double newP10 = (p[1][0] - k1 * phiTP0) / LAMBDA;
                final double newP11 = (p[1][1] - k1 * phiTP1) / LAMBDA;
                final double offDiagonal = 0.5 * (newP01 + newP10);
                p[0][0] = clamp(newP00, 1.0e-5, 1.0e4);
                p[0][1] = offDiagonal;
                p[1][0] = offDiagonal;
                p[1][1] = clamp(newP11, 1.0e-5, 1.0e4);
                adaptationEnabled[k] = 1;
            } else {
                p[0][0] = clamp(p[0][0] / LAMBDA, 1.0e-5, 1.0e4);
                p[1][1] = clamp(p[1][1] / LAMBDA, 1.0e-5, 1.0e4);
                p[0][1] /= LAMBDA;
                p[1][0] = p[0][1];
            }

            aHat[k + 1] = theta[0];
            bHat[k + 1] = theta[1];
            previousU = commandedU;
        }

        reference[N - 1] = setpoint(N - 1);
        heater[N - 1] = heater[N - 2];
        final double[] finalPlant = trueParameters(N - 1);
        aTrue[N - 1] = finalPlant[0];
        bTrue[N - 1] = finalPlant[1];
        blendLog[N - 1] = blendLog[N - 2];

        try (BufferedWriter writer = new BufferedWriter(
                new FileWriter("Chapter27_Lesson5_results_java.csv"))) {
            writer.write("time,setpoint,temperature,heater,a_true,b_true,a_hat,b_hat,innovation,"
                + "saturation,adaptation_enabled,blend\n");
            for (int k = 0; k < N; k++) {
                writer.write(String.format(Locale.US,
                    "%.6f,%.6f,%.6f,%.6f,%.6f,%.6f,%.6f,%.6f,%.6f,%d,%d,%.6f%n",
                    k * TS, reference[k], y[k], heater[k], aTrue[k], bTrue[k],
                    aHat[k], bHat[k], innovationLog[k], saturation[k],
                    adaptationEnabled[k], blendLog[k]));
            }
        }

        final int start = 120;
        double sumSquared = 0.0;
        double sumAbsolute = 0.0;
        double maxAbsolute = 0.0;
        int saturationCount = 0;
        for (int k = start; k < N; k++) {
            final double error = reference[k] - y[k];
            sumSquared += error * error;
            sumAbsolute += Math.abs(error);
            maxAbsolute = Math.max(maxAbsolute, Math.abs(error));
            saturationCount += saturation[k];
        }
        final double count = N - start;
        System.out.printf(Locale.US, "RMSE after commissioning: %.3f deg C%n",
            Math.sqrt(sumSquared / count));
        System.out.printf(Locale.US, "Mean absolute error:       %.3f deg C%n",
            sumAbsolute / count);
        System.out.printf(Locale.US, "Maximum absolute error:    %.3f deg C%n", maxAbsolute);
        System.out.printf(Locale.US, "Saturation fraction:       %.2f %%%n",
            100.0 * saturationCount / count);
        System.out.printf(Locale.US, "Final estimates: a_hat=%.5f, b_hat=%.5f%n",
            theta[0], theta[1]);
    }
}

The MATLAB script reproduces the experiment, exports the trace, and documents a direct Simulink mapping. A production Simulink design should separate estimation, control, supervision, and limiting into independently testable subsystems.

Chapter27_Lesson5.m

% Chapter27_Lesson5.m
% Industrial adaptive-control case study for a steam-heated thermal process.
% The script uses base MATLAB. The final comments map the algorithm to Simulink.

clear; close all; clc;
rng(27);

N = 900;
Ts = 1.0;
ambient = 25.0;
lambda = 0.995;
uMin = 0.0;
uMax = 100.0;
duMax = 3.0;
desiredPole = exp(-Ts / 20.0);

y = zeros(N,1); y(1) = ambient;
r = zeros(N,1); u = zeros(N,1);
aTrue = zeros(N,1); bTrue = zeros(N,1);
aHat = zeros(N,1); bHat = zeros(N,1);
innovationLog = zeros(N,1);
saturation = false(N,1);
adaptationEnabled = false(N,1);
blendLog = zeros(N,1);

theta = [0.900; 0.030];
P = 500.0 * eye(2);
aHat(1) = theta(1); bHat(1) = theta(2);

kp = 2.0;
ki = 0.04;
piIntegral = 20.0;
previousU = 40.0;

for index = 1:N-1
    k = index - 1;
    if k < 300
        r(index) = 55.0;
    elseif k < 600
        r(index) = 65.0;
    else
        r(index) = 58.0;
    end

    if k < 450
        a = 0.940; b = 0.035;
    else
        a = 0.965; b = 0.025;
    end
    aTrue(index) = a; bTrue(index) = b;

    measuredY = y(index) + 0.08 * randn;
    error = r(index) - measuredY;

    integralCandidate = piIntegral + ki * error;
    piRaw = 20.0 + kp * error + integralCandidate;
    piOutput = min(uMax, max(uMin, piRaw));
    safeIntegration = abs(piOutput - piRaw) < 1e-12 || ...
        (piOutput >= uMax && error < 0.0) || ...
        (piOutput <= uMin && error > 0.0);
    if safeIntegration
        piIntegral = integralCandidate;
    end

    x = measuredY - ambient;
    xReference = r(index) - ambient;
    effectiveB = max(theta(2), 0.005);
    adaptiveRaw = ((desiredPole - theta(1)) * x + ...
        (1.0 - desiredPole) * xReference) / effectiveB;

    if k >= 80 && abs(error) < 2.0 && adaptiveRaw > 3.0 && adaptiveRaw < 97.0
        if mod(floor(k / 15), 2) == 0
            adaptiveRaw = adaptiveRaw + 0.8;
        else
            adaptiveRaw = adaptiveRaw - 0.8;
        end
    end

    if k < 80
        blend = 0.0;
        desiredU = piOutput;
    else
        blend = min(1.0, (k - 80) / 40.0);
        desiredU = (1.0 - blend) * piOutput + blend * adaptiveRaw;
    end
    blendLog(index) = blend;

    rateLimitedU = min(previousU + duMax, max(previousU - duMax, desiredU));
    commandedU = min(uMax, max(uMin, rateLimitedU));
    u(index) = commandedU;
    saturation(index) = abs(commandedU - desiredU) > 1e-9 || ...
        commandedU < 0.5 || commandedU > 99.5;

    if k >= 620 && k < 645
        feedDisturbance = -0.45;
    else
        feedDisturbance = 0.0;
    end
    y(index + 1) = a * y(index) + b * commandedU + ...
        (1.0 - a) * ambient + feedDisturbance + 0.03 * randn;

    phi = [measuredY - ambient; commandedU];
    target = y(index + 1) - ambient;
    innovation = target - theta' * phi;
    innovationLog(index) = innovation;
    estimatorValid = k >= 20 && commandedU > 1.0 && ...
        commandedU < 99.0 && abs(innovation) < 5.0;

    if estimatorValid
        gain = P * phi / (lambda + phi' * P * phi);
        theta = theta + gain * innovation;
        theta(1) = min(0.995, max(0.80, theta(1)));
        theta(2) = min(0.10, max(0.005, theta(2)));
        P = (P - gain * phi' * P) / lambda;
        P = 0.5 * (P + P');
        [vectors, values] = eig(P);
        values = diag(min(1e4, max(1e-5, diag(values))));
        P = vectors * values * vectors';
        adaptationEnabled(index) = true;
    else
        P = P / lambda;
        [vectors, values] = eig(P);
        values = diag(min(1e4, max(1e-5, diag(values))));
        P = vectors * values * vectors';
    end

    aHat(index + 1) = theta(1);
    bHat(index + 1) = theta(2);
    previousU = commandedU;
end

r(N) = r(N-1); u(N) = u(N-1);
if N-1 < 450
    aTrue(N) = 0.940; bTrue(N) = 0.035;
else
    aTrue(N) = 0.965; bTrue(N) = 0.025;
end
blendLog(N) = blendLog(N-1);
time = (0:N-1)' * Ts;

results = table(time, r, y, u, aTrue, bTrue, aHat, bHat, ...
    innovationLog, saturation, adaptationEnabled, blendLog, ...
    'VariableNames', {'time','setpoint','temperature','heater', ...
    'a_true','b_true','a_hat','b_hat','innovation','saturation', ...
    'adaptation_enabled','blend'});
writetable(results, 'Chapter27_Lesson5_results_matlab.csv');

startIndex = 121;
error = r(startIndex:end) - y(startIndex:end);
fprintf('RMSE after commissioning: %.3f deg C\n', sqrt(mean(error.^2)));
fprintf('Mean absolute error:       %.3f deg C\n', mean(abs(error)));
fprintf('Maximum absolute error:    %.3f deg C\n', max(abs(error)));
fprintf('Saturation fraction:       %.2f %%\n', ...
    100.0 * mean(saturation(startIndex:end)));
fprintf('Final estimates: a_hat=%.5f, b_hat=%.5f\n', theta(1), theta(2));

figure('Name','Chapter 27 Lesson 5');
tiledlayout(3,1);
nexttile;
plot(time, r, '--', time, y, 'LineWidth', 1.1); grid on;
ylabel('Temperature (deg C)'); legend('Setpoint','Temperature');
nexttile;
plot(time, u, 'LineWidth', 1.1); grid on;
ylabel('Heater command (%)');
nexttile;
plot(time, aTrue, '--', time, aHat, time, bTrue, '--', time, bHat, ...
    'LineWidth', 1.1); grid on;
xlabel('Time (s)'); ylabel('Parameters');
legend('true a','estimated a','true b','estimated b','Location','best');
exportgraphics(gcf, 'Chapter27_Lesson5_results_matlab.png', 'Resolution', 160);

% Simulink deployment mapping:
% 1. Put the thermal plant or digital twin in a Discrete State-Space,
%    MATLAB Function, or identified-model block.
% 2. Implement projected RLS in a MATLAB Function block with persistent
%    theta and P states, and gate its enable port with quality/saturation logic.
% 3. Implement the certainty-equivalent equation in a second MATLAB Function
%    block. Follow it with Rate Limiter and Saturation blocks.
% 4. Use a Manual Switch or Stateflow supervisory chart for PI/adaptive mode,
%    and a Ramp block or state variable for the bumpless blend.
% 5. Log y, r, u, theta, innovation, gate status, and fallback events using
%    Signal Logging or To Workspace blocks for acceptance-test reports.

The downloadable package contains both Chapter27_Lesson5.nb and a plain Chapter27_Lesson5.wl source file. The following code is the executable Wolfram Language content represented in the notebook.

Chapter27_Lesson5.nb

(* Chapter27_Lesson5.nb
   Wolfram Language implementation of the industrial adaptive-control case study.
   Evaluate the notebook cells in order. *)

ClearAll["Global`*"];
SeedRandom[27];

n = 900;
ts = 1.0;
ambient = 25.0;
lambda = 0.995;
uMin = 0.0;
uMax = 100.0;
duMax = 3.0;
desiredPole = Exp[-ts/20.0];
clip[value_, lower_, upper_] := Min[upper, Max[lower, value]];
reference[k_] := Piecewise[{ {55.0, k < 300}, {65.0, k < 600} }, 58.0];
plantParameters[k_] := If[k < 450, {0.940, 0.035}, {0.965, 0.025}];

y = ConstantArray[ambient, n];
r = ConstantArray[0.0, n];
u = ConstantArray[0.0, n];
aTrue = ConstantArray[0.0, n];
bTrue = ConstantArray[0.0, n];
aHat = ConstantArray[0.0, n];
bHat = ConstantArray[0.0, n];
innovationLog = ConstantArray[0.0, n];
saturation = ConstantArray[0, n];
adaptationEnabled = ConstantArray[0, n];
blendLog = ConstantArray[0.0, n];

theta = {0.900, 0.030};
pMatrix = 500.0 IdentityMatrix[2];
aHat[[1]] = theta[[1]];
bHat[[1]] = theta[[2]];
kp = 2.0;
ki = 0.04;
piIntegral = 20.0;
previousU = 40.0;

Do[
  k = index - 1;
  r[[index]] = reference[k];
  {a, b} = plantParameters[k];
  aTrue[[index]] = a;
  bTrue[[index]] = b;

  measuredY = y[[index]] + RandomVariate[NormalDistribution[0.0, 0.08]];
  error = r[[index]] - measuredY;

  integralCandidate = piIntegral + ki error;
  piRaw = 20.0 + kp error + integralCandidate;
  piOutput = clip[piRaw, uMin, uMax];
  safeIntegration = Abs[piOutput - piRaw] < 10^-12 ||
    (piOutput >= uMax && error < 0.0) ||
    (piOutput <= uMin && error > 0.0);
  If[safeIntegration, piIntegral = integralCandidate];

  x = measuredY - ambient;
  xReference = r[[index]] - ambient;
  effectiveB = Max[theta[[2]], 0.005];
  adaptiveRaw = ((desiredPole - theta[[1]]) x +
      (1.0 - desiredPole) xReference)/effectiveB;

  If[k >= 80 && Abs[error] < 2.0 && 3.0 < adaptiveRaw < 97.0,
    adaptiveRaw += If[EvenQ[Floor[k/15]], 0.8, -0.8]
  ];

  If[k < 80,
    blend = 0.0; desiredU = piOutput,
    blend = Min[1.0, (k - 80)/40.0];
    desiredU = (1.0 - blend) piOutput + blend adaptiveRaw
  ];
  blendLog[[index]] = blend;

  rateLimitedU = clip[desiredU, previousU - duMax, previousU + duMax];
  commandedU = clip[rateLimitedU, uMin, uMax];
  u[[index]] = commandedU;
  saturation[[index]] = Boole[
    Abs[commandedU - desiredU] > 10^-9 || commandedU < 0.5 || commandedU > 99.5];

  feedDisturbance = If[620 <= k < 645, -0.45, 0.0];
  y[[index + 1]] = a y[[index]] + b commandedU +
    (1.0 - a) ambient + feedDisturbance +
    RandomVariate[NormalDistribution[0.0, 0.03]];

  phi = {measuredY - ambient, commandedU};
  target = y[[index + 1]] - ambient;
  innovation = target - theta.phi;
  innovationLog[[index]] = innovation;
  estimatorValid = k >= 20 && 1.0 < commandedU < 99.0 && Abs[innovation] < 5.0;

  If[estimatorValid,
    gain = pMatrix.phi/(lambda + phi.pMatrix.phi);
    theta = theta + gain innovation;
    theta[[1]] = clip[theta[[1]], 0.80, 0.995];
    theta[[2]] = clip[theta[[2]], 0.005, 0.10];
    pMatrix = (pMatrix - Outer[Times, gain, phi].pMatrix)/lambda;
    pMatrix = (pMatrix + Transpose[pMatrix])/2.0;
    {eigenvalues, eigenvectors} = Eigensystem[pMatrix];
    eigenvalues = clip[#, 10^-5, 10^4] & /@ eigenvalues;
    pMatrix = Transpose[eigenvectors].DiagonalMatrix[eigenvalues].eigenvectors;
    adaptationEnabled[[index]] = 1,
    pMatrix = pMatrix/lambda
  ];

  aHat[[index + 1]] = theta[[1]];
  bHat[[index + 1]] = theta[[2]];
  previousU = commandedU;
, {index, 1, n - 1}];

r[[-1]] = reference[n - 1];
u[[-1]] = u[[-2]];
{aTrue[[-1]], bTrue[[-1]]} = plantParameters[n - 1];
blendLog[[-1]] = blendLog[[-2]];
time = Range[0, n - 1] ts;

startIndex = 121;
trackingError = r[[startIndex ;;]] - y[[startIndex ;;]];
Print["RMSE after commissioning: ", NumberForm[Sqrt[Mean[trackingError^2]], {6, 3}], " deg C"];
Print["Mean absolute error:       ", NumberForm[Mean[Abs[trackingError]], {6, 3}], " deg C"];
Print["Maximum absolute error:    ", NumberForm[Max[Abs[trackingError]], {6, 3}], " deg C"];
Print["Saturation fraction:       ", NumberForm[100.0 Mean[saturation[[startIndex ;;]]], {5, 2}], " %"];
Print["Final estimates: a_hat=", NumberForm[theta[[1]], {7, 5}],
  ", b_hat=", NumberForm[theta[[2]], {7, 5}]];

rows = Transpose[{time, r, y, u, aTrue, bTrue, aHat, bHat,
    innovationLog, saturation, adaptationEnabled, blendLog}];
Export["Chapter27_Lesson5_results_wolfram.csv",
  Prepend[rows, {"time", "setpoint", "temperature", "heater", "a_true", "b_true",
    "a_hat", "b_hat", "innovation", "saturation", "adaptation_enabled", "blend"}]];

trackingPlot = ListLinePlot[{Transpose[{time, r}], Transpose[{time, y}]},
  PlotLegends -> {"Setpoint", "Temperature"}, Frame -> True,
  FrameLabel -> {"Time (s)", "Temperature (deg C)"}, ImageSize -> Large];
inputPlot = ListLinePlot[Transpose[{time, u}], Frame -> True,
  FrameLabel -> {"Time (s)", "Heater command (%)"}, ImageSize -> Large];
parameterPlot = ListLinePlot[
  {Transpose[{time, aTrue}], Transpose[{time, aHat}],
   Transpose[{time, bTrue}], Transpose[{time, bHat}]},
  PlotLegends -> {"true a", "estimated a", "true b", "estimated b"},
  Frame -> True, FrameLabel -> {"Time (s)", "Parameters"}, ImageSize -> Large];
Export["Chapter27_Lesson5_results_wolfram.png",
  GraphicsColumn[{trackingPlot, inputPlot, parameterPlot}], ImageResolution -> 160];

The first 80 samples represent shadow identification while the PI controller retains full authority. The next 40 samples blend to the adaptive command. At sample 450, the plant becomes slower and less responsive, representing fouling. The forgetting factor permits the estimator to move toward the new regime. At samples 620–644, a colder-feed disturbance tests rejection.

The bounded probe improves information content near steady state, but it is deliberately small and gated. The estimator can still show parameter bias because closed-loop data are correlated and the first-order model is approximate. The central deployment result is therefore not “the estimates equal the physical parameters.” It is that the supervised loop maintains acceptable tracking and actuator behavior while all estimates, transitions, and interventions remain bounded, observable, and reversible.

Problem 1 (Derivation of the adaptive command): For the estimated deviation model \( x(k+1)=\widehat{a}x(k)+\widehat{b}u(k) \), derive a certainty-equivalent command that imposes \( x(k+1)=p_mx(k)+(1-p_m)x_r(k) \).

Solution: Equate the estimated plant prediction to the desired model:

\[ \widehat{a}x(k)+\widehat{b}u(k)=p_mx(k)+(1-p_m)x_r(k). \]

Solving for the command gives

\[ u(k)=\frac{(p_m-\widehat{a})x(k)+(1-p_m)x_r(k)}{\widehat{b} }. \]

The division is valid only if \( \widehat{b} \) is separated from zero. This is why the deployment uses a known positive control direction and a projected lower bound.

Problem 2 (Effect of fouling on discrete parameters): Suppose a continuous thermal model has \( a=e^{-T_s/\tau} \) and \( b=K(1-e^{-T_s/\tau}) \). Show qualitatively how \( a \) and \( b \) change when fouling increases \( \tau \) and decreases \( K \).

Solution: Since

\[ \frac{\partial a}{\partial\tau}= e^{-T_s/\tau}\frac{T_s}{\tau^2}>0, \]

increasing \( \tau \) increases \( a \) toward one. Moreover, \( b=K(1-a) \). Fouling decreases \( K \) and also decreases \( 1-a \), so both effects reduce \( b \).

Problem 3 (Ultimate tracking bound): Given \( e(k+1)=p_me(k)+w_e(k) \), with \( |p_m|<1 \) and \( |w_e(k)|\le\bar{w} \), derive a bound on \( \limsup |e(k)| \).

Solution: Repeated substitution gives

\[ |e(k)|\le |p_m|^k|e(0)|+ \bar{w}\sum_{j=0}^{k-1}|p_m|^j. \]

Using the geometric-series limit,

\[ \limsup_{k\to\infty}|e(k)|\le \frac{\bar{w} }{1-|p_m|}. \]

Problem 4 (Why adaptation is gated during saturation): Explain why RLS updates based on saturated operation can produce misleading controller parameters.

Solution: During saturation, the controller-requested input and the plant-applied input differ. If the estimator uses the requested value, its regression is factually wrong. Even if it uses the applied value, prolonged saturation can make the data poorly informative because the input has little variation. In addition, unmodeled nonlinear valve behavior is strongest near limits. Freezing or cautiously weighting updates prevents the estimator from interpreting a constraint as a change in plant physics.

Problem 5 (Forgetting-factor memory): For \( \lambda=0.995 \), estimate the number of samples after which the weight of an old datum decays to one half.

Solution: Solve \( \lambda^n=1/2 \):

\[ n=\frac{\ln(1/2)}{\ln(0.995)}\approx138.3. \]

With a one-second sample period, the half-life is about 138 seconds. This is only an interpretation of weighting memory; actual parameter-tracking speed also depends on excitation and covariance.

Problem 6 (Bumpless blend rate): The authority factor increases linearly from zero to one over 40 samples. If the PI and adaptive commands differ by 12 percentage points and are otherwise constant, what command change per sample is introduced by blending alone?

Solution: The blend increment is \( \Delta\alpha=1/40 \). Therefore

\[ |\Delta u_{blend}|=\Delta\alpha|u_a-u_{PI}|= \frac{12}{40}=0.3 \text{ percentage point per sample}. \]

The independent rate limiter must still be applied because the two commands are not constant in a real transition.

Industrial adaptive control is a supervised lifecycle, not merely an estimator connected to a controller formula. In the thermal-process case study, projected RLS tracks a changing control-oriented model, while a certainty-equivalent law targets a conservative reference pole. Safe deployment requires a validated fallback controller, shadow operation, excitation governance, parameter and covariance safeguards, bumpless transfer, actuator limits, explicit state transitions, acceptance metrics, and auditable rollback. These mechanisms convert the adaptive-control theory developed throughout the course into an implementation that can be tested, monitored, and responsibly operated.