Lesson 1 distinguished direct adaptation from indirect self-tuning: a direct controller updates its own gains without first constructing an explicit estimate of the plant. The central question is therefore not “Which plant parameters should be estimated?” but rather: How should the controller be written so that its adjustable gains are mathematically meaningful and can be updated from measurable signals?

After this lesson, students should be able to:

• express a direct controller in the linear-in-adjustable-parameters form \(u=\hat{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}\);
• construct first- and second-order reference models from desired bandwidth, damping ratio, settling time, and steady-state gain;
• derive model-matching equations for state-feedback and command gains;
• identify structural conditions under which ideal controller gains exist;
• distinguish controller parameterization, plant parameterization, and reference-model parameterization;
• implement the same parameterized direct controller in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

Let \(\hat{\boldsymbol{\theta} }(t)\) denote the adjustable controller-gain vector. The controller receives a regressor \(\boldsymbol{\phi}(t)\) assembled from measured states, the command, and possibly stable filtered signals. The reference model produces the desired state \(\mathbf{x}_m(t)\). The tracking error drives a tuner that changes the controller gains directly.

flowchart TD
  R["Command r"] --> RM["Reference model"]
  RM --> XM["Desired state xm"]
  R --> REG["Regressor phi from state and command"]
  X["Measured plant state x"] --> REG
  REG --> PC["Parameterized controller: u = theta_hat^T phi"]
  TH["Adjustable controller gains theta_hat"] --> PC
  PC --> P["Unknown plant"]
  P --> X
  X --> ERR["Tracking error e = x - xm"]
  XM --> ERR
  ERR --> T["Direct gain tuner"]
  REG --> T
  T --> TH
        

No block computes an explicit estimate such as \(\hat a_0,\hat a_1,\hat b\). Nevertheless, the theoretical analysis introduces an ideal controller vector \(\boldsymbol{\theta}^{*}\). It is a mathematical matching value, not an online plant estimate.

A controller is linearly parameterized when the control input is affine or linear in its adjustable coefficients:

\[ u(t)=\hat{\boldsymbol{\theta} }^{T}(t)\boldsymbol{\phi}(t) =\sum_{i=1}^{p}\hat{\theta}_{i}(t)\phi_i(t). \]

The regressor may depend nonlinearly on measured signals, but the adjustable gains must enter linearly. For example, \(\phi_3=x_1^3\) is allowed because the gain multiplying it is still linear:

\[ u=\hat\theta_1 x_1+\hat\theta_2 x_2+\hat\theta_3 x_1^3 =\hat{\boldsymbol{\theta} }^T \begin{bmatrix}x_1&x_2&x_1^3\end{bmatrix}^{T}. \]

3.1 State-feedback plus command-feedforward gains

For a SISO state-space plant, the most important parameterization in this lesson is

\[ u(t)=\hat{\mathbf{k} }_x^{T}(t)\mathbf{x}(t) +\hat k_r(t)r(t). \]

Define the combined parameter vector and regressor as

\[ \hat{\boldsymbol{\theta} } = \begin{bmatrix} \hat{\mathbf{k} }_x\\ \hat k_r \end{bmatrix}, \qquad \boldsymbol{\phi} = \begin{bmatrix} \mathbf{x}\\ r \end{bmatrix}, \qquad u=\hat{\boldsymbol{\theta} }^T\boldsymbol{\phi}. \]

3.2 PID-like gain parameterization

A direct self-tuning PID-like law can also be written linearly in its gains. With tracking error \(e_y=r-y\) and integral state \(\dot \eta=e_y\),

\[ u=\hat k_p e_y+\hat k_i\eta+\hat k_d\dot e_y = \begin{bmatrix}\hat k_p&\hat k_i&\hat k_d\end{bmatrix} \begin{bmatrix}e_y&\eta&\dot e_y\end{bmatrix}^{T}. \]

In practice, an ideal derivative is replaced by a stable filtered derivative to avoid severe noise amplification. Detailed robustness and noise modifications appear in later chapters.

3.3 Why linear parameterization is valuable

If \(\tilde{\boldsymbol{\theta} } =\hat{\boldsymbol{\theta} }-\boldsymbol{\theta}^{*}\), then

\[ u-u^{*} = \left(\hat{\boldsymbol{\theta} } -\boldsymbol{\theta}^{*}\right)^T\boldsymbol{\phi} = \tilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}. \]

Thus the control mismatch is bilinear in parameter error and regressor. This factorization is the key algebraic property used by the Lyapunov proof in Lesson 3.

A reference model converts the command into a desired closed-loop trajectory. In state-space form,

\[ \dot{\mathbf{x} }_m=A_m\mathbf{x}_m+B_m r, \qquad y_m=C_m\mathbf{x}_m. \]

The matrix \(A_m\) must be Hurwitz. Its eigenvalues encode the desired speed and damping, while \(B_m\) and \(C_m\) determine command scaling and the model output.

4.1 First-order reference model

\[ \dot y_m=-a_m y_m+a_m r, \qquad G_m(s)=\frac{Y_m(s)}{R(s)}=\frac{a_m}{s+a_m}, \qquad a_m>0. \]

The model has unit DC gain and approximate two-percent settling time \(T_s\approx 4/a_m\).

4.2 Second-order reference model

For natural frequency \(\omega_n>0\) and damping ratio \(\zeta>0\), choose

\[ A_m= \begin{bmatrix} 0&1\\ -\omega_n^2&-2\zeta\omega_n \end{bmatrix}, \qquad B_m= \begin{bmatrix} 0\\ \omega_n^2 \end{bmatrix}, \qquad C_m= \begin{bmatrix} 1&0 \end{bmatrix}. \]

The associated transfer function is

\[ G_m(s)= \frac{\omega_n^2} {s^2+2\zeta\omega_n s+\omega_n^2}. \]

Since \(G_m(0)=1\), a constant command is reproduced without steady-state scaling in the model. For \(0<\zeta<1\), the poles are

\[ s_{1,2} = -\zeta\omega_n \pm j\omega_n\sqrt{1-\zeta^2}. \]

Approximate time-domain relations familiar from linear control remain useful:

\[ T_s\approx\frac{4}{\zeta\omega_n}, \qquad M_p\approx \exp\left( -\frac{\pi\zeta}{\sqrt{1-\zeta^2} } \right), \qquad 0<\zeta<1. \]

4.3 A reference model is a feasible target, not an arbitrary wish

Increasing \(\omega_n\) demands faster plant motion and usually larger control effort. A reference model should respect the known relative degree, actuator authority, neglected high-frequency dynamics, and measurement bandwidth. An unrealistically fast model can produce large transient gains even when the ideal matching equations are formally solvable.

Consider the unknown linear plant

\[ \dot{\mathbf{x} }=A_p\mathbf{x}+B_pu \]

and the direct controller

\[ u=\hat{\mathbf{K} }_x^{T}\mathbf{x} +\hat{\mathbf{K} }_r r. \]

For a SISO system, \(\mathbf{K}_r\) is scalar. If ideal gains \(\mathbf{K}_x^{*}\) and \(\mathbf{K}_r^{*}\) exist, substituting them into the plant gives

