A fixed PI or PID controller is normally tuned around one operating point. In a process plant, however, the apparent process gain, time constant, delay, and disturbance sensitivity change with throughput, composition, valve position, heat-transfer coefficient, vessel inventory, ambient conditions, and equipment aging. The same temperature controller may therefore be conservative at one production rate and oscillatory at another.

The adaptive strategy in this lesson is an indirect self-tuning regulator. A low-order dynamic model is estimated online, and the estimated model is converted into controller parameters. The adaptation layer does not replace ordinary feedback. Instead, it updates the feedback controller when the process dynamics move.

flowchart TD
  R["Setpoint r(k)"] --> C["Adaptive PI controller"]
  C --> A["Actuator and rate limits"]
  A --> P["Temperature / level / pressure process"]
  P --> Y["Measured output y(k)"]
  Y --> C
  A --> ID["Online first-order model estimator"]
  Y --> ID
  ID --> M["Estimated gain and time constant"]
  M --> C
  S["Supervisory logic: projection, \nexcitation, freeze, alarms"] --> ID
  S --> C
        

The architecture separates four functions: regulation, identification, controller redesign, and supervision. This separation is essential in industrial systems because an unconstrained estimator may produce a mathematically valid but physically impossible model, and a controller calculated from that model may demand unsafe actuator motion.

2.1 Well-Mixed Thermal Vessel

Consider a constant-volume, well-mixed vessel with volumetric flow \( F \), density \( \rho \), heat capacity \( c_p \), inlet temperature \( T_i \), vessel temperature \( T \), and manipulated heat rate \( Q \). Neglecting heat loss initially, the energy balance is

\[ \rho c_p V \dot{T} = \rho c_p F(T_i-T)+Q. \]

Dividing by \( \rho c_p F \) gives

\[ \frac{V}{F}\dot{T}+T = T_i+\frac{1}{\rho c_p F}Q. \]

Around a steady operating point, define deviation variables \( y_T=T-T_0 \) and \( u_T=Q-Q_0 \). The local model becomes

\[ \tau_T \dot{y}_T+y_T=K_Tu_T+d_T, \qquad \tau_T=\frac{V}{F_0}, \qquad K_T=\frac{1}{\rho c_p F_0}. \]

Thus, an increase in flow decreases both the thermal time constant and the steady-state heater-to-temperature gain. Fouling changes the model further when wall heat transfer is included.

2.2 Gravity-Drained Tank

Let \( h \) be liquid height, \( A(h) \) the cross-sectional area, \( q_i \) the inlet flow, and \( q_o=C_v\sqrt{2gh} \) the outlet flow. The mass balance for constant density is

\[ A(h)\dot{h}=q_i-C_v\sqrt{2gh}. \]

For a constant-area tank and an equilibrium \( (h_0,q_{i0}) \), first-order Taylor expansion of the outlet flow yields

\[ \delta\dot{h} = -\frac{C_v\sqrt{2g}}{2A\sqrt{h_0}}\delta h +\frac{1}{A}\delta q_i. \]

Therefore,

\[ \tau_h\delta\dot{h}+\delta h=K_h\delta q_i, \qquad \tau_h= \frac{2A\sqrt{h_0}}{C_v\sqrt{2g}}, \qquad K_h= \frac{2\sqrt{h_0}}{C_v\sqrt{2g}}. \]

Both parameters depend on the operating level. A controller tuned at a low level is consequently not dynamically equivalent to one tuned at a high level.

2.3 Isothermal Gas Vessel

For an isothermal ideal gas in a rigid vessel, \( PV=mRT \). With inlet and outlet mass flows,

\[ \frac{V}{RT}\dot{P}=\dot{m}_i-\dot{m}_o. \]

Suppose the local outlet-flow relation is \( \delta\dot{m}_o=k_P\delta P-k_u\delta u \), where positive \( u \) is defined in the direction that increases vessel pressure. Then

\[ \tau_P\delta\dot{P}+\delta P = K_P\delta u+d_P, \qquad \tau_P=\frac{V}{RTk_P}, \qquad K_P=\frac{k_u}{k_P}. \]

Gas compressibility, valve characteristics, upstream pressure, and downstream restrictions make \( k_P \) and \( k_u \) operating-point dependent.

The three balances motivate the local model

\[ G(s)=\frac{K}{\tau s+1}e^{-Ls}, \]

where \( K \) is the local process gain, \( \tau \) is the dominant time constant, and \( L \) is an effective delay due to transport, sensing, filtering, or computation.

With sampling period \( T_s \), a zero-order hold, and an integer delay \( d \), the disturbance-augmented sampled model is

\[ y_{k+1} = ay_k+bu_{k-d}+c+w_{k+1}, \]

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

The constant term \( c \) absorbs a slowly varying bias, unmeasured load, or imperfect deviation-variable origin. The white sequence \( w_k \) represents equation error; measurement noise may be included separately.

3.1 Exact Discretization Derivation

For \( \tau\dot{y}+y=Ku \) and constant input over one sampling interval,

\[ y((k+1)T_s) = e^{-T_s/\tau}y(kT_s) + \int_0^{T_s} \frac{K}{\tau}e^{-(T_s-\xi)/\tau}u_k\,d\xi. \]

Evaluating the integral gives

\[ \int_0^{T_s} \frac{K}{\tau}e^{-(T_s-\xi)/\tau}\,d\xi = K(1-e^{-T_s/\tau}), \]

which proves the expressions for \( a \) and \( b \). A physically stable nonintegrating process with positive control direction satisfies \( 0<a<1 \) and \( b>0 \).

Write the model in linear regression form:

\[ y_{k+1}=\boldsymbol{\phi}_k^T\boldsymbol{\theta}+w_{k+1}, \qquad \boldsymbol{\phi}_k= \begin{bmatrix} y_k & u_{k-d} & 1 \end{bmatrix}^T, \qquad \boldsymbol{\theta}= \begin{bmatrix} a & b & c \end{bmatrix}^T. \]

Exponentially weighted least squares minimizes

\[ J_k(\boldsymbol{\vartheta}) = \sum_{i=0}^k \lambda_f^{k-i} \left( y_{i+1}-\boldsymbol{\phi}_i^T\boldsymbol{\vartheta} \right)^2, \qquad 0<\lambda_f\le 1. \]

The recursive equations are

\[ \varepsilon_k = y_{k+1}- \boldsymbol{\phi}_k^T\hat{\boldsymbol{\theta}}_{k-1}, \]

\[ \mathbf{L}_k = \frac{ \mathbf{P}_{k-1}\boldsymbol{\phi}_k }{ \lambda_f+ \boldsymbol{\phi}_k^T\mathbf{P}_{k-1} \boldsymbol{\phi}_k }, \]

\[ \hat{\boldsymbol{\theta}}_k = \hat{\boldsymbol{\theta}}_{k-1} + \mathbf{L}_k\varepsilon_k, \]

\[ \mathbf{P}_k = \frac{1}{\lambda_f} \left( \mathbf{I}-\mathbf{L}_k\boldsymbol{\phi}_k^T \right) \mathbf{P}_{k-1}. \]

4.1 Positive-Definiteness of the Covariance Matrix

Applying the matrix inversion lemma to the covariance recursion gives the information-form identity

\[ \mathbf{P}_k^{-1} = \lambda_f\mathbf{P}_{k-1}^{-1} + \boldsymbol{\phi}_k\boldsymbol{\phi}_k^T. \]

Let \( \mathbf{x}\ne\mathbf{0} \). If \( \mathbf{P}_{k-1}\succ0 \) and \( \lambda_f>0 \), then

\[ \mathbf{x}^T\mathbf{P}_k^{-1}\mathbf{x} = \lambda_f \mathbf{x}^T\mathbf{P}_{k-1}^{-1}\mathbf{x} + \left( \boldsymbol{\phi}_k^T\mathbf{x} \right)^2 >0. \]

Hence \( \mathbf{P}_k^{-1}\succ0 \) and \( \mathbf{P}_k\succ0 \). In floating-point implementations, explicit symmetrization and diagonal bounds are still useful because roundoff can slowly destroy exact symmetry.

4.2 Forgetting Factor and Effective Memory

The weight of a sample \( j \) steps old is \( \lambda_f^j \). The total weight of an infinite history is

\[ \sum_{j=0}^\infty\lambda_f^j = \frac{1}{1-\lambda_f}, \]

so \( N_{\mathrm{eff}}\approx1/(1-\lambda_f) \) is a useful memory-length approximation. For \( \lambda_f=0.995 \), the effective memory is about 200 samples. Smaller values track changes faster but amplify noise and covariance growth.

Accurate tracking does not automatically imply accurate parameter estimation. Once a feedback loop reaches a nearly constant setpoint, the regressor may contain insufficient independent variation. A standard finite-window persistent-excitation condition is

\[ \sum_{i=k}^{k+N-1} \boldsymbol{\phi}_i\boldsymbol{\phi}_i^T \succeq \alpha\mathbf{I}, \qquad \alpha>0. \]

For a three-parameter model, constant \( y_k \) and constant \( u_k \) cannot span all parameter directions. Practical remedies include small bounded probing signals, naturally occurring setpoint changes, scheduled identification tests, multiple operating-point models, or adaptation only during informative transients.

Probing must be subordinate to product quality and safety. The examples use a small square-wave dither during identification and remove it before evaluating disturbance rejection. In a real plant, excitation should be approved by the supervisory layer and disabled near constraints.

Consider the identified continuous approximation

\[ \hat{G}(s)=\frac{\hat{K}}{\hat{\tau}s+1}. \]

Use the PI controller

\[ C(s)=K_c\left(1+\frac{1}{T_i s}\right). \]

Choose \( T_i=\hat{\tau} \). The controller zero then cancels the estimated first-order pole:

\[ C(s)\hat{G}(s) = K_c\frac{\hat{\tau}s+1}{\hat{\tau}s} \frac{\hat{K}}{\hat{\tau}s+1} = \frac{K_c\hat{K}}{\hat{\tau}s}. \]

The resulting complementary sensitivity is

\[ T(s) = \frac{K_c\hat{K}}{\hat{\tau}s+K_c\hat{K}}. \]

To match the desired first-order response \( T_d(s)=1/(\lambda_c s+1) \), select

\[ T_i=\hat{\tau}, \qquad K_c= \frac{\hat{\tau}}{\hat{K}\lambda_c}. \]

The design parameter \( \lambda_c \) is the desired closed-loop time constant. A larger value gives a slower but more robust response. The implementations choose \( \lambda_c=0.8\hat{\tau} \) and then impose lower and upper bounds.

6.1 Discrete Velocity Form

The digital PI controller is implemented incrementally:

\[ \Delta u_k = K_c(e_k-e_{k-1}) + K_c\frac{T_s}{T_i}e_k, \qquad u_k=u_{k-1}+\Delta u_k. \]

Incremental implementation supports direct rate limiting and avoids a separate large integral state. Conditional integration removes the integral contribution whenever it would push the actuator farther into saturation.

The estimator is projected onto a physically plausible set:

\[ a_{\min}\le\hat{a}_k\le a_{\max}<1, \qquad 0<b_{\min}\le\hat{b}_k\le b_{\max}, \qquad |\hat{c}_k|\le c_{\max}. \]

This preserves the assumed stable positive process direction and prevents division by a nearly zero \( 1-\hat{a} \). The derived controller parameters are also bounded:

\[ K_{c,\min}\le K_c\le K_{c,\max}, \qquad T_{i,\min}\le T_i\le T_{i,\max}. \]

flowchart TD
  S["Sample y(k), applied input, and setpoint"] --> V["Validate sensors and actuator status"]
  V -->|invalid| F["Freeze adaptation and \nuse safe fallback gains"]
  V -->|valid| R["RLS update with forgetting"]
  R --> P["Project model parameters \ninto physical bounds"]
  P --> I["Convert a and b \nto estimated gain \nand time constant"]
  I --> G["Compute bounded PI gains"]
  G --> E["Calculate proportional \nand integral increments"]
  E --> L["Apply input-rate limit"]
  L --> A["Apply saturation and \nconditional integration"]
  A --> X["Optional bounded \nexcitation when permitted"]
  X --> U["Send command and \nlog diagnostics"]
  U --> S
        

A production implementation should additionally freeze or reset adaptation after sensor failure, actuator stiction, mode transfer, maintenance, manual operation, or prolonged saturation. The continuously adapting controller should never be the only protection against high pressure, overtemperature, overflow, or dry running; independent interlocks and safety-instrumented functions remain necessary.

For fixed exact estimates, no delay, and exact pole-zero cancellation, the nominal closed-loop pole is

\[ s_c=-\frac{1}{\lambda_c}, \]

which is stable for \( \lambda_c>0 \). With slowly varying estimates, the controller can be interpreted as a family of frozen-time stable PI loops. A useful engineering time-scale condition is

\[ \left\|\hat{\boldsymbol{\theta}}_{k+1} -\hat{\boldsymbol{\theta}}_k\right\| \ll \left\|\hat{\boldsymbol{\theta}}_k\right\|, \]

meaning that controller parameters should change slowly compared with the closed-loop response. Rate limits on gains or parameters can enforce this condition.

This argument is not a global adaptive-stability proof. Exact cancellation is degraded by delay, higher-order dynamics, nonlinear valves, noise, and estimation error. Robust implementation therefore uses conservative \( \lambda_c \), projection, bounded adaptation, anti-windup, and fallback gains. Chapter 22 addresses robust modifications in greater depth.

9.1 Temperature Loop

In a heating or cooling loop, changes in throughput alter residence time, while fouling and utility conditions alter process gain. Temperature sensors are usually low-noise but may be strongly filtered, producing additional apparent delay. Adaptation should be slower than the sensor filter and should be frozen during batch transitions unless the model explicitly represents them.

9.2 Level Loop

Level processes may be self-regulating, integrating, or interacting. The first-order model used here applies to a self-regulating tank with gravity outflow. A vessel with nearly constant pumped outflow behaves approximately as an integrator and requires a different controller mapping. The supervisory layer should therefore classify the process structure before enabling this adaptive PI law.

9.3 Pressure Loop

Pressure dynamics are often faster than temperature dynamics and may be limited by valve and compressor rate constraints. Pressure is also a safety-critical variable. Conservative gain bounds, fast independent trips, and bumpless transfer are mandatory. Probing signals may be unacceptable in tightly constrained pressure systems; natural load variation or scheduled tests should then supply excitation.

The Python program uses only the standard library for simulation and CSV export. Matplotlib is optional for plots. It runs the same adaptive algorithm for temperature, level, and pressure profiles, introduces a plant-parameter change, and later applies a load disturbance.

Chapter27_Lesson1.py

"""
Chapter27_Lesson1.py
Adaptive temperature/level/pressure control using:
  1) recursive least squares (RLS) with exponential forgetting,
  2) online conversion of an identified first-order model to PI gains,
  3) projection, saturation, rate limiting, and conditional integration.

The simulated process is:
    y[k+1] = a[k] y[k] + (1-a[k]) K[k] u[k-d] + disturbance + noise
where a[k] = exp(-Ts/tau[k]).

Run:
    python Chapter27_Lesson1.py
Optional:
    python Chapter27_Lesson1.py --plot --csv Chapter27_Lesson1_results.csv
"""

from __future__ import annotations

import argparse
import csv
import math
import random
from collections import deque
from dataclasses import dataclass
from typing import Deque, Dict, List, Tuple


@dataclass(frozen=True)
class ProcessConfig:
    name: str
    unit: str
    setpoint_1: float
    setpoint_2: float
    gain_1: float
    gain_2: float
    tau_1: float
    tau_2: float
    delay_steps: int
    disturbance: float
    noise_std: float


CONFIGS: Tuple[ProcessConfig, ...] = (
    ProcessConfig(
        "Temperature", "degC deviation", 4.0, 7.0,
        0.85, 0.58, 90.0, 145.0, 3, -0.7, 0.035
    ),
    ProcessConfig(
        "Level", "cm deviation", 2.5, 4.5,
        0.62, 0.92, 55.0, 78.0, 2, -0.45, 0.025
    ),
    ProcessConfig(
        "Pressure", "kPa deviation", 3.0, 5.5,
        1.15, 0.78, 34.0, 52.0, 1, -0.6, 0.04
    ),
)


def clamp(value: float, low: float, high: float) -> float:
    return max(low, min(high, value))


class RLS:
    """Three-parameter RLS estimator for y[k]=a*y[k-1]+b*u_applied[k-1]+c."""

    def __init__(self, forgetting: float = 0.995) -> None:
        if not 0.90 <= forgetting <= 1.0:
            raise ValueError("forgetting must be in [0.90, 1.0]")
        self.forgetting = forgetting
        self.theta = [0.92, 0.045, 0.0]  # [a, b, c]
        self.P = [
            [800.0, 0.0, 0.0],
            [0.0, 800.0, 0.0],
            [0.0, 0.0, 200.0],
        ]

    def update(self, y: float, phi: Tuple[float, float, float]) -> float:
        p_phi = [
            sum(self.P[i][j] * phi[j] for j in range(3))
            for i in range(3)
        ]
        denominator = self.forgetting + sum(
            phi[i] * p_phi[i] for i in range(3)
        )
        denominator = max(denominator, 1.0e-12)
        gain = [v / denominator for v in p_phi]

        prediction = sum(self.theta[i] * phi[i] for i in range(3))
        innovation = y - prediction
        self.theta = [
            self.theta[i] + gain[i] * innovation for i in range(3)
        ]

        old_p = [row[:] for row in self.P]
        for i in range(3):
            for j in range(3):
                correction = gain[i] * sum(
                    phi[m] * old_p[m][j] for m in range(3)
                )
                self.P[i][j] = (old_p[i][j] - correction) / self.forgetting

        # Numerical symmetry and mild covariance bounds.
        for i in range(3):
            self.P[i][i] = clamp(self.P[i][i], 1.0e-9, 1.0e7)
            for j in range(i + 1, 3):
                sym = 0.5 * (self.P[i][j] + self.P[j][i])
                self.P[i][j] = sym
                self.P[j][i] = sym

        # Projection: stable first-order pole, known positive process direction,
        # and bounded constant disturbance/bias term.
        self.theta[0] = clamp(self.theta[0], 0.02, 0.9995)
        self.theta[1] = clamp(self.theta[1], 1.0e-4, 3.0)
        self.theta[2] = clamp(self.theta[2], -3.0, 3.0)
        return innovation

    def continuous_parameters(self, sample_time: float) -> Tuple[float, float]:
        a_hat, b_hat, _ = self.theta
        tau_hat = -sample_time / math.log(a_hat)
        process_gain_hat = b_hat / max(1.0 - a_hat, 1.0e-6)
        return (
            clamp(process_gain_hat, 0.05, 5.0),
            clamp(tau_hat, 2.0 * sample_time, 500.0),
        )