\[ \dot{\mathbf{x} } = \left(A_p+B_p\mathbf{K}_x^{*T}\right)\mathbf{x} +B_p\mathbf{K}_r^{*}r. \]

Exact matching with \(\dot{\mathbf{x} }_m=A_m\mathbf{x}_m+B_mr\) requires

\[ A_p+B_p\mathbf{K}_x^{*T}=A_m, \qquad B_p\mathbf{K}_r^{*}=B_m. \]

5.1 Matching condition as a range-space condition

The first equation is solvable only if every column of \(A_m-A_p\) lies in the column space of \(B_p\). Equivalently,

\[ \left(I-B_pB_p^{\dagger}\right)(A_m-A_p)=0, \qquad \left(I-B_pB_p^{\dagger}\right)B_m=0, \]

where \(B_p^{\dagger}\) is the Moore-Penrose pseudoinverse. For a nonzero SISO input vector,

\[ B_p^{\dagger} = \frac{B_p^T}{B_p^TB_p}. \]

These conditions express a structural fact: state feedback can alter the system matrix only through directions actuated by \(B_p\).

5.2 Proof of exact matching under ideal gains

Suppose the matching equations hold and the plant and model begin from the same initial state. Define \(\mathbf{e}=\mathbf{x}-\mathbf{x}_m\). Then

\[ \begin{aligned} \dot{\mathbf{e} } &= A_p\mathbf{x} +B_p\left(\mathbf{K}_x^{*T}\mathbf{x} +\mathbf{K}_r^{*}r\right) -A_m\mathbf{x}_m-B_mr\\ &= A_m\mathbf{x}-A_m\mathbf{x}_m\\ &=A_m\mathbf{e}. \end{aligned} \]

Because \(A_m\) is Hurwitz, \(\mathbf{e}(t)→\mathbf{0}\). If \(\mathbf{e}(0)=\mathbf{0}\), uniqueness of solutions gives \(\mathbf{e}(t)=\mathbf{0}\) for all time.

Consider

\[ \dot x_1=x_2, \qquad \dot x_2=-a_0x_1-a_1x_2+bu, \qquad b>0, \]

where \(a_0,a_1,b\) are unknown to the controller, but the sign of \(b\) is known. Use

\[ u=\hat k_1x_1+\hat k_2x_2+\hat k_rr. \]

The closed-loop second equation under ideal gains is

\[ \dot x_2= \left(-a_0+bk_1^{*}\right)x_1 + \left(-a_1+bk_2^{*}\right)x_2 +bk_r^{*}r. \]

Matching it to the second-order reference model from Section 4 gives

\[ -a_0+bk_1^{*}=-\omega_n^2, \qquad -a_1+bk_2^{*}=-2\zeta\omega_n, \qquad bk_r^{*}=\omega_n^2. \]

Therefore,

\[ k_1^{*}=\frac{a_0-\omega_n^2}{b}, \qquad k_2^{*}=\frac{a_1-2\zeta\omega_n}{b}, \qquad k_r^{*}=\frac{\omega_n^2}{b}. \]

These formulas are not evaluated by the direct controller. Their role is to prove that a constant gain vector capable of producing the desired model exists.

Let \(\tilde{\boldsymbol{\theta} } =\hat{\boldsymbol{\theta} }-\boldsymbol{\theta}^{*}\), with

\[ \hat{\boldsymbol{\theta} } = \begin{bmatrix}\hat k_1&\hat k_2&\hat k_r\end{bmatrix}^{T}, \qquad \boldsymbol{\phi} = \begin{bmatrix}x_1&x_2&r\end{bmatrix}^{T}. \]

Define \(\mathbf{e}=\mathbf{x}-\mathbf{x}_m\) and \(b_0=[0\;\;1]^T\), so that \(B_p=bb_0\). The actual closed-loop system can be decomposed as

\[ \dot{\mathbf{x} } = A_m\mathbf{x}+B_mr +bb_0\tilde{\boldsymbol{\theta} }^T\boldsymbol{\phi}. \]

Subtracting the reference model yields the central error model:

\[ \boxed{ \dot{\mathbf{e} } = A_m\mathbf{e} +bb_0\tilde{\boldsymbol{\theta} }^T\boldsymbol{\phi} }. \]

7.1 Derivation

\[ \begin{aligned} \dot{\mathbf{x} } &=A_p\mathbf{x} +bb_0\hat{\boldsymbol{\theta} }^T\boldsymbol{\phi}\\ &=A_p\mathbf{x} +bb_0\boldsymbol{\theta}^{*T}\boldsymbol{\phi} +bb_0\tilde{\boldsymbol{\theta} }^T\boldsymbol{\phi}\\ &=A_m\mathbf{x}+B_mr +bb_0\tilde{\boldsymbol{\theta} }^T\boldsymbol{\phi}. \end{aligned} \]

The final equality follows from the matching equations. This representation isolates all controller mismatch in one scalar term \(\tilde{\boldsymbol{\theta} }^T\boldsymbol{\phi}\).

7.2 A tuning-error direction

Let \(Q=Q^T>0\) and let \(P=P^T>0\) solve

\[ A_m^TP+PA_m=-Q. \]

The scalar \(\varepsilon_T=\mathbf{e}^TPb_0\) measures the component of the model-following error aligned with the input channel. A direct update law compatible with this parameterization is

\[ \dot{\hat{\boldsymbol{\theta} } } = -\Gamma\boldsymbol{\phi}\varepsilon_T, \qquad \Gamma=\Gamma^T>0. \]

The complete stability proof and the exact cancellation of the state-parameter cross term are developed in Lesson 3. Here, the law is introduced to make the implementation of the parameterization concrete.

The parameterized-controller concept is more general than static state feedback. Let stable linear filters generate signals \(\omega_i(t)\) from measured inputs and outputs. A dynamic controller may then be written as

\[ u(t)=\sum_{i=1}^{p}\hat\theta_i(t)\omega_i(t) =\hat{\boldsymbol{\theta} }^T(t)\boldsymbol{\omega}(t). \]

For example, with a stable monic polynomial \(\Lambda(s)\), filtered basis signals can be generated by

\[ \omega_{y,i}(s)= \frac{s^{i-1} }{\Lambda(s)}Y(s), \qquad \omega_{r,i}(s)= \frac{s^{i-1} }{\Lambda(s)}R(s). \]

Then a controller of the form

\[ U(s)= \hat{\boldsymbol{\theta} }_y^T\boldsymbol{\omega}_y(s) + \hat{\boldsymbol{\theta} }_r^T\boldsymbol{\omega}_r(s) \]

remains linear in the adjustable coefficients. This idea supports higher-order and output-feedback adaptive controllers, which are treated later in the course. The stable filters are part of the controller realization; they are not online plant models.

Regressor components can have very different physical units and magnitudes. Let \(D\) be a nonsingular diagonal scaling matrix and define

\[ \bar{\boldsymbol{\phi} }=D^{-1}\boldsymbol{\phi}, \qquad \bar{\boldsymbol{\theta} }=D^T\boldsymbol{\theta}. \]

Then the control signal is unchanged:

\[ \bar{\boldsymbol{\theta} }^T\bar{\boldsymbol{\phi} } = \boldsymbol{\theta}^T D D^{-1}\boldsymbol{\phi} = \boldsymbol{\theta}^T\boldsymbol{\phi}. \]

Thus the same physical controller can have different numerical parameter vectors. Adaptation gains must be interpreted in the chosen coordinates.

9.1 Transformation of the adaptation gain

Suppose \(\dot{\boldsymbol{\theta} } =-\Gamma\boldsymbol{\phi}\varepsilon_T\). Since \(\bar{\boldsymbol{\theta} }=D^T\boldsymbol{\theta}\),

\[ \dot{\bar{\boldsymbol{\theta} } } = -D^T\Gamma D\, \bar{\boldsymbol{\phi} }\varepsilon_T. \]

Therefore the equivalent scaled adaptation matrix is \(\bar\Gamma=D^T\Gamma D\).

9.2 Normalized update

To reduce gain changes caused only by a large regressor magnitude, a normalized form may be used:

\[ \dot{\hat{\boldsymbol{\theta} } } = -\Gamma \frac{\boldsymbol{\phi}\varepsilon_T} {m^2}, \qquad m^2=1+\boldsymbol{\phi}^T\boldsymbol{\phi}. \]

Normalization does not create information absent from the regressor. It only changes the update magnitude and numerical conditioning.

A parameterization is redundant if two different parameter vectors produce the same control signal for all admissible regressors. Suppose a nonzero vector \(\mathbf{v}\) satisfies

\[ \mathbf{v}^T\boldsymbol{\phi}(t)=0 \quad\text{for all }t. \]

Then, for any scalar \(\alpha\),

\[ \left(\boldsymbol{\theta}+\alpha\mathbf{v}\right)^T \boldsymbol{\phi} = \boldsymbol{\theta}^T\boldsymbol{\phi}. \]

Hence controller behavior cannot distinguish motion along \(\mathbf{v}\). The finite-time regressor Gramian

\[ G_\phi(t_0,T) = \int_{t_0}^{t_0+T} \boldsymbol{\phi}(\tau)\boldsymbol{\phi}^T(\tau)\,d\tau \]

reveals this issue. If \(G_\phi\) is singular, some parameter directions are unexcited. Good tracking may still be possible, but unique gain convergence cannot be expected. Persistent excitation and parameter convergence were developed in Chapter 10.

flowchart TD
  A["State the tracking objective and available measurements"] --> B["Choose a feasible reference-model order"]
  B --> C["Select desired poles, damping, bandwidth, and DC gain"]
  C --> D["Choose controller basis signals phi"]
  D --> E["Write u = theta_hat^T phi"]
  E --> F["Derive ideal model-matching equations"]
  F --> G{"Matching equations \nstructurally solvable?"}
  G -->|No| H["Revise model, controller structure, \nor measured signals"]
  G -->|Yes| I["Scale regressor components and \nchoose adaptation coordinates"]
  I --> J["Derive parameter-error contribution \nto tracking-error dynamics"]
  J --> K["Proceed to stability-law \ndesign and implementation"]
        

The matching check is conceptual in a direct controller because the true plant matrices are unknown online. It is justified through prior structural knowledge: known order, controllability pattern, relative degree, input direction, and known sign of the high-frequency gain.

The implementations use the unknown-to-controller simulation plant

\[ \dot x_1=x_2, \qquad \dot x_2=-1.6x_1-1.1x_2+1.4u. \]

Select \(\omega_n=1.8\) and \(\zeta=0.9\). Then

\[ A_m= \begin{bmatrix} 0&1\\ -3.24&-3.24 \end{bmatrix}, \qquad B_m= \begin{bmatrix} 0\\ 3.24 \end{bmatrix}. \]

The ideal matching gains, used only for post-simulation comparison, are

\[ \boldsymbol{\theta}^{*} = \begin{bmatrix} -1.1714286&-1.5285714&2.3142857 \end{bmatrix}^{T}. \]

With \(Q=I\), the Lyapunov equation gives

\[ P= \begin{bmatrix} 1.15432099&0.15432099\\ 0.15432099&0.20195092 \end{bmatrix}. \]

The demonstration uses \(\Gamma=\operatorname{diag}(6,3,6)\) and

\[ \dot{\hat{\boldsymbol{\theta} } } = -\Gamma \begin{bmatrix}x_1&x_2&r\end{bmatrix}^{T} \mathbf{e}^TPb_0. \]

The controller does not use the numerical values \(1.6,1.1,1.4\); those values appear only in the plant simulator and in the optional ideal-gain comparison.

File: Chapter13_Lesson2.py

Recommended libraries are NumPy for vector operations, SciPy for Lyapunov equations and general ODE solvers, Matplotlib for visualization, and the python-control package for linear-system analysis. The listing uses a from-scratch RK4 integrator so the adaptive dynamics are explicit.

"""
Chapter13_Lesson2.py

Direct self-tuning regulator example for a second-order plant.
The adjustable controller is linearly parameterized:

    u(t) = theta_hat(t)^T [x1(t), x2(t), r(t)]^T

The controller adapts gains directly; it does not estimate a0, a1, or b.
"""

from __future__ import annotations

import csv
from dataclasses import dataclass
from typing import Callable

import numpy as np

try:
    import matplotlib.pyplot as plt
except ImportError as exc:
    raise SystemExit(
        "Matplotlib is required for plotting. Install with: pip install matplotlib"
    ) from exc


@dataclass(frozen=True)
class Configuration:
    # Unknown-to-controller plant parameters used only by the simulator.
    a0: float = 1.6
    a1: float = 1.1
    b: float = 1.4

    # Reference-model parameters.
    natural_frequency: float = 1.8
    damping_ratio: float = 0.9

    # Fixed-step integration.
    dt: float = 0.002
    final_time: float = 20.0
    save_every: int = 10


CFG = Configuration()

# For A_m = [[0, 1], [-3.24, -3.24]], this P solves
# A_m^T P + P A_m = -I.
P = np.array(
    [
        [1.154320987654321, 0.154320987654321],
        [0.154320987654321, 0.201950922115531],
    ],
    dtype=float,
)
GAMMA = np.diag([6.0, 3.0, 6.0])
B0 = np.array([0.0, 1.0])


def reference_command(t: float) -> float:
    """Bounded command with several frequency components."""
    if t < 2.0:
        return 0.0
    if t < 10.0:
        return 1.0 + 0.35 * np.sin(0.8 * t)
    if t < 16.0:
        return -0.5 + 0.40 * np.sin(1.2 * t)
    return 0.8 + 0.25 * np.sin(0.55 * t)


def augmented_dynamics(t: float, z: np.ndarray) -> np.ndarray:
    """Dynamics of plant, reference model, and direct controller gains."""
    x = z[0:2]
    x_m = z[2:4]
    theta_hat = z[4:7]

    r = reference_command(t)
    phi = np.array([x[0], x[1], r])
    u = float(theta_hat @ phi)

    wn = CFG.natural_frequency
    zeta = CFG.damping_ratio

    x_dot = np.array(
        [
            x[1],
            -CFG.a0 * x[0] - CFG.a1 * x[1] + CFG.b * u,
        ]
    )
    x_m_dot = np.array(
        [
            x_m[1],
            -(wn**2) * x_m[0] - 2.0 * zeta * wn * x_m[1] + (wn**2) * r,
        ]
    )

    e = x - x_m
    tuning_error = float(e @ P @ B0)
    theta_dot = -(GAMMA @ phi) * tuning_error

    return np.concatenate((x_dot, x_m_dot, theta_dot))