class AdaptivePI:
    """PI controller retuned from identified K and tau."""

    def __init__(
        self,
        sample_time: float,
        u_min: float = 0.0,
        u_max: float = 12.0,
        du_max: float = 0.20,
    ) -> None:
        self.sample_time = sample_time
        self.u_min = u_min
        self.u_max = u_max
        self.du_max = du_max
        self.u = 0.0
        self.previous_error = 0.0
        self.last_kc = 1.0
        self.last_ti = 20.0

    def step(
        self,
        reference: float,
        measurement: float,
        gain_hat: float,
        tau_hat: float,
        excitation: float,
    ) -> Tuple[float, float, float]:
        # Lambda tuning: lambda increases with the identified time constant.
        desired_closed_loop_time = clamp(0.80 * tau_hat, 8.0, 180.0)
        kc = tau_hat / max(gain_hat * desired_closed_loop_time, 1.0e-6)
        ti = tau_hat
        kc = clamp(kc, 0.05, 8.0)
        ti = clamp(ti, 2.0 * self.sample_time, 500.0)

        error = reference - measurement
        proportional_increment = kc * (error - self.previous_error)
        integral_increment = kc * self.sample_time * error / ti
        requested_increment = proportional_increment + integral_increment
        requested_increment = clamp(
            requested_increment, -self.du_max, self.du_max
        )

        unsaturated = self.u + requested_increment
        drives_further_high = unsaturated > self.u_max and error > 0.0
        drives_further_low = unsaturated < self.u_min and error < 0.0

        # Conditional integration: omit the integral term if it pushes farther
        # into saturation.
        if drives_further_high or drives_further_low:
            requested_increment = clamp(
                proportional_increment, -self.du_max, self.du_max
            )

        self.u = clamp(
            self.u + requested_increment + excitation,
            self.u_min,
            self.u_max,
        )
        self.previous_error = error
        self.last_kc = kc
        self.last_ti = ti
        return self.u, kc, ti


def reference_at(k: int, cfg: ProcessConfig) -> float:
    if k < 120:
        return 0.0
    if k < 650:
        return cfg.setpoint_1
    return cfg.setpoint_2


def simulate(
    cfg: ProcessConfig,
    seed: int,
    sample_time: float = 1.0,
    steps: int = 1500,
) -> Tuple[Dict[str, float], List[Dict[str, float]]]:
    rng = random.Random(seed)
    estimator = RLS(forgetting=0.995)
    controller = AdaptivePI(sample_time)
    delay_line: Deque[float] = deque(
        [0.0] * (cfg.delay_steps + 1),
        maxlen=cfg.delay_steps + 1,
    )

    y = 0.0
    y_previous_measured = 0.0
    rows: List[Dict[str, float]] = []
    iae = 0.0
    squared_error = 0.0
    max_delta_u = 0.0
    previous_u = 0.0

    for k in range(steps):
        changed = k >= 760
        process_gain = cfg.gain_2 if changed else cfg.gain_1
        tau = cfg.tau_2 if changed else cfg.tau_1
        a_true = math.exp(-sample_time / tau)

        u_applied = delay_line[0]
        load = cfg.disturbance if k >= 1080 else 0.0
        process_noise = rng.gauss(0.0, 0.25 * cfg.noise_std)
        y = (
            a_true * y
            + (1.0 - a_true) * process_gain * u_applied
            + (1.0 - a_true) * load
            + process_noise
        )
        measured = y + rng.gauss(0.0, cfg.noise_std)

        phi = (y_previous_measured, u_applied, 1.0)
        innovation = estimator.update(measured, phi)
        gain_hat, tau_hat = estimator.continuous_parameters(sample_time)

        reference = reference_at(k, cfg)
        # Small deterministic probing signal. It is disabled after 1050 s,
        # when disturbance rejection is evaluated.
        excitation = 0.025 * (1.0 if (k // 25) % 2 == 0 else -1.0)
        if k >= 1050:
            excitation = 0.0

        command, kc, ti = controller.step(
            reference, measured, gain_hat, tau_hat, excitation
        )
        delay_line.append(command)

        error = reference - measured
        if k >= 120:
            iae += abs(error) * sample_time
            squared_error += error * error
        max_delta_u = max(max_delta_u, abs(command - previous_u))
        previous_u = command
        y_previous_measured = measured

        rows.append(
            {
                "time": k * sample_time,
                "reference": reference,
                "measurement": measured,
                "command": command,
                "applied_input": u_applied,
                "a_hat": estimator.theta[0],
                "b_hat": estimator.theta[1],
                "bias_hat": estimator.theta[2],
                "gain_hat": gain_hat,
                "tau_hat": tau_hat,
                "Kc": kc,
                "Ti": ti,
                "innovation": innovation,
            }
        )

    evaluated_samples = max(steps - 120, 1)
    metrics = {
        "IAE": iae,
        "RMSE": math.sqrt(squared_error / evaluated_samples),
        "max_delta_u": max_delta_u,
        "gain_hat_final": rows[-1]["gain_hat"],
        "tau_hat_final": rows[-1]["tau_hat"],
    }
    return metrics, rows


def write_csv(path: str, all_rows: Dict[str, List[Dict[str, float]]]) -> None:
    fieldnames = ["process"] + list(next(iter(all_rows.values()))[0].keys())
    with open(path, "w", newline="", encoding="utf-8") as handle:
        writer = csv.DictWriter(handle, fieldnames=fieldnames)
        writer.writeheader()
        for process_name, rows in all_rows.items():
            for row in rows:
                writer.writerow({"process": process_name, **row})


def plot_results(all_rows: Dict[str, List[Dict[str, float]]]) -> None:
    try:
        import matplotlib.pyplot as plt
    except ImportError as exc:
        raise SystemExit(
            "Plotting requires matplotlib: python -m pip install matplotlib"
        ) from exc

    for process_name, rows in all_rows.items():
        t = [row["time"] for row in rows]
        r = [row["reference"] for row in rows]
        y = [row["measurement"] for row in rows]
        u = [row["command"] for row in rows]
        k_hat = [row["gain_hat"] for row in rows]
        tau_hat = [row["tau_hat"] for row in rows]

        plt.figure()
        plt.plot(t, r, label="reference")
        plt.plot(t, y, label="measurement")
        plt.xlabel("time (s)")
        plt.ylabel(process_name)
        plt.title(f"{process_name}: adaptive tracking")
        plt.grid(True)
        plt.legend()

        plt.figure()
        plt.plot(t, u, label="command")
        plt.xlabel("time (s)")
        plt.ylabel("manipulated input")
        plt.title(f"{process_name}: control effort")
        plt.grid(True)
        plt.legend()

        plt.figure()
        plt.plot(t, k_hat, label="estimated process gain")
        plt.plot(t, tau_hat, label="estimated time constant")
        plt.xlabel("time (s)")
        plt.title(f"{process_name}: online estimates")
        plt.grid(True)
        plt.legend()

    plt.show()


def main() -> None:
    parser = argparse.ArgumentParser()
    parser.add_argument("--plot", action="store_true")
    parser.add_argument("--csv", default="")
    args = parser.parse_args()

    all_rows: Dict[str, List[Dict[str, float]]] = {}
    for index, cfg in enumerate(CONFIGS):
        metrics, rows = simulate(cfg, seed=2027 + index)
        all_rows[cfg.name] = rows
        print(
            f"{cfg.name:11s} | "
            f"IAE={metrics['IAE']:.3f} | "
            f"RMSE={metrics['RMSE']:.4f} | "
            f"max|du|={metrics['max_delta_u']:.4f} | "
            f"K_hat(final)={metrics['gain_hat_final']:.3f} | "
            f"tau_hat(final)={metrics['tau_hat_final']:.2f}"
        )

    if args.csv:
        write_csv(args.csv, all_rows)
        print(f"Wrote {args.csv}")
    if args.plot:
        plot_results(all_rows)


if __name__ == "__main__":
    main()

The C++17 version uses fixed-size arrays for the three-parameter RLS estimator and the standard library for simulation and random noise.

Chapter27_Lesson1.cpp

// Chapter27_Lesson1.cpp
// Adaptive temperature/level/pressure control with RLS and online PI retuning.
//
// Build:
//   g++ -std=c++17 -O2 -Wall -Wextra -pedantic Chapter27_Lesson1.cpp -o Chapter27_Lesson1
// Run:
//   ./Chapter27_Lesson1

#include <algorithm>
#include <array>
#include <cmath>
#include <deque>
#include <iomanip>
#include <iostream>
#include <random>
#include <stdexcept>
#include <string>
#include <vector>

struct ProcessConfig {
    std::string name;
    std::string unit;
    double setpoint1;
    double setpoint2;
    double gain1;
    double gain2;
    double tau1;
    double tau2;
    int delaySteps;
    double disturbance;
    double noiseStd;
};

static double clampValue(double value, double low, double high) {
    return std::max(low, std::min(high, value));
}

class RLS {
public:
    explicit RLS(double forgetting = 0.995)
        : forgetting_(forgetting),
          theta_{0.92, 0.045, 0.0},
          P_{{{800.0, 0.0, 0.0},
              {0.0, 800.0, 0.0},
              {0.0, 0.0, 200.0}}} {
        if (forgetting_ < 0.90 || forgetting_ > 1.0) {
            throw std::invalid_argument("forgetting must be in [0.90, 1.0]");
        }
    }

    double update(double y, const std::array<double, 3>& phi) {
        std::array<double, 3> pPhi{};
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                pPhi[i] += P_[i][j] * phi[j];
            }
        }

        double denominator = forgetting_;
        for (int i = 0; i < 3; ++i) {
            denominator += phi[i] * pPhi[i];
        }
        denominator = std::max(denominator, 1.0e-12);

        std::array<double, 3> gain{};
        for (int i = 0; i < 3; ++i) {
            gain[i] = pPhi[i] / denominator;
        }

        double prediction = 0.0;
        for (int i = 0; i < 3; ++i) {
            prediction += theta_[i] * phi[i];
        }
        const double innovation = y - prediction;

        for (int i = 0; i < 3; ++i) {
            theta_[i] += gain[i] * innovation;
        }

        const auto oldP = P_;
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                double rowProduct = 0.0;
                for (int m = 0; m < 3; ++m) {
                    rowProduct += phi[m] * oldP[m][j];
                }
                P_[i][j] =
                    (oldP[i][j] - gain[i] * rowProduct) / forgetting_;
            }
        }

        for (int i = 0; i < 3; ++i) {
            P_[i][i] = clampValue(P_[i][i], 1.0e-9, 1.0e7);
            for (int j = i + 1; j < 3; ++j) {
                const double symmetric = 0.5 * (P_[i][j] + P_[j][i]);
                P_[i][j] = symmetric;
                P_[j][i] = symmetric;
            }
        }

        theta_[0] = clampValue(theta_[0], 0.02, 0.9995);
        theta_[1] = clampValue(theta_[1], 1.0e-4, 3.0);
        theta_[2] = clampValue(theta_[2], -3.0, 3.0);
        return innovation;
    }

    std::pair<double, double> continuousParameters(double sampleTime) const {
        const double aHat = theta_[0];
        const double bHat = theta_[1];
        const double tauHat = -sampleTime / std::log(aHat);
        const double gainHat = bHat / std::max(1.0 - aHat, 1.0e-6);
        return {
            clampValue(gainHat, 0.05, 5.0),
            clampValue(tauHat, 2.0 * sampleTime, 500.0)
        };
    }

    const std::array<double, 3>& theta() const {
        return theta_;
    }

private:
    double forgetting_;
    std::array<double, 3> theta_;
    std::array<std::array<double, 3>, 3> P_;
};

class AdaptivePI {
public:
    explicit AdaptivePI(
        double sampleTime,
        double uMin = 0.0,
        double uMax = 12.0,
        double duMax = 0.20
    )
        : sampleTime_(sampleTime),
          uMin_(uMin),
          uMax_(uMax),
          duMax_(duMax) {}

    std::array<double, 3> step(
        double reference,
        double measurement,
        double gainHat,
        double tauHat,
        double excitation
    ) {
        const double lambda =
            clampValue(0.80 * tauHat, 8.0, 180.0);
        double kc = tauHat / std::max(gainHat * lambda, 1.0e-6);
        double ti = tauHat;
        kc = clampValue(kc, 0.05, 8.0);
        ti = clampValue(ti, 2.0 * sampleTime_, 500.0);

        const double error = reference - measurement;
        const double proportionalIncrement =
            kc * (error - previousError_);
        const double integralIncrement =
            kc * sampleTime_ * error / ti;

        double requestedIncrement =
            clampValue(
                proportionalIncrement + integralIncrement,
                -duMax_,
                duMax_
            );

        const double unsaturated = u_ + requestedIncrement;
        const bool drivesFurtherHigh =
            unsaturated > uMax_ && error > 0.0;
        const bool drivesFurtherLow =
            unsaturated < uMin_ && error < 0.0;

        if (drivesFurtherHigh || drivesFurtherLow) {
            requestedIncrement =
                clampValue(proportionalIncrement, -duMax_, duMax_);
        }

        u_ = clampValue(
            u_ + requestedIncrement + excitation,
            uMin_,
            uMax_
        );
        previousError_ = error;
        return {u_, kc, ti};
    }

private:
    double sampleTime_;
    double uMin_;
    double uMax_;
    double duMax_;
    double u_ = 0.0;
    double previousError_ = 0.0;
};

struct Metrics {
    double iae = 0.0;
    double rmse = 0.0;
    double maxDeltaU = 0.0;
    double gainHatFinal = 0.0;
    double tauHatFinal = 0.0;
};

static double referenceAt(int k, const ProcessConfig& cfg) {
    if (k < 120) {
        return 0.0;
    }
    if (k < 650) {
        return cfg.setpoint1;
    }
    return cfg.setpoint2;
}

static Metrics simulate(
    const ProcessConfig& cfg,
    unsigned int seed,
    double sampleTime = 1.0,
    int steps = 1500
) {
    std::mt19937 rng(seed);
    std::normal_distribution<double> unitNormal(0.0, 1.0);

    RLS estimator(0.995);
    AdaptivePI controller(sampleTime);
    std::deque<double> delayLine(
        static_cast<std::size_t>(cfg.delaySteps + 1),
        0.0
    );

    double y = 0.0;
    double previousMeasured = 0.0;
    double previousCommand = 0.0;
    double squaredError = 0.0;
    Metrics metrics;

    for (int k = 0; k < steps; ++k) {
        const bool changed = k >= 760;
        const double processGain = changed ? cfg.gain2 : cfg.gain1;
        const double tau = changed ? cfg.tau2 : cfg.tau1;
        const double aTrue = std::exp(-sampleTime / tau);

        const double appliedInput = delayLine.front();
        const double load = k >= 1080 ? cfg.disturbance : 0.0;
        const double processNoise =
            unitNormal(rng) * 0.25 * cfg.noiseStd;
        y = aTrue * y
            + (1.0 - aTrue) * processGain * appliedInput
            + (1.0 - aTrue) * load
            + processNoise;

        const double measured =
            y + unitNormal(rng) * cfg.noiseStd;

        const std::array<double, 3> phi{
            previousMeasured,
            appliedInput,
            1.0
        };
        estimator.update(measured, phi);
        const auto [gainHat, tauHat] =
            estimator.continuousParameters(sampleTime);

        const double reference = referenceAt(k, cfg);
        double excitation =
            0.025 * (((k / 25) % 2 == 0) ? 1.0 : -1.0);
        if (k >= 1050) {
            excitation = 0.0;
        }

        const auto controllerOutput =
            controller.step(
                reference,
                measured,
                gainHat,
                tauHat,
                excitation
            );
        const double command = controllerOutput[0];

        delayLine.pop_front();
        delayLine.push_back(command);

        const double error = reference - measured;
        if (k >= 120) {
            metrics.iae += std::abs(error) * sampleTime;
            squaredError += error * error;
        }
        metrics.maxDeltaU = std::max(
            metrics.maxDeltaU,
            std::abs(command - previousCommand)
        );
        previousCommand = command;
        previousMeasured = measured;
        metrics.gainHatFinal = gainHat;
        metrics.tauHatFinal = tauHat;
    }

    const int evaluatedSamples = std::max(steps - 120, 1);
    metrics.rmse =
        std::sqrt(squaredError / static_cast<double>(evaluatedSamples));
    return metrics;
}

int main() {
    const std::vector<ProcessConfig> configs{
        {
            "Temperature", "degC deviation", 4.0, 7.0,
            0.85, 0.58, 90.0, 145.0, 3, -0.7, 0.035
        },
        {
            "Level", "cm deviation", 2.5, 4.5,
            0.62, 0.92, 55.0, 78.0, 2, -0.45, 0.025
        },
        {
            "Pressure", "kPa deviation", 3.0, 5.5,
            1.15, 0.78, 34.0, 52.0, 1, -0.6, 0.04
        }
    };

    std::cout << std::fixed << std::setprecision(4);
    for (std::size_t i = 0; i < configs.size(); ++i) {
        const Metrics metrics =
            simulate(configs[i], 2027U + static_cast<unsigned int>(i));
        std::cout
            << std::left << std::setw(11) << configs[i].name
            << " | IAE=" << std::setw(10) << metrics.iae
            << " | RMSE=" << std::setw(8) << metrics.rmse
            << " | max|du|=" << std::setw(8) << metrics.maxDeltaU
            << " | K_hat(final)=" << std::setw(7)
            << metrics.gainHatFinal
            << " | tau_hat(final)=" << metrics.tauHatFinal
            << '\n';
    }
    return 0;
}

The Java version provides the same estimator, PI redesign, projection, delay queue, and evaluation metrics using only the Java standard library.

Chapter27_Lesson1.java

// Chapter27_Lesson1.java
// Adaptive temperature/level/pressure control with RLS and online PI retuning.
//
// Build:
//   javac Chapter27_Lesson1.java
// Run:
//   java Chapter27_Lesson1

import java.util.ArrayDeque;
import java.util.Deque;
import java.util.Locale;
import java.util.Random;

public final class Chapter27_Lesson1 {
    private Chapter27_Lesson1() {}

    private static final class ProcessConfig {
        final String name;
        final String unit;
        final double setpoint1;
        final double setpoint2;
        final double gain1;
        final double gain2;
        final double tau1;
        final double tau2;
        final int delaySteps;
        final double disturbance;
        final double noiseStd;

        ProcessConfig(
            String name,
            String unit,
            double setpoint1,
            double setpoint2,
            double gain1,
            double gain2,
            double tau1,
            double tau2,
            int delaySteps,
            double disturbance,
            double noiseStd
        ) {
            this.name = name;
            this.unit = unit;
            this.setpoint1 = setpoint1;
            this.setpoint2 = setpoint2;
            this.gain1 = gain1;
            this.gain2 = gain2;
            this.tau1 = tau1;
            this.tau2 = tau2;
            this.delaySteps = delaySteps;
            this.disturbance = disturbance;
            this.noiseStd = noiseStd;
        }
    }

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

    private static final class RLS {
        private final double forgetting;
        private final double[] theta = {0.92, 0.045, 0.0};
        private final double[][] p = {
            {800.0, 0.0, 0.0},
            {0.0, 800.0, 0.0},
            {0.0, 0.0, 200.0}
        };

        RLS(double forgetting) {
            if (forgetting < 0.90 || forgetting > 1.0) {
                throw new IllegalArgumentException(
                    "forgetting must be in [0.90, 1.0]"
                );
            }
            this.forgetting = forgetting;
        }