def rk4_step(
    dynamics: Callable[[float, np.ndarray], np.ndarray],
    t: float,
    z: np.ndarray,
    dt: float,
) -> np.ndarray:
    """One classical fourth-order Runge-Kutta step."""
    k1 = dynamics(t, z)
    k2 = dynamics(t + 0.5 * dt, z + 0.5 * dt * k1)
    k3 = dynamics(t + 0.5 * dt, z + 0.5 * dt * k2)
    k4 = dynamics(t + dt, z + dt * k3)
    return z + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def simulate() -> dict[str, np.ndarray]:
    steps = int(round(CFG.final_time / CFG.dt))
    z = np.zeros(7, dtype=float)

    records: list[list[float]] = []
    for step in range(steps + 1):
        t = step * CFG.dt
        if step % CFG.save_every == 0:
            r = reference_command(t)
            phi = np.array([z[0], z[1], r])
            u = float(z[4:7] @ phi)
            records.append([t, r, *z.tolist(), u])

        if step < steps:
            z = rk4_step(augmented_dynamics, t, z, CFG.dt)
            if not np.all(np.isfinite(z)):
                raise FloatingPointError("Simulation became non-finite.")

    data = np.asarray(records)
    keys = [
        "time",
        "reference",
        "x1",
        "x2",
        "xm1",
        "xm2",
        "theta1",
        "theta2",
        "theta_r",
        "control",
    ]
    return {key: data[:, index] for index, key in enumerate(keys)}


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


def plot_results(data: dict[str, np.ndarray]) -> None:
    t = data["time"]

    plt.figure()
    plt.plot(t, data["x1"], label="Plant output x1")
    plt.plot(t, data["xm1"], "--", label="Reference-model output xm1")
    plt.plot(t, data["reference"], ":", label="Command r")
    plt.xlabel("Time (s)")
    plt.ylabel("Output")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure()
    plt.plot(t, data["theta1"], label="theta1")
    plt.plot(t, data["theta2"], label="theta2")
    plt.plot(t, data["theta_r"], label="theta_r")
    plt.xlabel("Time (s)")
    plt.ylabel("Direct controller gains")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()

    plt.figure()
    plt.plot(t, data["control"])
    plt.xlabel("Time (s)")
    plt.ylabel("Control input u")
    plt.grid(True)
    plt.tight_layout()

    plt.show()


def main() -> None:
    data = simulate()
    error = data["x1"] - data["xm1"]
    rms_error = float(np.sqrt(np.mean(error**2)))

    ideal = np.array(
        [
            (CFG.a0 - CFG.natural_frequency**2) / CFG.b,
            (CFG.a1 - 2.0 * CFG.damping_ratio * CFG.natural_frequency) / CFG.b,
            CFG.natural_frequency**2 / CFG.b,
        ]
    )

    print(f"RMS output-model error: {rms_error:.6f}")
    print("Final direct gains:", data["theta1"][-1], data["theta2"][-1], data["theta_r"][-1])
    print("Ideal matching gains (simulation knowledge only):", ideal)

    save_csv(data, "Chapter13_Lesson2_results.csv")
    plot_results(data)


if __name__ == "__main__":
    main()

File: Chapter13_Lesson2.cpp

Eigen is a common C++ library for matrix-based control implementations, while Boost.Odeint can provide general numerical integration. This version intentionally uses only the C++ standard library and implements the seven-state RK4 update directly.

// Chapter13_Lesson2.cpp
//
// Direct self-tuning regulator with controller-gain parameterization.
// Standard-library-only implementation using fixed-step RK4.

#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <stdexcept>

namespace {

constexpr double A0 = 1.6;
constexpr double A1 = 1.1;
constexpr double B = 1.4;
constexpr double WN = 1.8;
constexpr double ZETA = 0.9;
constexpr double DT = 0.002;
constexpr double FINAL_TIME = 20.0;
constexpr int SAVE_EVERY = 10;

using State = std::array<double, 7>;

double referenceCommand(double t) {
    if (t < 2.0) {
        return 0.0;
    }
    if (t < 10.0) {
        return 1.0 + 0.35 * std::sin(0.8 * t);
    }
    if (t < 16.0) {
        return -0.5 + 0.40 * std::sin(1.2 * t);
    }
    return 0.8 + 0.25 * std::sin(0.55 * t);
}

State addScaled(const State& a, const State& b, double scale) {
    State result{};
    for (std::size_t i = 0; i < result.size(); ++i) {
        result[i] = a[i] + scale * b[i];
    }
    return result;
}

State dynamics(double t, const State& z) {
    const double x1 = z[0];
    const double x2 = z[1];
    const double xm1 = z[2];
    const double xm2 = z[3];
    const double theta1 = z[4];
    const double theta2 = z[5];
    const double thetaR = z[6];

    const double r = referenceCommand(t);
    const double u = theta1 * x1 + theta2 * x2 + thetaR * r;

    const double e1 = x1 - xm1;
    const double e2 = x2 - xm2;

    // e^T P b0 with b0 = [0, 1]^T.
    const double tuningError =
        0.154320987654321 * e1 + 0.201950922115531 * e2;

    State dz{};
    dz[0] = x2;
    dz[1] = -A0 * x1 - A1 * x2 + B * u;
    dz[2] = xm2;
    dz[3] = -(WN * WN) * xm1 - 2.0 * ZETA * WN * xm2 + (WN * WN) * r;
    dz[4] = -6.0 * x1 * tuningError;
    dz[5] = -3.0 * x2 * tuningError;
    dz[6] = -6.0 * r * tuningError;
    return dz;
}

State rk4Step(double t, const State& z) {
    const State k1 = dynamics(t, z);
    const State k2 = dynamics(t + 0.5 * DT, addScaled(z, k1, 0.5 * DT));
    const State k3 = dynamics(t + 0.5 * DT, addScaled(z, k2, 0.5 * DT));
    const State k4 = dynamics(t + DT, addScaled(z, k3, DT));

    State next{};
    for (std::size_t i = 0; i < next.size(); ++i) {
        next[i] = z[i] + (DT / 6.0) *
            (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
        if (!std::isfinite(next[i])) {
            throw std::runtime_error("Simulation became non-finite.");
        }
    }
    return next;
}

}  // namespace

int main() {
    State z{};
    const int steps = static_cast<int>(std::lround(FINAL_TIME / DT));

    std::ofstream csv("Chapter13_Lesson2_results_cpp.csv");
    if (!csv) {
        std::cerr << "Could not open output CSV.\n";
        return 1;
    }

    csv << "time,reference,x1,x2,xm1,xm2,theta1,theta2,theta_r,control\n";
    csv << std::setprecision(12);

    double squaredErrorSum = 0.0;
    std::size_t savedCount = 0;

    try {
        for (int step = 0; step <= steps; ++step) {
            const double t = step * DT;
            if (step % SAVE_EVERY == 0) {
                const double r = referenceCommand(t);
                const double u = z[4] * z[0] + z[5] * z[1] + z[6] * r;
                csv << t << ',' << r << ','
                    << z[0] << ',' << z[1] << ','
                    << z[2] << ',' << z[3] << ','
                    << z[4] << ',' << z[5] << ',' << z[6] << ','
                    << u << '\n';

                const double error = z[0] - z[2];
                squaredErrorSum += error * error;
                ++savedCount;
            }

            if (step < steps) {
                z = rk4Step(t, z);
            }
        }
    } catch (const std::exception& exception) {
        std::cerr << exception.what() << '\n';
        return 1;
    }

    const double rmsError = std::sqrt(squaredErrorSum / savedCount);
    std::cout << "RMS output-model error: " << rmsError << '\n';
    std::cout << "Final direct gains: "
              << z[4] << ", " << z[5] << ", " << z[6] << '\n';
    std::cout << "Ideal matching gains (simulation knowledge only): "
              << (A0 - WN * WN) / B << ", "
              << (A1 - 2.0 * ZETA * WN) / B << ", "
              << (WN * WN) / B << '\n';

    return 0;
}

File: Chapter13_Lesson2.java

EJML and Apache Commons Math are useful Java libraries for matrix algebra, eigenvalue calculations, and numerical methods. The following implementation is dependency-free and writes the trajectory to CSV.

// Chapter13_Lesson2.java
//
// Direct self-tuning regulator with linearly parameterized controller gains.
// Uses only the Java standard library and fixed-step RK4 integration.

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

public final class Chapter13_Lesson2 {
    private static final double A0 = 1.6;
    private static final double A1 = 1.1;
    private static final double B = 1.4;
    private static final double WN = 1.8;
    private static final double ZETA = 0.9;
    private static final double DT = 0.002;
    private static final double FINAL_TIME = 20.0;
    private static final int SAVE_EVERY = 10;