        double update(double y, double[] phi) {
            double[] pPhi = new double[3];
            for (int i = 0; i < 3; ++i) {
                for (int j = 0; j < 3; ++j) {
                    pPhi[i] += p[i][j] * phi[j];
                }
            }

            double denominator = forgetting;
            for (int i = 0; i < 3; ++i) {
                denominator += phi[i] * pPhi[i];
            }
            denominator = Math.max(denominator, 1.0e-12);

            double[] gain = new double[3];
            for (int i = 0; i < 3; ++i) {
                gain[i] = pPhi[i] / denominator;
            }

            double prediction = 0.0;
            for (int i = 0; i < 3; ++i) {
                prediction += theta[i] * phi[i];
            }
            double innovation = y - prediction;

            for (int i = 0; i < 3; ++i) {
                theta[i] += gain[i] * innovation;
            }

            double[][] oldP = new double[3][3];
            for (int i = 0; i < 3; ++i) {
                System.arraycopy(p[i], 0, oldP[i], 0, 3);
            }

            for (int i = 0; i < 3; ++i) {
                for (int j = 0; j < 3; ++j) {
                    double rowProduct = 0.0;
                    for (int m = 0; m < 3; ++m) {
                        rowProduct += phi[m] * oldP[m][j];
                    }
                    p[i][j] =
                        (oldP[i][j] - gain[i] * rowProduct) / forgetting;
                }
            }

            for (int i = 0; i < 3; ++i) {
                p[i][i] = clamp(p[i][i], 1.0e-9, 1.0e7);
                for (int j = i + 1; j < 3; ++j) {
                    double symmetric = 0.5 * (p[i][j] + p[j][i]);
                    p[i][j] = symmetric;
                    p[j][i] = symmetric;
                }
            }

            theta[0] = clamp(theta[0], 0.02, 0.9995);
            theta[1] = clamp(theta[1], 1.0e-4, 3.0);
            theta[2] = clamp(theta[2], -3.0, 3.0);
            return innovation;
        }

        double[] continuousParameters(double sampleTime) {
            double aHat = theta[0];
            double bHat = theta[1];
            double tauHat = -sampleTime / Math.log(aHat);
            double gainHat = bHat / Math.max(1.0 - aHat, 1.0e-6);
            return new double[] {
                clamp(gainHat, 0.05, 5.0),
                clamp(tauHat, 2.0 * sampleTime, 500.0)
            };
        }
    }

    private static final class AdaptivePI {
        private final double sampleTime;
        private final double uMin;
        private final double uMax;
        private final double duMax;
        private double u = 0.0;
        private double previousError = 0.0;

        AdaptivePI(
            double sampleTime,
            double uMin,
            double uMax,
            double duMax
        ) {
            this.sampleTime = sampleTime;
            this.uMin = uMin;
            this.uMax = uMax;
            this.duMax = duMax;
        }

        double[] step(
            double reference,
            double measurement,
            double gainHat,
            double tauHat,
            double excitation
        ) {
            double lambda = clamp(0.80 * tauHat, 8.0, 180.0);
            double kc =
                tauHat / Math.max(gainHat * lambda, 1.0e-6);
            double ti = tauHat;
            kc = clamp(kc, 0.05, 8.0);
            ti = clamp(ti, 2.0 * sampleTime, 500.0);

            double error = reference - measurement;
            double proportionalIncrement =
                kc * (error - previousError);
            double integralIncrement =
                kc * sampleTime * error / ti;
            double requestedIncrement = clamp(
                proportionalIncrement + integralIncrement,
                -duMax,
                duMax
            );

            double unsaturated = u + requestedIncrement;
            boolean drivesFurtherHigh =
                unsaturated > uMax && error > 0.0;
            boolean drivesFurtherLow =
                unsaturated < uMin && error < 0.0;

            if (drivesFurtherHigh || drivesFurtherLow) {
                requestedIncrement = clamp(
                    proportionalIncrement,
                    -duMax,
                    duMax
                );
            }

            u = clamp(
                u + requestedIncrement + excitation,
                uMin,
                uMax
            );
            previousError = error;
            return new double[] {u, kc, ti};
        }
    }

    private static final class Metrics {
        double iae;
        double rmse;
        double maxDeltaU;
        double gainHatFinal;
        double tauHatFinal;
    }

    private static double referenceAt(int k, ProcessConfig cfg) {
        if (k < 120) {
            return 0.0;
        }
        if (k < 650) {
            return cfg.setpoint1;
        }
        return cfg.setpoint2;
    }

    private static Metrics simulate(
        ProcessConfig cfg,
        long seed,
        double sampleTime,
        int steps
    ) {
        Random random = new Random(seed);
        RLS estimator = new RLS(0.995);
        AdaptivePI controller =
            new AdaptivePI(sampleTime, 0.0, 12.0, 0.20);

        Deque<Double> delayLine = new ArrayDeque<>();
        for (int i = 0; i < cfg.delaySteps + 1; ++i) {
            delayLine.addLast(0.0);
        }

        double y = 0.0;
        double previousMeasured = 0.0;
        double previousCommand = 0.0;
        double squaredError = 0.0;
        Metrics metrics = new Metrics();

        for (int k = 0; k < steps; ++k) {
            boolean changed = k >= 760;
            double processGain = changed ? cfg.gain2 : cfg.gain1;
            double tau = changed ? cfg.tau2 : cfg.tau1;
            double aTrue = Math.exp(-sampleTime / tau);

            double appliedInput = delayLine.removeFirst();
            double load = k >= 1080 ? cfg.disturbance : 0.0;
            double processNoise =
                random.nextGaussian() * 0.25 * cfg.noiseStd;

            y = aTrue * y
                + (1.0 - aTrue) * processGain * appliedInput
                + (1.0 - aTrue) * load
                + processNoise;

            double measured =
                y + random.nextGaussian() * cfg.noiseStd;

            double[] phi = {
                previousMeasured,
                appliedInput,
                1.0
            };
            estimator.update(measured, phi);
            double[] estimated =
                estimator.continuousParameters(sampleTime);
            double gainHat = estimated[0];
            double tauHat = estimated[1];

            double reference = referenceAt(k, cfg);
            double excitation =
                0.025 * (((k / 25) % 2 == 0) ? 1.0 : -1.0);
            if (k >= 1050) {
                excitation = 0.0;
            }

            double[] controllerOutput = controller.step(
                reference,
                measured,
                gainHat,
                tauHat,
                excitation
            );
            double command = controllerOutput[0];
            delayLine.addLast(command);

            double error = reference - measured;
            if (k >= 120) {
                metrics.iae += Math.abs(error) * sampleTime;
                squaredError += error * error;
            }
            metrics.maxDeltaU = Math.max(
                metrics.maxDeltaU,
                Math.abs(command - previousCommand)
            );

            previousCommand = command;
            previousMeasured = measured;
            metrics.gainHatFinal = gainHat;
            metrics.tauHatFinal = tauHat;
        }

        int evaluatedSamples = Math.max(steps - 120, 1);
        metrics.rmse =
            Math.sqrt(squaredError / (double) evaluatedSamples);
        return metrics;
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.ROOT);

        ProcessConfig[] configs = {
            new ProcessConfig(
                "Temperature", "degC deviation", 4.0, 7.0,
                0.85, 0.58, 90.0, 145.0, 3, -0.7, 0.035
            ),
            new ProcessConfig(
                "Level", "cm deviation", 2.5, 4.5,
                0.62, 0.92, 55.0, 78.0, 2, -0.45, 0.025
            ),
            new ProcessConfig(
                "Pressure", "kPa deviation", 3.0, 5.5,
                1.15, 0.78, 34.0, 52.0, 1, -0.6, 0.04
            )
        };

        for (int i = 0; i < configs.length; ++i) {
            Metrics metrics =
                simulate(configs[i], 2027L + i, 1.0, 1500);
            System.out.printf(
                "%-11s | IAE=%10.4f | RMSE=%8.4f | "
                    + "max|du|=%8.4f | K_hat(final)=%7.4f | "
                    + "tau_hat(final)=%.4f%n",
                configs[i].name,
                metrics.iae,
                metrics.rmse,
                metrics.maxDeltaU,
                metrics.gainHatFinal,
                metrics.tauHatFinal
            );
        }
    }
}

The MATLAB script simulates all three process profiles, prints a summary table, and plots tracking, control effort, and parameter estimates.

Chapter27_Lesson1.m

%% Chapter27_Lesson1.m
% Adaptive temperature/level/pressure control using RLS with forgetting and
% online PI retuning from an identified first-order model.
%
% Required products:
%   MATLAB
% Optional:
%   Simulink (see Chapter27_Lesson1_Simulink.m)
%
% Run:
%   Chapter27_Lesson1

clear; clc; close all;
rng(2027, "twister");

configs(1) = struct( ...
    "name", "Temperature", ...
    "unit", "degC deviation", ...
    "setpoint1", 4.0, "setpoint2", 7.0, ...
    "gain1", 0.85, "gain2", 0.58, ...
    "tau1", 90.0, "tau2", 145.0, ...
    "delaySteps", 3, "disturbance", -0.7, "noiseStd", 0.035);

configs(2) = struct( ...
    "name", "Level", ...
    "unit", "cm deviation", ...
    "setpoint1", 2.5, "setpoint2", 4.5, ...
    "gain1", 0.62, "gain2", 0.92, ...
    "tau1", 55.0, "tau2", 78.0, ...
    "delaySteps", 2, "disturbance", -0.45, "noiseStd", 0.025);

configs(3) = struct( ...
    "name", "Pressure", ...
    "unit", "kPa deviation", ...
    "setpoint1", 3.0, "setpoint2", 5.5, ...
    "gain1", 1.15, "gain2", 0.78, ...
    "tau1", 34.0, "tau2", 52.0, ...
    "delaySteps", 1, "disturbance", -0.6, "noiseStd", 0.040);

sampleTime = 1.0;
steps = 1500;
results = cell(numel(configs), 1);
metrics = zeros(numel(configs), 5);

for i = 1:numel(configs)
    [results{i}, metrics(i, :)] = simulateProcess( ...
        configs(i), sampleTime, steps, 2027 + i - 1);
end

summaryTable = table( ...
    string({configs.name})', ...
    metrics(:, 1), metrics(:, 2), metrics(:, 3), ...
    metrics(:, 4), metrics(:, 5), ...
    'VariableNames', { ...
        'Process', 'IAE', 'RMSE', 'MaxDeltaU', ...
        'FinalGainEstimate', 'FinalTauEstimate'});
disp(summaryTable);

for i = 1:numel(configs)
    data = results{i};

    figure('Name', configs(i).name + " adaptive tracking");
    plot(data.time, data.reference, '--', 'LineWidth', 1.2);
    hold on;
    plot(data.time, data.measurement, 'LineWidth', 1.0);
    grid on;
    xlabel('Time (s)');
    ylabel(configs(i).unit);
    title(configs(i).name + ": adaptive tracking");
    legend('Reference', 'Measurement', 'Location', 'best');

    figure('Name', configs(i).name + " controller and estimates");
    tiledlayout(3, 1);

    nexttile;
    plot(data.time, data.command, 'LineWidth', 1.0);
    grid on;
    ylabel('u');
    title(configs(i).name + ": manipulated input");

    nexttile;
    plot(data.time, data.gainHat, 'LineWidth', 1.0);
    grid on;
    ylabel('K estimate');

    nexttile;
    plot(data.time, data.tauHat, 'LineWidth', 1.0);
    grid on;
    ylabel('tau estimate');
    xlabel('Time (s)');
end