    private Chapter13_Lesson2() {
        // Utility class.
    }

    private static double referenceCommand(double t) {
        if (t < 2.0) {
            return 0.0;
        }
        if (t < 10.0) {
            return 1.0 + 0.35 * Math.sin(0.8 * t);
        }
        if (t < 16.0) {
            return -0.5 + 0.40 * Math.sin(1.2 * t);
        }
        return 0.8 + 0.25 * Math.sin(0.55 * t);
    }

    private static double[] dynamics(double t, double[] z) {
        double x1 = z[0];
        double x2 = z[1];
        double xm1 = z[2];
        double xm2 = z[3];
        double theta1 = z[4];
        double theta2 = z[5];
        double thetaR = z[6];

        double r = referenceCommand(t);
        double u = theta1 * x1 + theta2 * x2 + thetaR * r;

        double e1 = x1 - xm1;
        double e2 = x2 - xm2;
        double tuningError =
            0.154320987654321 * e1 + 0.201950922115531 * e2;

        return new double[] {
            x2,
            -A0 * x1 - A1 * x2 + B * u,
            xm2,
            -(WN * WN) * xm1 - 2.0 * ZETA * WN * xm2 + (WN * WN) * r,
            -6.0 * x1 * tuningError,
            -3.0 * x2 * tuningError,
            -6.0 * r * tuningError
        };
    }

    private static double[] addScaled(double[] a, double[] b, double scale) {
        double[] result = new double[a.length];
        for (int i = 0; i < a.length; i++) {
            result[i] = a[i] + scale * b[i];
        }
        return result;
    }