function [data, metrics] = simulateProcess(cfg, sampleTime, steps, seed)
    rng(seed, "twister");

    forgetting = 0.995;
    theta = [0.92; 0.045; 0.0];
    P = diag([800.0, 800.0, 200.0]);

    uMin = 0.0;
    uMax = 12.0;
    duMax = 0.20;
    command = 0.0;
    previousCommand = 0.0;
    previousError = 0.0;
    previousMeasured = 0.0;
    y = 0.0;

    delayLine = zeros(1, cfg.delaySteps + 1);

    time = (0:steps-1)' * sampleTime;
    reference = zeros(steps, 1);
    measurement = zeros(steps, 1);
    commandLog = zeros(steps, 1);
    appliedInputLog = zeros(steps, 1);
    gainHatLog = zeros(steps, 1);
    tauHatLog = zeros(steps, 1);
    KcLog = zeros(steps, 1);
    TiLog = zeros(steps, 1);
    innovationLog = zeros(steps, 1);

    iae = 0.0;
    squaredError = 0.0;
    maxDeltaU = 0.0;

    for k = 1:steps
        changed = (k - 1) >= 760;
        if changed
            processGain = cfg.gain2;
            tau = cfg.tau2;
        else
            processGain = cfg.gain1;
            tau = cfg.tau1;
        end

        aTrue = exp(-sampleTime / tau);
        appliedInput = delayLine(1);
        if (k - 1) >= 1080
            load = cfg.disturbance;
        else
            load = 0.0;
        end

        processNoise = 0.25 * cfg.noiseStd * randn();
        y = aTrue * y ...
            + (1.0 - aTrue) * processGain * appliedInput ...
            + (1.0 - aTrue) * load ...
            + processNoise;
        measured = y + cfg.noiseStd * randn();

        phi = [previousMeasured; appliedInput; 1.0];
        pPhi = P * phi;
        denominator = max(forgetting + phi' * pPhi, 1.0e-12);
        estimatorGain = pPhi / denominator;
        prediction = theta' * phi;
        innovation = measured - prediction;
        theta = theta + estimatorGain * innovation;
        P = (P - estimatorGain * phi' * P) / forgetting;
        P = 0.5 * (P + P');

        P(1, 1) = clampValue(P(1, 1), 1.0e-9, 1.0e7);
        P(2, 2) = clampValue(P(2, 2), 1.0e-9, 1.0e7);
        P(3, 3) = clampValue(P(3, 3), 1.0e-9, 1.0e7);

        theta(1) = clampValue(theta(1), 0.02, 0.9995);
        theta(2) = clampValue(theta(2), 1.0e-4, 3.0);
        theta(3) = clampValue(theta(3), -3.0, 3.0);

        tauHat = -sampleTime / log(theta(1));
        gainHat = theta(2) / max(1.0 - theta(1), 1.0e-6);
        tauHat = clampValue(tauHat, 2.0 * sampleTime, 500.0);
        gainHat = clampValue(gainHat, 0.05, 5.0);

        r = referenceAt(k - 1, cfg);
        desiredClosedLoopTime = clampValue(0.80 * tauHat, 8.0, 180.0);
        Kc = tauHat / max(gainHat * desiredClosedLoopTime, 1.0e-6);
        Ti = tauHat;
        Kc = clampValue(Kc, 0.05, 8.0);
        Ti = clampValue(Ti, 2.0 * sampleTime, 500.0);

        error = r - measured;
        proportionalIncrement = Kc * (error - previousError);
        integralIncrement = Kc * sampleTime * error / Ti;
        requestedIncrement = clampValue( ...
            proportionalIncrement + integralIncrement, -duMax, duMax);

        unsaturated = command + requestedIncrement;
        drivesFurtherHigh = unsaturated > uMax && error > 0.0;
        drivesFurtherLow = unsaturated < uMin && error < 0.0;
        if drivesFurtherHigh || drivesFurtherLow
            requestedIncrement = clampValue( ...
                proportionalIncrement, -duMax, duMax);
        end

        excitation = 0.025;
        if mod(floor((k - 1) / 25), 2) == 1
            excitation = -excitation;
        end
        if (k - 1) >= 1050
            excitation = 0.0;
        end

        command = clampValue( ...
            command + requestedIncrement + excitation, uMin, uMax);
        delayLine = [delayLine(2:end), command]; %#ok<AGROW>

        if (k - 1) >= 120
            iae = iae + abs(error) * sampleTime;
            squaredError = squaredError + error^2;
        end
        maxDeltaU = max(maxDeltaU, abs(command - previousCommand));

        reference(k) = r;
        measurement(k) = measured;
        commandLog(k) = command;
        appliedInputLog(k) = appliedInput;
        gainHatLog(k) = gainHat;
        tauHatLog(k) = tauHat;
        KcLog(k) = Kc;
        TiLog(k) = Ti;
        innovationLog(k) = innovation;

        previousCommand = command;
        previousError = error;
        previousMeasured = measured;
    end

    evaluatedSamples = max(steps - 120, 1);
    metrics = [ ...
        iae, ...
        sqrt(squaredError / evaluatedSamples), ...
        maxDeltaU, ...
        gainHatLog(end), ...
        tauHatLog(end)];

    data = table( ...
        time, reference, measurement, commandLog, appliedInputLog, ...
        gainHatLog, tauHatLog, KcLog, TiLog, innovationLog, ...
        'VariableNames', { ...
            'time', 'reference', 'measurement', 'command', ...
            'appliedInput', 'gainHat', 'tauHat', 'Kc', 'Ti', ...
            'innovation'});
end

function r = referenceAt(k, cfg)
    if k < 120
        r = 0.0;
    elseif k < 650
        r = cfg.setpoint1;
    else
        r = cfg.setpoint2;
    end
end

function value = clampValue(value, low, high)
    value = max(low, min(high, value));
end

The following MATLAB script constructs a discrete Simulink model with a MATLAB Function block containing the RLS and adaptive PI logic. It saves the generated model as Chapter27_Lesson1_AdaptiveProcess.slx.

Chapter27_Lesson1_Simulink.m

%% Chapter27_Lesson1_Simulink.m
% Programmatically build a discrete Simulink demonstration of the adaptive
% RLS + PI controller used in Chapter 27, Lesson 1.
%
% Requirements:
%   MATLAB, Simulink, Stateflow (the MATLAB Function block infrastructure)
%
% Run:
%   Chapter27_Lesson1_Simulink
%
% The generated model is saved as Chapter27_Lesson1_AdaptiveProcess.slx.

clear; clc;

modelName = "Chapter27_Lesson1_AdaptiveProcess";
if bdIsLoaded(modelName)
    close_system(modelName, 0);
end
if isfile(modelName + ".slx")
    delete(modelName + ".slx");
end

new_system(modelName);
open_system(modelName);

set_param(modelName, ...
    "Solver", "FixedStepDiscrete", ...
    "FixedStep", "1", ...
    "StopTime", "1200", ...
    "SaveTime", "on", ...
    "TimeSaveName", "tout");

add_block("simulink/Sources/Step", modelName + "/Reference", ...
    "Position", [40 65 80 95], ...
    "Time", "100", ...
    "Before", "0", ...
    "After", "5", ...
    "SampleTime", "1");

add_block("simulink/User-Defined Functions/MATLAB Function", ...
    modelName + "/Adaptive RLS PI", ...
    "Position", [145 45 310 120]);

add_block("simulink/Discrete/Discrete Transfer Fcn", ...
    modelName + "/Nominal Process", ...
    "Position", [385 50 545 105], ...
    "Numerator", "[0 0.0105]", ...
    "Denominator", "[1 -0.9876]", ...
    "SampleTime", "1");

add_block("simulink/Sources/Step", modelName + "/Load Disturbance", ...
    "Position", [385 155 425 185], ...
    "Time", "800", ...
    "Before", "0", ...
    "After", "-0.6", ...
    "SampleTime", "1");

add_block("simulink/Math Operations/Sum", modelName + "/Output Sum", ...
    "Position", [600 64 630 101], ...
    "Inputs", "++");

add_block("simulink/Sinks/Scope", modelName + "/Tracking Scope", ...
    "Position", [760 40 810 90], ...
    "NumInputPorts", "2");

add_block("simulink/Signal Routing/Mux", modelName + "/Parameter Mux", ...
    "Position", [610 165 635 225], ...
    "Inputs", "2");

add_block("simulink/Sinks/Scope", modelName + "/Parameter Scope", ...
    "Position", [760 175 810 225]);

add_block("simulink/Sinks/To Workspace", modelName + "/y_log", ...
    "Position", [680 85 740 115], ...
    "VariableName", "y_log", ...
    "SaveFormat", "Structure With Time");

add_block("simulink/Sinks/To Workspace", modelName + "/u_log", ...
    "Position", [325 125 385 155], ...
    "VariableName", "u_log", ...
    "SaveFormat", "Structure With Time");

add_line(modelName, "Reference/1", "Adaptive RLS PI/1");
add_line(modelName, "Adaptive RLS PI/1", "Nominal Process/1");
add_line(modelName, "Nominal Process/1", "Output Sum/1");
add_line(modelName, "Load Disturbance/1", "Output Sum/2");
add_line(modelName, "Output Sum/1", "Adaptive RLS PI/2", "autorouting", "on");
add_line(modelName, "Reference/1", "Tracking Scope/1", "autorouting", "on");
add_line(modelName, "Output Sum/1", "Tracking Scope/2", "autorouting", "on");
add_line(modelName, "Adaptive RLS PI/2", "Parameter Mux/1");
add_line(modelName, "Adaptive RLS PI/3", "Parameter Mux/2");
add_line(modelName, "Parameter Mux/1", "Parameter Scope/1");
add_line(modelName, "Output Sum/1", "y_log/1", "autorouting", "on");
add_line(modelName, "Adaptive RLS PI/1", "u_log/1", "autorouting", "on");

root = sfroot;
chart = root.find("-isa", "Stateflow.EMChart", ...
    "Path", modelName + "/Adaptive RLS PI");

chart.Script = sprintf([ ...
"function [u,Kc,Ti] = fcn(r,y)\n" ...
"%%#codegen\n" ...
"persistent theta P uPrev ePrev yPrev uModelPrev\n" ...
"if isempty(theta)\n" ...
"    theta = [0.92; 0.045; 0.0];\n" ...
"    P = diag([800.0 800.0 200.0]);\n" ...
"    uPrev = 0.0;\n" ...
"    ePrev = 0.0;\n" ...
"    yPrev = 0.0;\n" ...
"    uModelPrev = 0.0;\n" ...
"end\n" ...
"Ts = 1.0;\n" ...
"lambdaF = 0.995;\n" ...
"phi = [yPrev; uModelPrev; 1.0];\n" ...
"pPhi = P*phi;\n" ...
"denom = max(lambdaF + phi'*pPhi, 1.0e-12);\n" ...
"L = pPhi/denom;\n" ...
"eps = y - theta'*phi;\n" ...
"theta = theta + L*eps;\n" ...
"P = (P - L*phi'*P)/lambdaF;\n" ...
"P = 0.5*(P + P');\n" ...
"theta(1) = min(0.9995,max(0.02,theta(1)));\n" ...
"theta(2) = min(3.0,max(1.0e-4,theta(2)));\n" ...
"theta(3) = min(3.0,max(-3.0,theta(3)));\n" ...
"tauHat = -Ts/log(theta(1));\n" ...
"gainHat = theta(2)/max(1.0-theta(1),1.0e-6);\n" ...
"tauHat = min(500.0,max(2.0,tauHat));\n" ...
"gainHat = min(5.0,max(0.05,gainHat));\n" ...
"lambdaC = min(180.0,max(8.0,0.80*tauHat));\n" ...
"Kc = tauHat/max(gainHat*lambdaC,1.0e-6);\n" ...
"Kc = min(8.0,max(0.05,Kc));\n" ...
"Ti = min(500.0,max(2.0,tauHat));\n" ...
"e = r-y;\n" ...
"duP = Kc*(e-ePrev);\n" ...
"duI = Kc*Ts*e/Ti;\n" ...
"du = min(0.20,max(-0.20,duP+duI));\n" ...
"uUnsat = uPrev + du;\n" ...
"if (uUnsat > 12.0 && e > 0.0) || (uUnsat < 0.0 && e < 0.0)\n" ...
"    du = min(0.20,max(-0.20,duP));\n" ...
"end\n" ...
"u = min(12.0,max(0.0,uPrev+du));\n" ...
"yPrev = y;\n" ...
"uModelPrev = u;\n" ...
"uPrev = u;\n" ...
"ePrev = e;\n" ...
"end\n"]);

save_system(modelName);
set_param(modelName, "SimulationCommand", "update");

disp("Created " + modelName + ".slx");
disp("Open the tracking and parameter scopes, then run the model.");

The notebook evaluates the following Wolfram Language program. It uses associations for process configurations, executes the RLS and adaptive PI loop, displays a metrics dataset, and generates tracking and parameter-estimate plots.

Chapter27_Lesson1.nb

ClearAll["Global`*"];

clamp[x_, low_, high_] := Clip[x, {low, high}];

configs = {
  <|
    "Name" -> "Temperature",
    "Unit" -> "degC deviation",
    "Setpoint1" -> 4.0,
    "Setpoint2" -> 7.0,
    "Gain1" -> 0.85,
    "Gain2" -> 0.58,
    "Tau1" -> 90.0,
    "Tau2" -> 145.0,
    "DelaySteps" -> 3,
    "Disturbance" -> -0.7,
    "NoiseStd" -> 0.035
  |>,
  <|
    "Name" -> "Level",
    "Unit" -> "cm deviation",
    "Setpoint1" -> 2.5,
    "Setpoint2" -> 4.5,
    "Gain1" -> 0.62,
    "Gain2" -> 0.92,
    "Tau1" -> 55.0,
    "Tau2" -> 78.0,
    "DelaySteps" -> 2,
    "Disturbance" -> -0.45,
    "NoiseStd" -> 0.025
  |>,
  <|
    "Name" -> "Pressure",
    "Unit" -> "kPa deviation",
    "Setpoint1" -> 3.0,
    "Setpoint2" -> 5.5,
    "Gain1" -> 1.15,
    "Gain2" -> 0.78,
    "Tau1" -> 34.0,
    "Tau2" -> 52.0,
    "DelaySteps" -> 1,
    "Disturbance" -> -0.6,
    "NoiseStd" -> 0.040
  |>
};

referenceAt[k_, cfg_] := Which[
  k < 120, 0.0,
  k < 650, cfg["Setpoint1"],
  True, cfg["Setpoint2"]
];

simulateProcess[cfg_, seed_, sampleTime_: 1.0, steps_: 1500] :=
 Module[
  {
    forgetting = 0.995,
    theta = {0.92, 0.045, 0.0},
    p = DiagonalMatrix[{800.0, 800.0, 200.0}],
    uMin = 0.0,
    uMax = 12.0,
    duMax = 0.20,
    command = 0.0,
    previousCommand = 0.0,
    previousError = 0.0,
    previousMeasured = 0.0,
    y = 0.0,
    delayLine,
    rows,
    iae = 0.0,
    squaredError = 0.0,
    maxDeltaU = 0.0,
    changed,
    processGain,
    tau,
    aTrue,
    appliedInput,
    load,
    processNoise,
    measured,
    phi,
    pPhi,
    denominator,
    estimatorGain,
    prediction,
    innovation,
    tauHat,
    gainHat,
    reference,
    desiredClosedLoopTime,
    kc,
    ti,
    error,
    proportionalIncrement,
    integralIncrement,
    requestedIncrement,
    unsaturated,
    drivesFurtherHigh,
    drivesFurtherLow,
    excitation,
    metrics
  },

  SeedRandom[seed];
  delayLine = ConstantArray[0.0, cfg["DelaySteps"] + 1];
  rows = ConstantArray[<||>, steps];

  Do[
    changed = k >= 760;
    processGain = If[changed, cfg["Gain2"], cfg["Gain1"]];
    tau = If[changed, cfg["Tau2"], cfg["Tau1"]];
    aTrue = Exp[-sampleTime/tau];

    appliedInput = First[delayLine];
    load = If[k >= 1080, cfg["Disturbance"], 0.0];
    processNoise =
      RandomVariate[NormalDistribution[0.0, 0.25 cfg["NoiseStd"]]];

    y =
      aTrue y +
      (1.0 - aTrue) processGain appliedInput +
      (1.0 - aTrue) load +
      processNoise;

    measured =
      y + RandomVariate[NormalDistribution[0.0, cfg["NoiseStd"]]];

    phi = {previousMeasured, appliedInput, 1.0};
    pPhi = p . phi;
    denominator = Max[forgetting + phi . pPhi, 1.0*^-12];
    estimatorGain = pPhi/denominator;
    prediction = theta . phi;
    innovation = measured - prediction;
    theta = theta + estimatorGain innovation;
    p = (p - Outer[Times, estimatorGain, phi] . p)/forgetting;
    p = 0.5 (p + Transpose[p]);

    p[[1, 1]] = clamp[p[[1, 1]], 1.0*^-9, 1.0*^7];
    p[[2, 2]] = clamp[p[[2, 2]], 1.0*^-9, 1.0*^7];
    p[[3, 3]] = clamp[p[[3, 3]], 1.0*^-9, 1.0*^7];

    theta[[1]] = clamp[theta[[1]], 0.02, 0.9995];
    theta[[2]] = clamp[theta[[2]], 1.0*^-4, 3.0];
    theta[[3]] = clamp[theta[[3]], -3.0, 3.0];

    tauHat = -sampleTime/Log[theta[[1]]];
    gainHat = theta[[2]]/Max[1.0 - theta[[1]], 1.0*^-6];
    tauHat = clamp[tauHat, 2.0 sampleTime, 500.0];
    gainHat = clamp[gainHat, 0.05, 5.0];

    reference = referenceAt[k, cfg];
    desiredClosedLoopTime = clamp[0.80 tauHat, 8.0, 180.0];
    kc = tauHat/Max[gainHat desiredClosedLoopTime, 1.0*^-6];
    ti = tauHat;
    kc = clamp[kc, 0.05, 8.0];
    ti = clamp[ti, 2.0 sampleTime, 500.0];

    error = reference - measured;
    proportionalIncrement = kc (error - previousError);
    integralIncrement = kc sampleTime error/ti;
    requestedIncrement =
      clamp[proportionalIncrement + integralIncrement, -duMax, duMax];

    unsaturated = command + requestedIncrement;
    drivesFurtherHigh = unsaturated > uMax && error > 0.0;
    drivesFurtherLow = unsaturated < uMin && error < 0.0;

    If[drivesFurtherHigh || drivesFurtherLow,
      requestedIncrement =
        clamp[proportionalIncrement, -duMax, duMax]
    ];

    excitation =
      If[EvenQ[Floor[k/25]], 0.025, -0.025];
    If[k >= 1050, excitation = 0.0];

    command =
      clamp[command + requestedIncrement + excitation, uMin, uMax];
    delayLine = Append[Rest[delayLine], command];

    If[k >= 120,
      iae += Abs[error] sampleTime;
      squaredError += error^2;
    ];
    maxDeltaU = Max[maxDeltaU, Abs[command - previousCommand]];

    rows[[k + 1]] = <|
      "Time" -> k sampleTime,
      "Reference" -> reference,
      "Measurement" -> measured,
      "Command" -> command,
      "AppliedInput" -> appliedInput,
      "GainHat" -> gainHat,
      "TauHat" -> tauHat,
      "Kc" -> kc,
      "Ti" -> ti,
      "Innovation" -> innovation
    |>;

    previousCommand = command;
    previousError = error;
    previousMeasured = measured;
    ,
    {k, 0, steps - 1}
  ];

  metrics = <|
    "IAE" -> iae,
    "RMSE" -> Sqrt[squaredError/Max[steps - 120, 1]],
    "MaxDeltaU" -> maxDeltaU,
    "FinalGainEstimate" -> rows[[-1]]["GainHat"],
    "FinalTauEstimate" -> rows[[-1]]["TauHat"]
  |>;

  <|"Config" -> cfg, "Rows" -> rows, "Metrics" -> metrics|>
];

results = MapIndexed[
  simulateProcess[#1, 2026 + First[#2]] &,
  configs
];

Dataset[
  Map[
    Join[
      <|"Process" -> #["Config"]["Name"]|>,
      #["Metrics"]
    ] &,
    results
  ]
]

Do[
  With[
    {
      name = result["Config"]["Name"],
      rows = result["Rows"]
    },
    Print[
      ListLinePlot[
        {
          Lookup[rows, {"Time", "Reference"}],
          Lookup[rows, {"Time", "Measurement"}]
        },
        PlotLegends -> {"Reference", "Measurement"},
        Frame -> True,
        FrameLabel -> {"Time (s)", name},
        PlotLabel -> name <> ": adaptive tracking",
        ImageSize -> Large
      ]
    ];

    Print[
      ListLinePlot[
        {
          Lookup[rows, {"Time", "GainHat"}],
          Lookup[rows, {"Time", "TauHat"}]
        },
        PlotLegends -> {"Estimated gain", "Estimated time constant"},
        Frame -> True,
        FrameLabel -> {"Time (s)", "Estimate"},
        PlotLabel -> name <> ": online parameter estimates",
        ImageSize -> Large
      ]
    ];
  ],
  {result, results}
];

Controller evaluation should separate tracking, disturbance rejection, control effort, and estimator behavior. Useful metrics are

\[ \mathrm{IAE} = \sum_{k=k_0}^{N-1}|e_k|T_s, \qquad \mathrm{RMSE} = \sqrt{ \frac{1}{N-k_0} \sum_{k=k_0}^{N-1}e_k^2 }, \]

\[ \Delta u_{\max} = \max_k|u_k-u_{k-1}|. \]

A defensible simulation or commissioning sequence is:

1. Start with conservative fixed fallback gains.
2. Verify sensor scaling, actuator direction, and saturation limits.
3. Enable estimation while retaining fixed control gains.
4. Check that projected estimates remain physically plausible.
5. Enable slow gain adaptation with strict gain-rate limits.
6. Apply approved setpoint or load tests across operating points.
7. Compare adaptive and fixed-gain baselines using identical tests.
8. Test sensor faults, saturation, manual transfer, and estimator reset.

Parameter convergence should not be judged solely by comparing estimates with nominal values. In closed loop, different low-order parameter combinations may produce similar input-output behavior. Control performance, constraint compliance, innovation statistics, and repeatability across tests are all relevant.

Problem 1 — Linearization of a Gravity Tank: Starting from \( A\dot{h}=q_i-C_v\sqrt{2gh} \), derive the small-signal time constant and gain from inlet flow to level.

Solution:

Define \( f(h,q_i)=(q_i-C_v\sqrt{2gh})/A \). At the equilibrium, \( q_{i0}=C_v\sqrt{2gh_0} \). The partial derivatives are

\[ \left.\frac{\partial f}{\partial h}\right|_0 = -\frac{C_v\sqrt{2g}}{2A\sqrt{h_0}}, \qquad \left.\frac{\partial f}{\partial q_i}\right|_0 = \frac{1}{A}. \]

Hence

\[ \delta\dot{h} = -\frac{1}{\tau_h}\delta h + \frac{K_h}{\tau_h}\delta q_i, \]

with \( \tau_h=2A\sqrt{h_0}/(C_v\sqrt{2g}) \) and \( K_h=2\sqrt{h_0}/(C_v\sqrt{2g}) \).

Problem 2 — Exact Sampled Model: Derive the zero-order-hold coefficients for \( \tau\dot{y}+y=Ku \).

Solution:

Solve the linear differential equation over one interval with constant \( u_k \):

\[ y_{k+1} = e^{-T_s/\tau}y_k + \frac{K}{\tau} \int_0^{T_s} e^{-(T_s-\xi)/\tau}d\xi\,u_k. \]

The integral equals \( \tau(1-e^{-T_s/\tau}) \). Therefore \( a=e^{-T_s/\tau} \) and \( b=K(1-a) \).

Problem 3 — RLS Covariance Positivity: Prove that the RLS covariance remains positive definite when \( \mathbf{P}_0\succ0 \) and \( \lambda_f>0 \).

Solution:

The information recursion is \( \mathbf{P}_k^{-1}= \lambda_f\mathbf{P}_{k-1}^{-1}+ \boldsymbol{\phi}_k\boldsymbol{\phi}_k^T \). For every nonzero vector \( \mathbf{x} \),

\[ \mathbf{x}^T\mathbf{P}_k^{-1}\mathbf{x} = \lambda_f \mathbf{x}^T\mathbf{P}_{k-1}^{-1}\mathbf{x} + (\boldsymbol{\phi}_k^T\mathbf{x})^2 >0. \]

Thus \( \mathbf{P}_k^{-1}\succ0 \), implying \( \mathbf{P}_k\succ0 \).

Problem 4 — PI Pole Matching: For \( G(s)=K/(\tau s+1) \), show that \( T_i=\tau \) and \( K_c=\tau/(K\lambda_c) \) produce \( T(s)=1/(\lambda_c s+1) \).

Solution:

With \( C(s)=K_c(\tau s+1)/(\tau s) \), \( C(s)G(s)=K_cK/(\tau s) \). Therefore

\[ T(s) = \frac{K_cK}{\tau s+K_cK}. \]

Substituting \( K_c=\tau/(K\lambda_c) \) gives \( T(s)=1/(\lambda_c s+1) \).

Problem 5 — Forgetting-Factor Selection: A temperature loop is sampled every 2 s. Select a constant forgetting factor whose approximate effective memory is 10 min.

Solution:

Ten minutes contains \( N_{\mathrm{eff}}=600/2=300 \) samples. Using \( N_{\mathrm{eff}}\approx1/(1-\lambda_f) \),

\[ \lambda_f \approx 1-\frac{1}{300} = 0.99667. \]

This is a starting value, not a guarantee. Noise, excitation, and the actual rate of parameter variation must be checked experimentally.

Problem 6 — Lack of Excitation: Suppose a level loop reaches a constant setpoint and both measured level and valve command become constant. Explain why all three parameters \( a,b,c \) cannot be uniquely estimated.

Solution:

The regressor becomes the constant vector \( \boldsymbol{\phi}=[\bar{y},\bar{u},1]^T \). Every outer product is proportional to \( \boldsymbol{\phi}\boldsymbol{\phi}^T \), which has rank one. Consequently, the windowed information matrix cannot be positive definite in three dimensions. Only one combination of the parameters is informed by the steady operating point. Setpoint changes, load changes, or safe probing are required to increase rank.

Problem 7 — Conditional Integration: The pressure controller is at its upper actuator limit and the pressure error is positive. Which PI term should be suppressed, and why?

Solution:

The integral increment should be suppressed because it would continue increasing the internal control demand while the actuator cannot move farther. Retaining only the bounded proportional increment prevents integral windup. Integration may resume when the error reverses or the actuator leaves saturation.

Temperature, level, and pressure loops can often be represented locally by first-order-plus-delay models, but their effective gain and time constant change with operating conditions. This lesson derived those models from energy, mass, and gas balances; formulated a forgetting-factor RLS estimator; proved covariance positivity; mapped estimated dynamics to PI gains by pole matching; and added projection, rate limits, saturation handling, excitation management, and fallback supervision. The supplied implementations demonstrate the same algorithm across five programming environments and provide a base for more advanced industrial deployment studies.