    private static double[] rk4Step(double t, double[] z) {
        double[] k1 = dynamics(t, z);
        double[] k2 = dynamics(t + 0.5 * DT, addScaled(z, k1, 0.5 * DT));
        double[] k3 = dynamics(t + 0.5 * DT, addScaled(z, k2, 0.5 * DT));
        double[] k4 = dynamics(t + DT, addScaled(z, k3, DT));

        double[] next = new double[z.length];
        for (int i = 0; i < z.length; i++) {
            next[i] = z[i] + (DT / 6.0)
                * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
            if (!Double.isFinite(next[i])) {
                throw new ArithmeticException("Simulation became non-finite.");
            }
        }
        return next;
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        double[] z = new double[7];
        int steps = (int) Math.round(FINAL_TIME / DT);

        double squaredErrorSum = 0.0;
        long savedCount = 0;

        try (PrintWriter csv = new PrintWriter(
                new BufferedWriter(new FileWriter("Chapter13_Lesson2_results_java.csv")))) {

            csv.println("time,reference,x1,x2,xm1,xm2,theta1,theta2,theta_r,control");

            for (int step = 0; step <= steps; step++) {
                double t = step * DT;

                if (step % SAVE_EVERY == 0) {
                    double r = referenceCommand(t);
                    double u = z[4] * z[0] + z[5] * z[1] + z[6] * r;
                    csv.printf(
                        "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f%n",
                        t, r, z[0], z[1], z[2], z[3], z[4], z[5], z[6], u
                    );

                    double error = z[0] - z[2];
                    squaredErrorSum += error * error;
                    savedCount++;
                }

                if (step < steps) {
                    z = rk4Step(t, z);
                }
            }
        } catch (IOException exception) {
            System.err.println("File error: " + exception.getMessage());
            System.exit(1);
        } catch (ArithmeticException exception) {
            System.err.println(exception.getMessage());
            System.exit(1);
        }

        double rmsError = Math.sqrt(squaredErrorSum / savedCount);
        System.out.printf("RMS output-model error: %.6f%n", rmsError);
        System.out.printf(
            "Final direct gains: %.6f, %.6f, %.6f%n",
            z[4], z[5], z[6]
        );
        System.out.printf(
            "Ideal matching gains (simulation knowledge only): %.6f, %.6f, %.6f%n",
            (A0 - WN * WN) / B,
            (A1 - 2.0 * ZETA * WN) / B,
            (WN * WN) / B
        );
    }
}

File: Chapter13_Lesson2.m

MATLAB provides ode45 for simulation and lyap in Control System Toolbox for solving the Lyapunov equation. The script simulates the augmented plant-model-gain dynamics, exports a table, and plots tracking, gains, and control effort.

% Chapter13_Lesson2.m
% Direct self-tuning regulator with controller-gain parameterization.
% The plant parameters below are used only to simulate the unknown plant.

clear; clc; close all;

cfg.a0 = 1.6;
cfg.a1 = 1.1;
cfg.b = 1.4;
cfg.wn = 1.8;
cfg.zeta = 0.9;
cfg.Gamma = diag([6.0, 3.0, 6.0]);

% P solves A_m' P + P A_m = -I.
% With Control System Toolbox, verify using:
% A_m = [0 1; -cfg.wn^2 -2*cfg.zeta*cfg.wn];
% P_check = lyap(A_m', eye(2));
cfg.P = [ ...
    1.154320987654321, 0.154320987654321; ...
    0.154320987654321, 0.201950922115531 ...
];
cfg.b0 = [0; 1];

z0 = zeros(7, 1);
tspan = linspace(0, 20, 4001);

options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
[t, z] = ode45(@(t, z) augmentedDynamics(t, z, cfg), tspan, z0, options);

r = arrayfun(@referenceCommand, t);
u = z(:, 5).*z(:, 1) + z(:, 6).*z(:, 2) + z(:, 7).*r;
trackingError = z(:, 1) - z(:, 3);
rmsError = sqrt(mean(trackingError.^2));

idealGains = [ ...
    (cfg.a0 - cfg.wn^2) / cfg.b; ...
    (cfg.a1 - 2*cfg.zeta*cfg.wn) / cfg.b; ...
    cfg.wn^2 / cfg.b ...
];

fprintf('RMS output-model error: %.6f\n', rmsError);
fprintf('Final direct gains: %.6f, %.6f, %.6f\n', z(end, 5:7));
fprintf('Ideal matching gains (simulation knowledge only): %.6f, %.6f, %.6f\n', idealGains);

results = table( ...
    t, r, z(:,1), z(:,2), z(:,3), z(:,4), z(:,5), z(:,6), z(:,7), u, ...
    'VariableNames', {'time','reference','x1','x2','xm1','xm2', ...
                      'theta1','theta2','theta_r','control'} ...
);
writetable(results, 'Chapter13_Lesson2_results_matlab.csv');

figure;
plot(t, z(:,1), 'LineWidth', 1.2); hold on;
plot(t, z(:,3), '--', 'LineWidth', 1.2);
plot(t, r, ':', 'LineWidth', 1.0);
grid on;
xlabel('Time (s)');
ylabel('Output');
legend('Plant output x_1', 'Reference output x_{m1}', 'Command r', ...
       'Location', 'best');
title('Direct STR Tracking');

figure;
plot(t, z(:,5:7), 'LineWidth', 1.2);
grid on;
xlabel('Time (s)');
ylabel('Direct controller gains');
legend('\theta_1', '\theta_2', '\theta_r', 'Location', 'best');
title('Controller-Gain Evolution');

figure;
plot(t, u, 'LineWidth', 1.2);
grid on;
xlabel('Time (s)');
ylabel('Control input u');
title('Control Effort');

function dz = augmentedDynamics(t, z, cfg)
    x = z(1:2);
    xm = z(3:4);
    thetaHat = z(5:7);

    r = referenceCommand(t);
    phi = [x(1); x(2); r];
    u = thetaHat' * phi;

    xDot = [ ...
        x(2); ...
        -cfg.a0*x(1) - cfg.a1*x(2) + cfg.b*u ...
    ];

    xmDot = [ ...
        xm(2); ...
        -cfg.wn^2*xm(1) - 2*cfg.zeta*cfg.wn*xm(2) + cfg.wn^2*r ...
    ];

    e = x - xm;
    tuningError = e' * cfg.P * cfg.b0;
    thetaDot = -cfg.Gamma * phi * tuningError;

    dz = [xDot; xmDot; thetaDot];
end

function r = referenceCommand(t)
    if t < 2.0
        r = 0.0;
    elseif t < 10.0
        r = 1.0 + 0.35*sin(0.8*t);
    elseif t < 16.0
        r = -0.5 + 0.40*sin(1.2*t);
    else
        r = 0.8 + 0.25*sin(0.55*t);
    end
end

File: Chapter13_Lesson2_Simulink.m

The builder creates separate continuous blocks for the unknown plant, reference model, controller-gain integrator, parameterized controller, and adaptation law. The controller block receives the plant state, command, and adaptive gains; the adaptation block returns the three gain derivatives.

% Chapter13_Lesson2_Simulink.m
% Programmatically builds a continuous-time Simulink realization of the
% parameterized direct self-tuning regulator.
%
% Requirements:
%   - MATLAB
%   - Simulink
%   - Stateflow (used internally by MATLAB Function blocks)

clear; clc;

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

new_system(model);
open_system(model);

% Unknown plant used by the simulation.
Ap = [0 1; -1.6 -1.1];
bp = [0; 1.4];

% Chosen reference model.
wn = 1.8;
zeta = 0.9;
Am = [0 1; -wn^2 -2*zeta*wn];
bm = [0; wn^2];

add_block('simulink/Sources/Step', [model '/Command'], ...
    'Time', '2', 'Before', '0', 'After', '1', ...
    'Position', [40 70 70 100]);

add_block('simulink/Continuous/State-Space', [model '/Reference Model'], ...
    'A', 'Am', 'B', 'bm', 'C', 'eye(2)', 'D', 'zeros(2,1)', ...
    'X0', '[0;0]', ...
    'Position', [150 40 280 100]);

add_block('simulink/Continuous/State-Space', [model '/Plant'], ...
    'A', 'Ap', 'B', 'bp', 'C', 'eye(2)', 'D', 'zeros(2,1)', ...
    'X0', '[0;0]', ...
    'Position', [520 180 650 240]);

add_block('simulink/Continuous/Integrator', [model '/Gain Integrator'], ...
    'InitialCondition', '[0;0;0]', ...
    'Position', [360 300 400 350]);

add_block('simulink/User-Defined Functions/MATLAB Function', ...
    [model '/Parameterized Controller'], ...
    'Position', [330 150 470 230]);

add_block('simulink/User-Defined Functions/MATLAB Function', ...
    [model '/Adaptation Law'], ...
    'Position', [140 260 300 350]);

add_block('simulink/Signal Routing/Mux', [model '/Output Mux'], ...
    'Inputs', '3', ...
    'Position', [710 80 715 220]);

add_block('simulink/Sinks/Scope', [model '/Scope'], ...
    'Position', [760 105 800 145]);

add_block('simulink/Sinks/To Workspace', [model '/Plant State Data'], ...
    'VariableName', 'x_data', 'SaveFormat', 'Structure With Time', ...
    'Position', [710 250 820 280]);

add_block('simulink/Sinks/To Workspace', [model '/Model State Data'], ...
    'VariableName', 'xm_data', 'SaveFormat', 'Structure With Time', ...
    'Position', [320 40 430 70]);

add_block('simulink/Sinks/To Workspace', [model '/Gain Data'], ...
    'VariableName', 'theta_data', 'SaveFormat', 'Structure With Time', ...
    'Position', [440 310 550 340]);

% Set MATLAB Function block scripts through the Stateflow API.
root = sfroot;

controllerChart = root.find('-isa', 'Stateflow.EMChart', ...
    'Path', [model '/Parameterized Controller']);
controllerChart.Script = sprintf([ ...
    'function u = fcn(x,r,theta)\n' ...
    '%% Controller gains are the adaptive parameters.\n' ...
    'phi = [x(1); x(2); r];\n' ...
    'u = theta'' * phi;\n' ...
    'end\n' ...
]);

adaptationChart = root.find('-isa', 'Stateflow.EMChart', ...
    'Path', [model '/Adaptation Law']);
adaptationChart.Script = sprintf([ ...
    'function thetaDot = fcn(x,xm,r)\n' ...
    'P = [1.154320987654321 0.154320987654321; ...\n' ...
    '     0.154320987654321 0.201950922115531];\n' ...
    'Gamma = diag([6.0 3.0 6.0]);\n' ...
    'b0 = [0;1];\n' ...
    'e = x - xm;\n' ...
    'phi = [x(1); x(2); r];\n' ...
    'tuningError = e'' * P * b0;\n' ...
    'thetaDot = -Gamma * phi * tuningError;\n' ...
    'end\n' ...
]);

% Connections.
add_line(model, 'Command/1', 'Reference Model/1', 'autorouting', 'on');
add_line(model, 'Command/1', 'Parameterized Controller/2', 'autorouting', 'on');
add_line(model, 'Command/1', 'Adaptation Law/3', 'autorouting', 'on');

add_line(model, 'Reference Model/1', 'Adaptation Law/2', 'autorouting', 'on');
add_line(model, 'Reference Model/1', 'Model State Data/1', 'autorouting', 'on');
add_line(model, 'Reference Model/1', 'Output Mux/2', 'autorouting', 'on');

add_line(model, 'Plant/1', 'Parameterized Controller/1', 'autorouting', 'on');
add_line(model, 'Plant/1', 'Adaptation Law/1', 'autorouting', 'on');
add_line(model, 'Plant/1', 'Plant State Data/1', 'autorouting', 'on');
add_line(model, 'Plant/1', 'Output Mux/1', 'autorouting', 'on');

add_line(model, 'Adaptation Law/1', 'Gain Integrator/1', 'autorouting', 'on');
add_line(model, 'Gain Integrator/1', 'Parameterized Controller/3', 'autorouting', 'on');
add_line(model, 'Gain Integrator/1', 'Gain Data/1', 'autorouting', 'on');
add_line(model, 'Gain Integrator/1', 'Output Mux/3', 'autorouting', 'on');

add_line(model, 'Parameterized Controller/1', 'Plant/1', 'autorouting', 'on');
add_line(model, 'Output Mux/1', 'Scope/1', 'autorouting', 'on');

set_param(model, ...
    'StopTime', '20', ...
    'Solver', 'ode45', ...
    'SaveTime', 'on', ...
    'TimeSaveName', 'tout');

save_system(model);
set_param(model, 'SimulationCommand', 'update');

disp(['Built ' model '.slx']);
disp('Run the model with: sim(''Chapter13_Lesson2_DirectSTR'')');

File: Chapter13_Lesson2.nb

The notebook uses symbolic Solve for the Lyapunov equation, NDSolveValue for the augmented nonlinear differential equations, and Export for simulation data.

Notebook[{
 Cell["Chapter 13, Lesson 2: Direct STR Gain Parameterization", "Title"],
 Cell["Chapter13_Lesson2.nb", "Subtitle"],
 Cell["This notebook simulates a direct adaptive state-feedback controller whose adjustable parameters are the controller gains rather than estimates of the plant coefficients.", "Text"],

 Cell[BoxData[
  RowBox[{
   RowBox[{"ClearAll", "[", "\"Global`*\"", "]"}], ";"}]], "Input"],

 Cell[BoxData[
  ToBoxes[
   a0 = 1.6;
   a1 = 1.1;
   b = 1.4;
   wn = 1.8;
   zeta = 0.9;
   am = { {0, 1}, {-wn^2, -2 zeta wn} };
   bm = {0, wn^2};
   b0 = {0, 1};
   gamma = DiagonalMatrix[{6.0, 3.0, 6.0}];
  ]], "Input"],

 Cell["Solve the Lyapunov equation am^T.p + p.am = -I.", "Text"],

 Cell[BoxData[
  ToBoxes[
   pRules = First@Solve[
      Transpose[am].{ {p11, p12}, {p12, p22} } +
        { {p11, p12}, {p12, p22} }.am == -IdentityMatrix[2],
      {p11, p12, p22}
    ];
   p = N[{ {p11, p12}, {p12, p22} } /. pRules];
   MatrixForm[p]
  ]], "Input"],

 Cell[BoxData[
  ToBoxes[
   r[t_?NumericQ] := Piecewise[{
      {0.0, t < 2.0},
      {1.0 + 0.35 Sin[0.8 t], t < 10.0},
      {-0.5 + 0.40 Sin[1.2 t], t < 16.0}
    }, 0.8 + 0.25 Sin[0.55 t]];
  ]], "Input"],

 Cell["The augmented state contains plant states, reference-model states, and the three direct controller gains.", "Text"],

 Cell[BoxData[
  ToBoxes[
   solution = NDSolveValue[
     {
      x1'[t] == x2[t],
      x2'[t] == -a0 x1[t] - a1 x2[t] +
        b (theta1[t] x1[t] + theta2[t] x2[t] + thetaR[t] r[t]),

      xm1'[t] == xm2[t],
      xm2'[t] == -wn^2 xm1[t] - 2 zeta wn xm2[t] + wn^2 r[t],

      theta1'[t] == -gamma[[1, 1]] x1[t] tuningError[t],
      theta2'[t] == -gamma[[2, 2]] x2[t] tuningError[t],
      thetaR'[t] == -gamma[[3, 3]] r[t] tuningError[t],

      tuningError[t] ==
        ({x1[t] - xm1[t], x2[t] - xm2[t]}.p.b0),

      x1[0] == 0, x2[0] == 0,
      xm1[0] == 0, xm2[0] == 0,
      theta1[0] == 0, theta2[0] == 0, thetaR[0] == 0
     },
     {x1, x2, xm1, xm2, theta1, theta2, thetaR},
     {t, 0, 20},
     Method -> {"EquationSimplification" -> "Residual"}
   ];
  ]], "Input"],

 Cell[BoxData[
  ToBoxes[
   {x1f, x2f, xm1f, xm2f, theta1f, theta2f, thetaRf} = solution;
   idealGains = {
     (a0 - wn^2)/b,
     (a1 - 2 zeta wn)/b,
     wn^2/b
   };
   finalGains = {theta1f[20], theta2f[20], thetaRf[20]};
   rmsError = Sqrt[
     NIntegrate[(x1f[t] - xm1f[t])^2, {t, 0, 20}]/20
   ];
   {rmsError, finalGains, idealGains}
  ]], "Input"],

 Cell[BoxData[
  ToBoxes[
   Plot[
    Evaluate[{x1f[t], xm1f[t], r[t]}],
    {t, 0, 20},
    PlotLegends -> {"Plant output x1", "Reference output xm1", "Command r"},
    AxesLabel -> {"Time (s)", "Output"},
    GridLines -> Automatic,
    PlotRange -> All
   ]
  ]], "Input"],

 Cell[BoxData[
  ToBoxes[
   Plot[
    Evaluate[{theta1f[t], theta2f[t], thetaRf[t]}],
    {t, 0, 20},
    PlotLegends -> {"theta1", "theta2", "thetaR"},
    AxesLabel -> {"Time (s)", "Direct controller gains"},
    GridLines -> Automatic,
    PlotRange -> All
   ]
  ]], "Input"],

 Cell[BoxData[
  ToBoxes[
   data = Table[
     {
      tt, r[tt], x1f[tt], x2f[tt], xm1f[tt], xm2f[tt],
      theta1f[tt], theta2f[tt], thetaRf[tt],
      theta1f[tt] x1f[tt] + theta2f[tt] x2f[tt] + thetaRf[tt] r[tt]
     },
     {tt, 0, 20, 0.01}
   ];
   Export[
    "Chapter13_Lesson2_results_mathematica.csv",
    Prepend[data, {"time", "reference", "x1", "x2", "xm1", "xm2",
      "theta1", "theta2", "theta_r", "control"}]
   ];
  ]], "Input"]
},
WindowTitle -> "Chapter13_Lesson2",
StyleDefinitions -> "Default.nb"
]

Several observations should be made when running the five implementations:

1. The plant output approaches the reference-model output even though the controller never estimates \(a_0,a_1,b\).
2. The adaptive gains need not converge rapidly to the displayed ideal gains. Tracking and unique parameter convergence are different objectives.
3. The command contains multiple frequency components to provide a richer regressor than a single constant step.
4. The same physical law is represented in every language: \(u=\hat{\boldsymbol{\theta} }^T\boldsymbol{\phi}\).
5. Numerical integration is not part of the adaptive-control theorem. It is an implementation approximation whose step size must resolve both plant and adaptation dynamics.

Problem 1: Ideal gains for a second-order plant

Consider \(\dot x_1=x_2\) and \(\dot x_2=-4x_1-3x_2+2u\). Let \(u=k_1x_1+k_2x_2+k_rr\). Find the ideal gains that match a reference model with \(\omega_n=3\) and \(\zeta=0.8\).

Solution:

The matching equations are

\[ -4+2k_1^{*}=-9, \qquad -3+2k_2^{*}=-4.8, \qquad 2k_r^{*}=9. \]

Thus

\[ k_1^{*}=-2.5, \qquad k_2^{*}=-0.9, \qquad k_r^{*}=4.5. \]

Problem 2: Reference model from settling-time and overshoot requirements

Design a second-order unit-DC-gain reference model with approximately five-percent overshoot and two-percent settling time \(T_s=2\) seconds.

Solution:

The overshoot relation is

\[ 0.05= \exp\left( -\frac{\pi\zeta}{\sqrt{1-\zeta^2} } \right). \]

Solving gives approximately \(\zeta=0.6901\). From

\[ T_s\approx\frac{4}{\zeta\omega_n}=2, \]

we obtain

\[ \omega_n\approx\frac{2}{0.6901}\approx2.898. \]

Therefore, a suitable model is

\[ G_m(s)= \frac{8.398} {s^2+4s+8.398}. \]

The coefficient of \(s\) is approximately \(2\zeta\omega_n=4\).

Problem 3: Structural feasibility of matching

Let

\[ B_p= \begin{bmatrix}0\\1\end{bmatrix}, \qquad A_m-A_p= \begin{bmatrix} 1&0\\ 2&3 \end{bmatrix}. \]

Determine whether a row vector \(\mathbf{K}_x^{*T}\) can satisfy \(B_p\mathbf{K}_x^{*T}=A_m-A_p\).

Solution:

Every matrix of the form \(B_p\mathbf{K}_x^{*T}\) has a zero first row:

\[ B_p \begin{bmatrix}k_1^{*}&k_2^{*}\end{bmatrix} = \begin{bmatrix} 0&0\\ k_1^{*}&k_2^{*} \end{bmatrix}. \]

The desired difference has first row \([1\;\;0]\), so it does not lie in the allowable range. Exact model matching is impossible with this input direction and controller structure.

Problem 4: Equivalent scaled parameterization

Let \(\boldsymbol{\phi}=[x_1\;\;x_2\;\;r]^T\) and choose \(D=\operatorname{diag}(10,0.1,2)\). Write the scaled regressor and parameter vector that preserve the control input.

Solution:

\[ \bar{\boldsymbol{\phi} } = D^{-1}\boldsymbol{\phi} = \begin{bmatrix} 0.1x_1\\ 10x_2\\ 0.5r \end{bmatrix}, \qquad \bar{\boldsymbol{\theta} } = D^T\boldsymbol{\theta} = \begin{bmatrix} 10\theta_1\\ 0.1\theta_2\\ 2\theta_r \end{bmatrix}. \]

Then

\[ \bar{\boldsymbol{\theta} }^T\bar{\boldsymbol{\phi} } = \theta_1x_1+\theta_2x_2+\theta_rr. \]

If the original adaptation matrix is \(\Gamma\), the equivalent scaled matrix is \(\bar\Gamma=D^T\Gamma D\).

Problem 5: Redundant regressor components

Suppose a controller uses \(\boldsymbol{\phi}=[x_1\;\;2x_1\;\;r]^T\). Show that its first two gains cannot be uniquely determined from the control signal.

Solution:

The control law is

\[ u=\theta_1x_1+2\theta_2x_1+\theta_3r = \left(\theta_1+2\theta_2\right)x_1+\theta_3r. \]

Only the combination \(\theta_1+2\theta_2\) affects the controller. Let \(\mathbf{v}=[-2\;\;1\;\;0]^T\). Then

\[ \mathbf{v}^T\boldsymbol{\phi} = -2x_1+2x_1=0. \]

Therefore, \(\boldsymbol{\theta}+\alpha\mathbf{v}\) produces the same control for every scalar \(\alpha\). The parameterization is redundant, and the regressor Gramian is singular.

Problem 6: Deriving the parameterized error equation

Starting from the general plant \(\dot{\mathbf{x} }=A_p\mathbf{x}+B_pu\), controller \(u=\hat{\mathbf{K} }_x^T\mathbf{x}+\hat K_rr\), and matching equations, derive the tracking-error dynamics.

Solution:

Write

\[ \hat{\mathbf{K} }_x = \mathbf{K}_x^{*}+\tilde{\mathbf{K} }_x, \qquad \hat K_r=K_r^{*}+\tilde K_r. \]

Then

\[ \begin{aligned} \dot{\mathbf{x} } &= \left(A_p+B_p\mathbf{K}_x^{*T}\right)\mathbf{x} +B_pK_r^{*}r\\ &\quad+ B_p\left( \tilde{\mathbf{K} }_x^T\mathbf{x} +\tilde K_rr \right)\\ &= A_m\mathbf{x}+B_mr +B_p\tilde{\boldsymbol{\theta} }^T\boldsymbol{\phi}. \end{aligned} \]

Subtracting \(\dot{\mathbf{x} }_m=A_m\mathbf{x}_m+B_mr\) gives

\[ \dot{\mathbf{e} } = A_m\mathbf{e} +B_p\tilde{\boldsymbol{\theta} }^T\boldsymbol{\phi}. \]

Direct self-tuning regulation requires a parameterized controller, not an explicit online plant model. The standard form \(u=\hat{\boldsymbol{\theta} }^T\boldsymbol{\phi}\) separates adjustable gains from measurable basis signals. A Hurwitz reference model specifies the desired closed-loop dynamics. Ideal gains are defined by model-matching equations, and subtraction of the reference-model dynamics yields an error equation in which parameter error appears through \(\tilde{\boldsymbol{\theta} }^T\boldsymbol{\phi}\). Scaling changes numerical coordinates but not the physical control signal, while redundant or insufficiently excited regressors prevent unique gain convergence. Lesson 3 uses these equations to construct and prove stable Lyapunov-based update laws.