After completing this lesson, the student should be able to:

1. distinguish parametric uncertainty from additive and multiplicative unmodeled dynamics;
2. compute frozen-time loop transfer functions, sensitivity functions, and classical gain/phase margins for an adaptive controller at a fixed parameter estimate;
3. explain why a positive frozen-time phase margin is not, by itself, a stability proof for the full adaptive system;
4. derive a scalar small-gain condition for multiplicative uncertainty;
5. relate leakage, projection, normalization, and adaptation-gain selection to boundedness and robustness reserve;
6. interpret simulation metrics that compare a robustified adaptive law with an aggressive unmodified law.

The lesson assumes knowledge of linear feedback, Bode plots, Nyquist stability, gain and phase margins, sensitivity functions, and Lyapunov stability. It uses adaptive-control concepts introduced in Chapters 1–21, particularly MRAC, projection, normalization, persistent excitation, composite adaptation, and the robust modifications reviewed in Lessons 1 and 2 of this chapter.

Adaptive control and robust control begin with different descriptions of uncertainty. In a linearly parameterized adaptive model,

\[ \dot{x}=A x+B\left(u+\theta^{*\mathsf T}\phi(x,t)\right)+d(t), \]

the unknown constant vector \( \theta^* \) has a known structural location. The adaptive law estimates or compensates for this structured uncertainty online. By contrast, robust analysis often describes an entire set of possible plants without trying to identify the exact member:

\[ \mathcal{P}=\left\{P(s,\Delta):\Delta\in\mathcal{D}\right\}. \]

The adaptive question is: Can the controller update its parameters so that tracking remains satisfactory despite unknown structured coefficients? The robust question is: Does the closed loop remain stable and performant for every admissible uncertainty?

A practical robust adaptive controller must answer both. Adaptation handles uncertainty that is learnable through the assumed regressor structure; robustification limits the consequences of uncertainty that is not captured, not sufficiently excited, time varying, noisy, delayed, or high frequency.

flowchart TD
  R["Reference r"] --> M["Reference model"]
  R --> C["Adaptive controller C(theta_hat)"]
  M --> E["Tracking error e"]
  P["Actual plant P0 plus uncertainty"] --> E
  E --> A["Adaptive law: \ngradient plus \nrobust modification"]
  A --> C
  C --> P
  U["Unmodeled dynamics, delay, \nnoise, disturbance"] --> P
  F["Frozen-time analysis: \nL, S, T, gain margin, \nphase margin"] --> C
  F --> P
        

3.1 Parametric uncertainty

For a state-space plant with uncertain coefficients,

\[ \dot{x}=\left(A_0+\sum_{i=1}^{p}\theta_i^* A_i\right)x+B u, \]

the uncertainty is finite dimensional and structured. When the matching conditions hold, the unknown coefficients can often be moved into a regressor:

\[ \dot{x}=A_0x+B\left(u+\theta^{*\mathsf T}\phi(x,t)\right). \]

This is the natural domain of MRAC and related adaptive laws. Parameter convergence still requires excitation, but bounded tracking may be proved without exact identification.

3.2 Additive dynamic uncertainty

A reduced-order nominal transfer function may omit stable high-frequency dynamics:

\[ P(s)=P_0(s)+W_a(s)\Delta_a(s),\qquad \|\Delta_a\|_\infty\leq 1. \]

The weight \( W_a(s) \) expresses the frequency-dependent size of the modeling error. Additive uncertainty is convenient when the absolute model discrepancy is meaningful.

3.3 Multiplicative dynamic uncertainty

Relative model error is represented by

\[ P(s)=P_0(s)\left[1+W_m(s)\Delta_m(s)\right],\qquad \|\Delta_m\|_\infty\leq 1. \]

Multiplicative uncertainty is especially natural for neglected actuator, sensor, structural, and transport dynamics whose relative effect grows with frequency. The uncertainty weight is usually small in the control band and larger near or above the neglected-dynamics frequencies.

3.4 Disturbances, noise, and delay

An output disturbance \( d_y \), measurement noise \( n \), and delay \( e^{-s\tau} \) are not equivalent to a constant parameter error. Persistent disturbances can drive a pure gradient adaptive law even when the tracking error cannot be eliminated. This mechanism motivates leakage, dead zones, projection, and bandwidth limitation.

Let the controller depend on an adaptive parameter vector:

\[ C(s,\hat{\theta}(t)). \]

At a fixed time \( t_0 \), freeze the parameter estimate:

\[ \bar{\theta}=\hat{\theta}(t_0),\qquad L_f(s;t_0)=P(s)C(s,\bar{\theta}). \]

The associated frozen sensitivity and complementary-sensitivity functions are

\[ S_f(s;t_0)=\frac{1}{1+L_f(s;t_0)},\qquad T_f(s;t_0)=\frac{L_f(s;t_0)}{1+L_f(s;t_0)}. \]

These retain the standard identity \( S_f+T_f=1 \). The usual interpretations remain useful:

• small \( |S_f| \) at low frequency implies disturbance rejection and reduced sensitivity to plant variations;
• small \( |T_f| \) at high frequency reduces amplification of measurement noise and multiplicative uncertainty;
• a large peak \( M_s=\|S_f\|_\infty \) signals fragile feedback geometry even when the nominal frozen loop is stable.

Gain and phase margins are also defined from \( L_f(j\omega;t_0) \). If \( \omega_c \) is the gain-crossover frequency,

\[ |L_f(j\omega_c;t_0)|=1,\qquad \mathrm{PM}(t_0)=180^\circ+\angle L_f(j\omega_c;t_0). \]

If \( \omega_\pi \) is the phase-crossover frequency,

\[ \angle L_f(j\omega_\pi;t_0)=-180^\circ,\qquad \mathrm{GM}(t_0)=\frac{1}{|L_f(j\omega_\pi;t_0)|}. \]

Because \( \hat{\theta}(t) \) moves, the quantities \( S_f,T_f,\mathrm{PM},\mathrm{GM} \) form trajectories rather than single constants. A robust adaptive design seeks acceptable frozen margins over the entire admissible parameter set, not merely at the final estimate.

Classical margins are exact properties of a linear time-invariant feedback loop. An adaptive loop is nonlinear and time varying because \( \hat{\theta}(t) \) is generated by the closed-loop signals. Freezing the parameters discards the dynamics

\[ \dot{\hat{\theta}}=\mathcal{A}\!\left(e,\phi,\hat{\theta}\right). \]

The frozen system can therefore be stable for every fixed \( \bar{\theta} \) while the time-varying interconnection is poorly damped or even unstable. The missing mechanism is the energy exchange between the plant state, tracking error, and adaptation state.

A useful decomposition is

\[ u(t)=C(\hat{\theta}(t))y(t) =C(\bar{\theta})y(t) +\left[C(\hat{\theta}(t))-C(\bar{\theta})\right]y(t). \]

The first term belongs to a frozen LTI loop. The second is an endogenous, signal-dependent perturbation. Its size depends on both parameter motion and the plant signal:

\[ \|\delta u(t)\| \leq \left\|\frac{\partial C}{\partial\theta}\right\| \|\hat{\theta}(t)-\bar{\theta}\|\,\|y(t)\| +\mathcal{O}\!\left(\|\hat{\theta}-\bar{\theta}\|^2\right). \]

Consequently, frozen margins should be used together with:

1. a Lyapunov or input-output boundedness argument;
2. bounds on the admissible parameter set;
3. bounds on the parameter rate;
4. frequency separation from neglected dynamics;
5. monitoring of control effort, noise amplification, and saturation.

Consider the frozen nominal loop

\[ L_0(s)=P_0(s)C(s,\bar{\theta}), \qquad T_0(s)=\frac{L_0(s)}{1+L_0(s)}. \]

Let the actual plant be

\[ P(s)=P_0(s)\left[1+W_m(s)\Delta_m(s)\right], \qquad \|\Delta_m\|_\infty\leq 1. \]

The uncertain characteristic equation is

\[ 1+P_0C+P_0CW_m\Delta_m=0. \]

Factor the nominal characteristic term:

\[ \left(1+L_0\right) \left(1+T_0W_m\Delta_m\right)=0. \]

If the nominal frozen loop is internally stable, a sufficient robust stability condition follows from the small-gain theorem:

\[ \|T_0W_m\|_\infty<1. \]

Pointwise, the condition asks that

\[ |T_0(j\omega)|\,|W_m(j\omega)|<1 \qquad \text{for all }\omega. \]

This equation provides a direct qualitative lesson for adaptive control: if adaptation increases the effective feedback gain and pushes the crossover frequency toward neglected actuator or structural modes, then \( |T_0| \) remains large where \( |W_m| \) is no longer small. Robustness reserve is lost even though low-frequency tracking may temporarily improve.

Define the frozen uncertainty index

\[ \rho(\bar{\theta}) = \sup_{\omega} |T_0(j\omega;\bar{\theta})W_m(j\omega)|. \]

Then \( \rho<1 \) is a sufficient frozen-loop certificate for the stated uncertainty set. In an adaptive implementation, monitor \( \rho(\hat{\theta}(t)) \) over the projected parameter set or along the actual parameter trajectory.

Consider the matched error model

\[ \dot{e}=A_m e+B\tilde{\theta}^{\mathsf T}\phi+B d(t), \]

where \( A_m \) is Hurwitz, \( \tilde{\theta}=\hat{\theta}-\theta^* \), and \( \|d(t)\|\leq d_{\max} \). Let

\[ A_m^{\mathsf T}P+PA_m=-Q,\qquad P=P^{\mathsf T}>0,\quad Q=Q^{\mathsf T}>0. \]

Use the robustified adaptive law

\[ \dot{\hat{\theta}} =-\Gamma\phi e^{\mathsf T}PB-\sigma\Gamma\hat{\theta}, \qquad \Gamma=\Gamma^{\mathsf T}>0,\quad \sigma>0. \]

With the Lyapunov function

\[ V=e^{\mathsf T}Pe+ \tilde{\theta}^{\mathsf T}\Gamma^{-1}\tilde{\theta}, \]

differentiation and cancellation of the matched cross term give

\[ \dot{V} = -e^{\mathsf T}Qe +2e^{\mathsf T}PBd -2\sigma\tilde{\theta}^{\mathsf T}\hat{\theta}. \]

Since \( \hat{\theta}=\tilde{\theta}+\theta^* \),

\[ -2\sigma\tilde{\theta}^{\mathsf T}\hat{\theta} = -2\sigma\|\tilde{\theta}\|^2 -2\sigma\tilde{\theta}^{\mathsf T}\theta^*. \]

Apply Young's inequality twice. For arbitrary \( \varepsilon_d>0 \) and \( \varepsilon_\theta>0 \),

\[ 2e^{\mathsf T}PBd \leq \varepsilon_d\|e\|^2+ \frac{\|PB\|^2}{\varepsilon_d}d_{\max}^2, \]

\[ 2\sigma|\tilde{\theta}^{\mathsf T}\theta^*| \leq \sigma\varepsilon_\theta\|\tilde{\theta}\|^2+ \frac{\sigma}{\varepsilon_\theta}\|\theta^*\|^2. \]

Choose \( 0<\varepsilon_d<\lambda_{\min}(Q) \) and \( 0<\varepsilon_\theta<2 \). Then

\[ \dot{V} \leq -c_e\|e\|^2-c_\theta\|\tilde{\theta}\|^2+c_0, \]

where

\[ c_e=\lambda_{\min}(Q)-\varepsilon_d,\qquad c_\theta=\sigma(2-\varepsilon_\theta), \]

\[ c_0= \frac{\|PB\|^2}{\varepsilon_d}d_{\max}^2+ \frac{\sigma}{\varepsilon_\theta}\|\theta^*\|^2. \]

Hence the state and parameter errors are uniformly ultimately bounded. A conservative residual-set description is

\[ c_e\|e\|^2+c_\theta\|\tilde{\theta}\|^2\leq c_0. \]

This result is the nonlinear adaptive counterpart of a robustness reserve. The controller does not promise exact rejection of all nonparametric uncertainty. It promises that bounded uncertainty and disturbances produce a bounded residual response whose size can be shaped by \( Q,\Gamma,\sigma \), signal normalization, and projection limits.

A common misconception is that increasing the adaptation gain \( \Gamma \) always improves performance. For the gradient law

\[ \dot{\hat{\theta}}=-\Gamma\phi e, \]

the parameter-rate bound is

\[ \|\dot{\hat{\theta}}\| \leq \|\Gamma\|\,\|\phi\|\,\|e\|. \]

Larger \( \Gamma \) therefore increases the possible bandwidth and amplitude of controller variation. In the presence of measurement noise or an unmodeled high-frequency mode, this can create:

• rapid parameter oscillation;
• high-frequency control activity;
• movement of the frozen crossover frequency;
• reduced phase margin and delay margin;
• saturation, which invalidates the nominal matching argument.

A useful engineering separation is

\[ \omega_{\mathrm{cmd}} \ll \omega_{\mathrm{cl}} \ll \omega_{\mathrm{unmodeled}}, \]

together with an adaptation bandwidth that does not inject substantial energy near \( \omega_{\mathrm{unmodeled}} \). This is not a universal theorem expressed by one scalar inequality; it is a design discipline checked through parameter-rate bounds, filtered regressors, actuator models, frozen-loop margins, and nonlinear simulation.

Projection gives a compact admissible set \( \Omega_\theta \). One can then require

\[ \inf_{\bar{\theta}\in\Omega_\theta}\mathrm{PM}(\bar{\theta}) \geq \mathrm{PM}_{\min}, \]

\[ \inf_{\bar{\theta}\in\Omega_\theta}\mathrm{GM}(\bar{\theta}) \geq \mathrm{GM}_{\min}, \qquad \sup_{\bar{\theta}\in\Omega_\theta}\rho(\bar{\theta})<1. \]

These are frozen-set design checks. The Lyapunov proof remains responsible for the full time-varying closed loop.

The implementations below use the scalar plant

\[ \dot{x}=-a x+z_2+d(t), \]

\[ \tau_1\dot{z}_1=-z_1+u,\qquad \tau_2\dot{z}_2=-z_2+z_1. \]

The adaptive controller is designed from the reduced model that assumes \( z_2=u \). The two stable actuator poles are therefore unmodeled dynamics. The reference model is

\[ \dot{x}_m=-a_mx_m+b_mr. \]

The control law and adaptive law are

\[ u=-\hat{k}x+\hat{k}_r r, \]

\[ \dot{\hat{k}}=\gamma x e-\sigma\hat{k},\qquad \dot{\hat{k}}_r=-\gamma r e-\sigma\hat{k}_r, \qquad e=x-x_m. \]

Both estimates are projected onto \( [0,10] \), and the control is limited to \( [-30,30] \). Two cases are compared:

1. Robustified: \( \gamma=8,\ \sigma=0.15 \).
2. Aggressive plain adaptation: \( \gamma=60,\ \sigma=0 \).

For frozen regulation feedback \( u=-k x \), the actual loop is

\[ L_f(s;k)= \frac{k} {(s+a)(1+\tau_1s)(1+\tau_2s)}. \]

The reduced nominal plant is \( P_0(s)=1/(s+a) \), and the exact relative error caused by the omitted poles is

\[ \Delta_m(s)= \frac{1}{(1+\tau_1s)(1+\tau_2s)}-1. \]

The programs estimate phase margin, gain margin, and \( \sup_\omega|\Delta_mT_0| \) at the largest adaptive feedback gain reached during each simulation.

flowchart TD
  A["Specify reduced model and matched parameters"] --> B["List nonparametric uncertainty: poles, delay, noise, saturation"]
  B --> C["Choose adaptive law \nand robust modification"]
  C --> D["Prove boundedness or \nultimate boundedness"]
  D --> E["Define projected parameter set"]
  E --> F["Freeze parameters over \nthe admissible set"]
  F --> G["Compute L, S, T, gain margin, \nphase margin, uncertainty index"]
  G --> H{"Margins and uncertainty \nreserve acceptable?"}
  H -->|"No"| I["Reduce adaptation rate, \nfilter signals, tighten projection, \nredesign baseline loop"]
  I --> C
  H -->|"Yes"| J["Run nonlinear simulations \nwith uncertainty and saturation"]
  J --> K["Monitor tracking, control effort, \nparameter rate, margins, \nand residual set"]
        

The order matters. Frequency-domain checks do not replace the adaptive stability argument, and the Lyapunov proof does not automatically guarantee useful classical margins or acceptable actuator activity.

This standard-library program simulates both cases and computes frozen-loop metrics by direct frequency gridding.

Chapter22_Lesson3.py

# Chapter22_Lesson3.py
# Qualitative connection between robust-control margins and robust adaptive control.
# Standard-library implementation: no third-party packages are required.

import cmath
import math
from dataclasses import dataclass


@dataclass
class CaseResult:
    name: str
    rms_error: float
    max_error: float
    rms_control: float
    max_k: float
    max_kr: float
    final_k: float
    final_kr: float
    phase_margin_deg_at_max_k: float
    gain_margin_db_at_max_k: float
    small_gain_index_at_max_k: float


def reference(t: float) -> float:
    return 1.0 + 0.30 * math.sin(0.50 * t) + 0.15 * math.sin(1.70 * t)


def disturbance(t: float) -> float:
    return 0.05 * math.sin(8.0 * t)


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


def frozen_metrics(k: float, a: float, tau1: float, tau2: float):
    """Estimate frozen-controller margins and a multiplicative-uncertainty index."""
    n = 24000
    log_w_min = math.log(1.0e-3)
    log_w_max = math.log(1.0e4)

    previous_w = None
    previous_mag_db = None
    previous_phase_deg = None
    gain_cross = None
    phase_cross = None
    max_small_gain = 0.0

    for i in range(n):
        alpha = i / (n - 1)
        w = math.exp(log_w_min + alpha * (log_w_max - log_w_min))
        s = 1j * w

        p0 = 1.0 / (s + a)
        parasitic = 1.0 / ((1.0 + tau1 * s) * (1.0 + tau2 * s))
        p = p0 * parasitic
        loop = k * p
        mag = abs(loop)
        mag_db = 20.0 * math.log10(max(mag, 1.0e-300))
        phase_deg = math.degrees(cmath.phase(loop))

        # For this stable three-pole loop the principal phase decreases from 0 to -270 deg.
        if phase_deg > 0.0:
            phase_deg -= 360.0

        delta_m = parasitic - 1.0
        loop0 = k * p0
        t0 = loop0 / (1.0 + loop0)
        max_small_gain = max(max_small_gain, abs(delta_m * t0))

        if previous_w is not None:
            if gain_cross is None and previous_mag_db >= 0.0 and mag_db < 0.0:
                fraction = previous_mag_db / (previous_mag_db - mag_db)
                log_wc = math.log(previous_w) + fraction * (math.log(w) - math.log(previous_w))
                phase_wc = previous_phase_deg + fraction * (phase_deg - previous_phase_deg)
                gain_cross = (math.exp(log_wc), 180.0 + phase_wc)

            if phase_cross is None and previous_phase_deg > -180.0 and phase_deg <= -180.0:
                fraction = (previous_phase_deg + 180.0) / (previous_phase_deg - phase_deg)
                log_wp = math.log(previous_w) + fraction * (math.log(w) - math.log(previous_w))
                mag_db_wp = previous_mag_db + fraction * (mag_db - previous_mag_db)
                phase_cross = (math.exp(log_wp), -mag_db_wp)

        previous_w = w
        previous_mag_db = mag_db
        previous_phase_deg = phase_deg

    phase_margin = float("nan") if gain_cross is None else gain_cross[1]
    gain_margin_db = float("inf") if phase_cross is None else phase_cross[1]
    return phase_margin, gain_margin_db, max_small_gain


def run_case(name: str, gamma: float, sigma: float) -> CaseResult:
    # True plant: x_dot = -a*x + z2 + d.
    # The controller is designed as if z2 = u, while z1 and z2 are unmodeled
    # actuator poles. This is a reduced-order adaptive-control benchmark.
    a = 1.0
    am = 2.0
    bm = 2.0
    tau1 = 0.05
    tau2 = 0.02

    dt = 5.0e-4
    final_time = 20.0
    steps = int(final_time / dt)

    x = 0.0
    z1 = 0.0
    z2 = 0.0
    xm = 0.0
    k = 0.0
    kr = 0.0

    error_energy = 0.0
    control_energy = 0.0
    max_error = 0.0
    max_k = 0.0
    max_kr = 0.0

    for step in range(steps):
        t = step * dt
        r = reference(t)
        d = disturbance(t)
        e = x - xm

        # u = -k*x + kr*r. Projection protects against parameter escape.
        u = clamp(-k * x + kr * r, -30.0, 30.0)

        dx = -a * x + z2 + d
        dz1 = (-z1 + u) / tau1
        dz2 = (-z2 + z1) / tau2
        dxm = -am * xm + bm * r

        # Gradient adaptation plus sigma-modification.
        dk = gamma * x * e - sigma * k
        dkr = -gamma * r * e - sigma * kr

        x += dt * dx
        z1 += dt * dz1
        z2 += dt * dz2
        xm += dt * dxm
        k = clamp(k + dt * dk, 0.0, 10.0)
        kr = clamp(kr + dt * dkr, 0.0, 10.0)

        error_energy += e * e * dt
        control_energy += u * u * dt
        max_error = max(max_error, abs(e))
        max_k = max(max_k, k)
        max_kr = max(max_kr, kr)

    phase_margin, gain_margin_db, small_gain_index = frozen_metrics(
        max(max_k, 1.0e-8), a, tau1, tau2
    )

    return CaseResult(
        name=name,
        rms_error=math.sqrt(error_energy / final_time),
        max_error=max_error,
        rms_control=math.sqrt(control_energy / final_time),
        max_k=max_k,
        max_kr=max_kr,
        final_k=k,
        final_kr=kr,
        phase_margin_deg_at_max_k=phase_margin,
        gain_margin_db_at_max_k=gain_margin_db,
        small_gain_index_at_max_k=small_gain_index,
    )


def format_number(value: float) -> str:
    if math.isnan(value):
        return "N/A"
    if math.isinf(value):
        return "infinity"
    return f"{value:.6f}"


def main() -> None:
    cases = [
        run_case("robustified", gamma=8.0, sigma=0.15),
        run_case("aggressive_plain", gamma=60.0, sigma=0.0),
    ]

    print(
        "case,rms_error,max_error,rms_control,max_k,max_kr,final_k,final_kr,"
        "phase_margin_deg_at_max_k,gain_margin_db_at_max_k,small_gain_index_at_max_k"
    )
    for result in cases:
        print(
            f"{result.name},{result.rms_error:.6f},{result.max_error:.6f},"
            f"{result.rms_control:.6f},{result.max_k:.6f},{result.max_kr:.6f},"
            f"{result.final_k:.6f},{result.final_kr:.6f},"
            f"{format_number(result.phase_margin_deg_at_max_k)},"
            f"{format_number(result.gain_margin_db_at_max_k)},"
            f"{result.small_gain_index_at_max_k:.6f}"
        )

    print("\nInterpretation:")
    print("- Classical margins are frozen-time diagnostics, not a proof for the adaptive loop.")
    print("- A smaller peak adaptive feedback gain generally preserves more phase/gain reserve.")
    print("- The small-gain index below one is a useful sufficient check for this uncertainty model.")
    print("- Sigma-modification and projection trade exact parameter convergence for boundedness.")


if __name__ == "__main__":
    main()

The C++17 implementation uses std::complex for the frequency response and explicit Euler integration for transparency.

Chapter22_Lesson3.cpp

// Chapter22_Lesson3.cpp
// Qualitative connection between robust-control margins and robust adaptive control.
// Build: g++ -std=c++17 -O2 Chapter22_Lesson3.cpp -o Chapter22_Lesson3

#include <algorithm>
#include <cmath>
#include <complex>
#include <iomanip>
#include <iostream>
#include <limits>
#include <string>
#include <sstream>
#include <vector>

struct CaseResult {
    std::string name;
    double rmsError;
    double maxError;
    double rmsControl;
    double maxK;
    double maxKr;
    double finalK;
    double finalKr;
    double phaseMarginDeg;
    double gainMarginDb;
    double smallGainIndex;
};

double referenceSignal(double t) {
    return 1.0 + 0.30 * std::sin(0.50 * t) + 0.15 * std::sin(1.70 * t);
}

double disturbance(double t) {
    return 0.05 * std::sin(8.0 * t);
}

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

struct FrozenMetrics {
    double phaseMarginDeg;
    double gainMarginDb;
    double smallGainIndex;
};

FrozenMetrics frozenMetrics(double k, double a, double tau1, double tau2) {
    const int n = 24000;
    const double logWMin = std::log(1.0e-3);
    const double logWMax = std::log(1.0e4);
    const std::complex<double> j(0.0, 1.0);

    bool havePrevious = false;
    bool foundGainCross = false;
    bool foundPhaseCross = false;
    double previousW = 0.0;
    double previousMagDb = 0.0;
    double previousPhaseDeg = 0.0;
    double phaseMargin = std::numeric_limits<double>::quiet_NaN();
    double gainMarginDb = std::numeric_limits<double>::infinity();
    double maxSmallGain = 0.0;

    for (int i = 0; i < n; ++i) {
        const double alpha = static_cast<double>(i) / static_cast<double>(n - 1);
        const double w = std::exp(logWMin + alpha * (logWMax - logWMin));
        const std::complex<double> s = j * w;
        const std::complex<double> p0 = 1.0 / (s + a);
        const std::complex<double> parasitic =
            1.0 / ((1.0 + tau1 * s) * (1.0 + tau2 * s));
        const std::complex<double> p = p0 * parasitic;
        const std::complex<double> loop = k * p;

        const double magnitude = std::abs(loop);
        const double magDb = 20.0 * std::log10(std::max(magnitude, 1.0e-300));
        double phaseDeg = std::arg(loop) * 180.0 / std::acos(-1.0);
        if (phaseDeg > 0.0) {
            phaseDeg -= 360.0;
        }

        const std::complex<double> deltaM = parasitic - 1.0;
        const std::complex<double> loop0 = k * p0;
        const std::complex<double> t0 = loop0 / (1.0 + loop0);
        maxSmallGain = std::max(maxSmallGain, std::abs(deltaM * t0));

        if (havePrevious) {
            if (!foundGainCross && previousMagDb >= 0.0 && magDb < 0.0) {
                const double fraction = previousMagDb / (previousMagDb - magDb);
                const double phaseAtCross =
                    previousPhaseDeg + fraction * (phaseDeg - previousPhaseDeg);
                phaseMargin = 180.0 + phaseAtCross;
                foundGainCross = true;
            }

            if (!foundPhaseCross && previousPhaseDeg > -180.0 && phaseDeg <= -180.0) {
                const double fraction =
                    (previousPhaseDeg + 180.0) / (previousPhaseDeg - phaseDeg);
                const double magDbAtCross =
                    previousMagDb + fraction * (magDb - previousMagDb);
                gainMarginDb = -magDbAtCross;
                foundPhaseCross = true;
            }
        }

        havePrevious = true;
        previousW = w;
        previousMagDb = magDb;
        previousPhaseDeg = phaseDeg;
        (void)previousW;
    }

    return {phaseMargin, gainMarginDb, maxSmallGain};
}

CaseResult runCase(const std::string& name, double gamma, double sigma) {
    const double a = 1.0;
    const double am = 2.0;
    const double bm = 2.0;
    const double tau1 = 0.05;
    const double tau2 = 0.02;
    const double dt = 5.0e-4;
    const double finalTime = 20.0;
    const int steps = static_cast<int>(finalTime / dt);

    double x = 0.0;
    double z1 = 0.0;
    double z2 = 0.0;
    double xm = 0.0;
    double k = 0.0;
    double kr = 0.0;

    double errorEnergy = 0.0;
    double controlEnergy = 0.0;
    double maxError = 0.0;
    double maxK = 0.0;
    double maxKr = 0.0;

    for (int step = 0; step < steps; ++step) {
        const double t = step * dt;
        const double r = referenceSignal(t);
        const double d = disturbance(t);
        const double e = x - xm;
        const double u = clampValue(-k * x + kr * r, -30.0, 30.0);

        const double dx = -a * x + z2 + d;
        const double dz1 = (-z1 + u) / tau1;
        const double dz2 = (-z2 + z1) / tau2;
        const double dxm = -am * xm + bm * r;
        const double dk = gamma * x * e - sigma * k;
        const double dkr = -gamma * r * e - sigma * kr;

        x += dt * dx;
        z1 += dt * dz1;
        z2 += dt * dz2;
        xm += dt * dxm;
        k = clampValue(k + dt * dk, 0.0, 10.0);
        kr = clampValue(kr + dt * dkr, 0.0, 10.0);

        errorEnergy += e * e * dt;
        controlEnergy += u * u * dt;
        maxError = std::max(maxError, std::abs(e));
        maxK = std::max(maxK, k);
        maxKr = std::max(maxKr, kr);
    }

    const FrozenMetrics metrics = frozenMetrics(std::max(maxK, 1.0e-8), a, tau1, tau2);
    return {
        name,
        std::sqrt(errorEnergy / finalTime),
        maxError,
        std::sqrt(controlEnergy / finalTime),
        maxK,
        maxKr,
        k,
        kr,
        metrics.phaseMarginDeg,
        metrics.gainMarginDb,
        metrics.smallGainIndex
    };
}

std::string formatNumber(double value) {
    if (std::isnan(value)) {
        return "N/A";
    }
    if (std::isinf(value)) {
        return "infinity";
    }
    std::ostringstream stream;
    stream << std::fixed << std::setprecision(6) << value;
    return stream.str();
}

int main() {
    const std::vector<CaseResult> cases = {
        runCase("robustified", 8.0, 0.15),
        runCase("aggressive_plain", 60.0, 0.0)
    };

    std::cout << std::fixed << std::setprecision(6);
    std::cout
        << "case,rms_error,max_error,rms_control,max_k,max_kr,final_k,final_kr,"
        << "phase_margin_deg_at_max_k,gain_margin_db_at_max_k,"
        << "small_gain_index_at_max_k\n";

    for (const CaseResult& result : cases) {
        std::cout << result.name << ","
                  << result.rmsError << ","
                  << result.maxError << ","
                  << result.rmsControl << ","
                  << result.maxK << ","
                  << result.maxKr << ","
                  << result.finalK << ","
                  << result.finalKr << ","
                  << formatNumber(result.phaseMarginDeg) << ","
                  << formatNumber(result.gainMarginDb) << ","
                  << result.smallGainIndex << "\n";
    }

    std::cout << "\nInterpretation:\n"
              << "- Classical margins are frozen-time diagnostics, not a proof for the adaptive loop.\n"
              << "- A smaller peak adaptive feedback gain generally preserves more phase/gain reserve.\n"
              << "- A small-gain index below one is sufficient for this uncertainty representation.\n"
              << "- Sigma-modification and projection exchange exact convergence for boundedness.\n";
    return 0;
}

Java has no complex-number type in the standard library, so a small immutable complex class is included from scratch.

Chapter22_Lesson3.java

// Chapter22_Lesson3.java
// Qualitative connection between robust-control margins and robust adaptive control.
// Build and run:
//   javac Chapter22_Lesson3.java
//   java Chapter22_Lesson3

public final class Chapter22_Lesson3 {
    private static final class Complex {
        final double re;
        final double im;

        Complex(double re, double im) {
            this.re = re;
            this.im = im;
        }

        Complex add(Complex other) {
            return new Complex(re + other.re, im + other.im);
        }

        Complex subtract(Complex other) {
            return new Complex(re - other.re, im - other.im);
        }

        Complex multiply(Complex other) {
            return new Complex(
                re * other.re - im * other.im,
                re * other.im + im * other.re
            );
        }

        Complex multiply(double scalar) {
            return new Complex(re * scalar, im * scalar);
        }

        Complex divide(Complex other) {
            double denominator = other.re * other.re + other.im * other.im;
            return new Complex(
                (re * other.re + im * other.im) / denominator,
                (im * other.re - re * other.im) / denominator
            );
        }

        double abs() {
            return Math.hypot(re, im);
        }

        double phaseDegrees() {
            return Math.toDegrees(Math.atan2(im, re));
        }
    }

    private static final class FrozenMetrics {
        final double phaseMarginDeg;
        final double gainMarginDb;
        final double smallGainIndex;

        FrozenMetrics(double phaseMarginDeg, double gainMarginDb, double smallGainIndex) {
            this.phaseMarginDeg = phaseMarginDeg;
            this.gainMarginDb = gainMarginDb;
            this.smallGainIndex = smallGainIndex;
        }
    }

    private static final class CaseResult {
        final String name;
        final double rmsError;
        final double maxError;
        final double rmsControl;
        final double maxK;
        final double maxKr;
        final double finalK;
        final double finalKr;
        final double phaseMarginDeg;
        final double gainMarginDb;
        final double smallGainIndex;

        CaseResult(
            String name,
            double rmsError,
            double maxError,
            double rmsControl,
            double maxK,
            double maxKr,
            double finalK,
            double finalKr,
            double phaseMarginDeg,
            double gainMarginDb,
            double smallGainIndex
        ) {
            this.name = name;
            this.rmsError = rmsError;
            this.maxError = maxError;
            this.rmsControl = rmsControl;
            this.maxK = maxK;
            this.maxKr = maxKr;
            this.finalK = finalK;
            this.finalKr = finalKr;
            this.phaseMarginDeg = phaseMarginDeg;
            this.gainMarginDb = gainMarginDb;
            this.smallGainIndex = smallGainIndex;
        }
    }

    private static double reference(double t) {
        return 1.0 + 0.30 * Math.sin(0.50 * t) + 0.15 * Math.sin(1.70 * t);
    }

    private static double disturbance(double t) {
        return 0.05 * Math.sin(8.0 * t);
    }

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

    private static FrozenMetrics frozenMetrics(double k, double a, double tau1, double tau2) {
        final int n = 24000;
        final double logWMin = Math.log(1.0e-3);
        final double logWMax = Math.log(1.0e4);

        boolean havePrevious = false;
        boolean foundGainCross = false;
        boolean foundPhaseCross = false;
        double previousMagDb = 0.0;
        double previousPhaseDeg = 0.0;
        double phaseMargin = Double.NaN;
        double gainMarginDb = Double.POSITIVE_INFINITY;
        double maxSmallGain = 0.0;

        for (int i = 0; i < n; ++i) {
            double alpha = (double) i / (double) (n - 1);
            double w = Math.exp(logWMin + alpha * (logWMax - logWMin));
            Complex s = new Complex(0.0, w);

            Complex p0 = new Complex(1.0, 0.0).divide(s.add(new Complex(a, 0.0)));
            Complex firstPole = new Complex(1.0, tau1 * w);
            Complex secondPole = new Complex(1.0, tau2 * w);
            Complex parasitic = new Complex(1.0, 0.0).divide(firstPole.multiply(secondPole));
            Complex plant = p0.multiply(parasitic);
            Complex loop = plant.multiply(k);

            double magnitude = loop.abs();
            double magDb = 20.0 * Math.log10(Math.max(magnitude, 1.0e-300));
            double phaseDeg = loop.phaseDegrees();
            if (phaseDeg > 0.0) {
                phaseDeg -= 360.0;
            }

            Complex deltaM = parasitic.subtract(new Complex(1.0, 0.0));
            Complex loop0 = p0.multiply(k);
            Complex t0 = loop0.divide(new Complex(1.0, 0.0).add(loop0));
            maxSmallGain = Math.max(maxSmallGain, deltaM.multiply(t0).abs());

            if (havePrevious) {
                if (!foundGainCross && previousMagDb >= 0.0 && magDb < 0.0) {
                    double fraction = previousMagDb / (previousMagDb - magDb);
                    double phaseAtCross =
                        previousPhaseDeg + fraction * (phaseDeg - previousPhaseDeg);
                    phaseMargin = 180.0 + phaseAtCross;
                    foundGainCross = true;
                }

                if (!foundPhaseCross && previousPhaseDeg > -180.0 && phaseDeg <= -180.0) {
                    double fraction =
                        (previousPhaseDeg + 180.0) / (previousPhaseDeg - phaseDeg);
                    double magDbAtCross =
                        previousMagDb + fraction * (magDb - previousMagDb);
                    gainMarginDb = -magDbAtCross;
                    foundPhaseCross = true;
                }
            }

            havePrevious = true;
            previousMagDb = magDb;
            previousPhaseDeg = phaseDeg;
        }

        return new FrozenMetrics(phaseMargin, gainMarginDb, maxSmallGain);
    }

    private static CaseResult runCase(String name, double gamma, double sigma) {
        final double a = 1.0;
        final double am = 2.0;
        final double bm = 2.0;
        final double tau1 = 0.05;
        final double tau2 = 0.02;
        final double dt = 5.0e-4;
        final double finalTime = 20.0;
        final int steps = (int) (finalTime / dt);

        double x = 0.0;
        double z1 = 0.0;
        double z2 = 0.0;
        double xm = 0.0;
        double k = 0.0;
        double kr = 0.0;

        double errorEnergy = 0.0;
        double controlEnergy = 0.0;
        double maxError = 0.0;
        double maxK = 0.0;
        double maxKr = 0.0;

        for (int step = 0; step < steps; ++step) {
            double t = step * dt;
            double r = reference(t);
            double d = disturbance(t);
            double e = x - xm;
            double u = clamp(-k * x + kr * r, -30.0, 30.0);

            double dx = -a * x + z2 + d;
            double dz1 = (-z1 + u) / tau1;
            double dz2 = (-z2 + z1) / tau2;
            double dxm = -am * xm + bm * r;
            double dk = gamma * x * e - sigma * k;
            double dkr = -gamma * r * e - sigma * kr;

            x += dt * dx;
            z1 += dt * dz1;
            z2 += dt * dz2;
            xm += dt * dxm;
            k = clamp(k + dt * dk, 0.0, 10.0);
            kr = clamp(kr + dt * dkr, 0.0, 10.0);

            errorEnergy += e * e * dt;
            controlEnergy += u * u * dt;
            maxError = Math.max(maxError, Math.abs(e));
            maxK = Math.max(maxK, k);
            maxKr = Math.max(maxKr, kr);
        }

        FrozenMetrics metrics = frozenMetrics(Math.max(maxK, 1.0e-8), a, tau1, tau2);
        return new CaseResult(
            name,
            Math.sqrt(errorEnergy / finalTime),
            maxError,
            Math.sqrt(controlEnergy / finalTime),
            maxK,
            maxKr,
            k,
            kr,
            metrics.phaseMarginDeg,
            metrics.gainMarginDb,
            metrics.smallGainIndex
        );
    }

    private static String formatNumber(double value) {
        if (Double.isNaN(value)) {
            return "N/A";
        }
        if (Double.isInfinite(value)) {
            return "infinity";
        }
        return String.format("%.6f", value);
    }

    public static void main(String[] args) {
        CaseResult[] cases = {
            runCase("robustified", 8.0, 0.15),
            runCase("aggressive_plain", 60.0, 0.0)
        };

        System.out.println(
            "case,rms_error,max_error,rms_control,max_k,max_kr,final_k,final_kr,"
            + "phase_margin_deg_at_max_k,gain_margin_db_at_max_k,"
            + "small_gain_index_at_max_k"
        );

        for (CaseResult result : cases) {
            System.out.printf(
                "%s,%.6f,%.6f,%.6f,%.6f,%.6f,%.6f,%.6f,%s,%s,%.6f%n",
                result.name,
                result.rmsError,
                result.maxError,
                result.rmsControl,
                result.maxK,
                result.maxKr,
                result.finalK,
                result.finalKr,
                formatNumber(result.phaseMarginDeg),
                formatNumber(result.gainMarginDb),
                result.smallGainIndex
            );
        }

        System.out.println("\nInterpretation:");
        System.out.println("- Classical margins are frozen-time diagnostics, not a proof for the adaptive loop.");
        System.out.println("- A smaller peak adaptive feedback gain generally preserves more phase/gain reserve.");
        System.out.println("- A small-gain index below one is sufficient for this uncertainty representation.");
        System.out.println("- Sigma-modification and projection exchange exact convergence for boundedness.");
    }
}

The MATLAB script uses only core language operations. The same differential equations map directly into Simulink Integrator blocks: four integrators for \( x,z_1,z_2,x_m \), two for \( \hat{k},\hat{k}_r \), Saturation blocks for control and projection, and MATLAB Function blocks for the update laws.

Chapter22_Lesson3.m

% Chapter22_Lesson3.m
% Qualitative connection between robust-control margins and robust adaptive control.
% This script uses only core MATLAB functions. Run it directly from the editor.

clear; clc;

robustified = runCase("robustified", 8.0, 0.15);
aggressive = runCase("aggressive_plain", 60.0, 0.0);
results = [robustified; aggressive];

fprintf(['case,rms_error,max_error,rms_control,max_k,max_kr,final_k,final_kr,' ...
         'phase_margin_deg_at_max_k,gain_margin_db_at_max_k,' ...
         'small_gain_index_at_max_k\n']);

for i = 1:numel(results)
    fprintf('%s,%.6f,%.6f,%.6f,%.6f,%.6f,%.6f,%.6f,%s,%s,%.6f\n', ...
        results(i).name, ...
        results(i).rmsError, ...
        results(i).maxError, ...
        results(i).rmsControl, ...
        results(i).maxK, ...
        results(i).maxKr, ...
        results(i).finalK, ...
        results(i).finalKr, ...
        formatNumber(results(i).phaseMarginDeg), ...
        formatNumber(results(i).gainMarginDb), ...
        results(i).smallGainIndex);
end

fprintf('\nInterpretation:\n');
fprintf('- Classical margins are frozen-time diagnostics, not a proof for the adaptive loop.\n');
fprintf('- A smaller peak adaptive feedback gain generally preserves more phase/gain reserve.\n');
fprintf('- A small-gain index below one is sufficient for this uncertainty representation.\n');
fprintf('- Sigma-modification and projection exchange exact convergence for boundedness.\n');

function result = runCase(name, gamma, sigma)
    a = 1.0;
    am = 2.0;
    bm = 2.0;
    tau1 = 0.05;
    tau2 = 0.02;
    dt = 5.0e-4;
    finalTime = 20.0;
    steps = floor(finalTime / dt);

    x = 0.0;
    z1 = 0.0;
    z2 = 0.0;
    xm = 0.0;
    k = 0.0;
    kr = 0.0;

    errorEnergy = 0.0;
    controlEnergy = 0.0;
    maxError = 0.0;
    maxK = 0.0;
    maxKr = 0.0;

    for step = 0:(steps - 1)
        t = step * dt;
        r = 1.0 + 0.30 * sin(0.50 * t) + 0.15 * sin(1.70 * t);
        d = 0.05 * sin(8.0 * t);
        e = x - xm;
        u = clampValue(-k * x + kr * r, -30.0, 30.0);

        dx = -a * x + z2 + d;
        dz1 = (-z1 + u) / tau1;
        dz2 = (-z2 + z1) / tau2;
        dxm = -am * xm + bm * r;
        dk = gamma * x * e - sigma * k;
        dkr = -gamma * r * e - sigma * kr;

        x = x + dt * dx;
        z1 = z1 + dt * dz1;
        z2 = z2 + dt * dz2;
        xm = xm + dt * dxm;
        k = clampValue(k + dt * dk, 0.0, 10.0);
        kr = clampValue(kr + dt * dkr, 0.0, 10.0);

        errorEnergy = errorEnergy + e^2 * dt;
        controlEnergy = controlEnergy + u^2 * dt;
        maxError = max(maxError, abs(e));
        maxK = max(maxK, k);
        maxKr = max(maxKr, kr);
    end

    [phaseMarginDeg, gainMarginDb, smallGainIndex] = ...
        frozenMetrics(max(maxK, 1.0e-8), a, tau1, tau2);

    result = struct( ...
        'name', name, ...
        'rmsError', sqrt(errorEnergy / finalTime), ...
        'maxError', maxError, ...
        'rmsControl', sqrt(controlEnergy / finalTime), ...
        'maxK', maxK, ...
        'maxKr', maxKr, ...
        'finalK', k, ...
        'finalKr', kr, ...
        'phaseMarginDeg', phaseMarginDeg, ...
        'gainMarginDb', gainMarginDb, ...
        'smallGainIndex', smallGainIndex);
end

function [phaseMarginDeg, gainMarginDb, maxSmallGain] = frozenMetrics(k, a, tau1, tau2)
    w = logspace(-3, 4, 24000);
    s = 1i * w;
    p0 = 1 ./ (s + a);
    parasitic = 1 ./ ((1 + tau1 * s) .* (1 + tau2 * s));
    plant = p0 .* parasitic;
    loop = k .* plant;

    magDb = 20 * log10(max(abs(loop), realmin));
    phaseDeg = unwrap(angle(loop)) * 180 / pi;

    gainIndex = find(magDb(1:end-1) >= 0 & magDb(2:end) < 0, 1, 'first');
    if isempty(gainIndex)
        phaseMarginDeg = NaN;
    else
        fraction = magDb(gainIndex) / ...
            (magDb(gainIndex) - magDb(gainIndex + 1));
        phaseAtCross = phaseDeg(gainIndex) + fraction * ...
            (phaseDeg(gainIndex + 1) - phaseDeg(gainIndex));
        phaseMarginDeg = 180 + phaseAtCross;
    end

    phaseIndex = find(phaseDeg(1:end-1) > -180 & phaseDeg(2:end) <= -180, 1, 'first');
    if isempty(phaseIndex)
        gainMarginDb = Inf;
    else
        fraction = (phaseDeg(phaseIndex) + 180) / ...
            (phaseDeg(phaseIndex) - phaseDeg(phaseIndex + 1));
        magDbAtCross = magDb(phaseIndex) + fraction * ...
            (magDb(phaseIndex + 1) - magDb(phaseIndex));
        gainMarginDb = -magDbAtCross;
    end

    deltaM = parasitic - 1;
    loop0 = k .* p0;
    t0 = loop0 ./ (1 + loop0);
    maxSmallGain = max(abs(deltaM .* t0));
end

function value = clampValue(value, lower, upper)
    value = max(lower, min(upper, value));
end

function text = formatNumber(value)
    if isnan(value)
        text = 'N/A';
    elseif isinf(value)
        text = 'infinity';
    else
        text = sprintf('%.6f', value);
    end
end

The notebook contains an executable Wolfram Language input cell implementing the same benchmark and returning a Dataset of metrics.

Chapter22_Lesson3.nb

Notebook[{
  Cell["Chapter 22, Lesson 3: Robust-Control Notions and Adaptive Loops", "Title"],
  Cell["Evaluate the input cell to run the reduced-order robust-adaptive benchmark.", "Text"],
  Cell[BoxData["ClearAll[\"Global`*\"];

reference[t_] := 1.0 + 0.30 Sin[0.50 t] + 0.15 Sin[1.70 t];
disturbance[t_] := 0.05 Sin[8.0 t];
clamp[value_, lower_, upper_] := Max[lower, Min[upper, value]];

frozenMetrics[k_, a_, tau1_, tau2_] := Module[
  {w, s, p0, parasitic, plant, loop, magDb, phaseDeg, gainIndex,
   phaseIndex, phaseMargin, gainMarginDb, deltaM, loop0, t0},
  w = Exp@Subdivide[Log[10^-3], Log[10^4], 23999];
  s = I w;
  p0 = 1/(s + a);
  parasitic = 1/((1 + tau1 s) (1 + tau2 s));
  plant = p0 parasitic;
  loop = k plant;
  magDb = 20 Log10[Map[Max[#, 10^-300] &, Abs[loop]]];
  phaseDeg = -(180/Pi) (ArcTan[w/a] + ArcTan[tau1 w] + ArcTan[tau2 w]);

  gainIndex = SelectFirst[
    Range[Length[w] - 1],
    magDb[[#]] >= 0 && magDb[[# + 1]] < 0 &,
    Missing[\"NotFound\"]
  ];
  phaseMargin = If[
    MissingQ[gainIndex],
    Indeterminate,
    With[
      {fraction = magDb[[gainIndex]]/
        (magDb[[gainIndex]] - magDb[[gainIndex + 1]])},
      180 + phaseDeg[[gainIndex]] +
        fraction (phaseDeg[[gainIndex + 1]] - phaseDeg[[gainIndex]])
    ]
  ];

  phaseIndex = SelectFirst[
    Range[Length[w] - 1],
    phaseDeg[[#]] > -180 && phaseDeg[[# + 1]] <= -180 &,
    Missing[\"NotFound\"]
  ];
  gainMarginDb = If[
    MissingQ[phaseIndex],
    Infinity,
    With[
      {fraction = (phaseDeg[[phaseIndex]] + 180)/
        (phaseDeg[[phaseIndex]] - phaseDeg[[phaseIndex + 1]])},
      -(magDb[[phaseIndex]] +
        fraction (magDb[[phaseIndex + 1]] - magDb[[phaseIndex]]))
    ]
  ];

  deltaM = parasitic - 1;
  loop0 = k p0;
  t0 = loop0/(1 + loop0);
  <|
    \"PhaseMarginDeg\" -> phaseMargin,
    \"GainMarginDb\" -> gainMarginDb,
    \"SmallGainIndex\" -> Max[Abs[deltaM t0]]
  |>
];

runCase[name_, gamma_, sigma_] := Module[
  {a = 1.0, am = 2.0, bm = 2.0, tau1 = 0.05, tau2 = 0.02,
   dt = 5.0*10^-4, finalTime = 20.0, steps, x = 0.0, z1 = 0.0,
   z2 = 0.0, xm = 0.0, k = 0.0, kr = 0.0, errorEnergy = 0.0,
   controlEnergy = 0.0, maxError = 0.0, maxK = 0.0, maxKr = 0.0,
   t, r, d, e, u, dx, dz1, dz2, dxm, dk, dkr, metrics},

  steps = Floor[finalTime/dt];

  Do[
    t = step dt;
    r = reference[t];
    d = disturbance[t];
    e = x - xm;
    u = clamp[-k x + kr r, -30.0, 30.0];

    dx = -a x + z2 + d;
    dz1 = (-z1 + u)/tau1;
    dz2 = (-z2 + z1)/tau2;
    dxm = -am xm + bm r;
    dk = gamma x e - sigma k;
    dkr = -gamma r e - sigma kr;

    x += dt dx;
    z1 += dt dz1;
    z2 += dt dz2;
    xm += dt dxm;
    k = clamp[k + dt dk, 0.0, 10.0];
    kr = clamp[kr + dt dkr, 0.0, 10.0];

    errorEnergy += e^2 dt;
    controlEnergy += u^2 dt;
    maxError = Max[maxError, Abs[e]];
    maxK = Max[maxK, k];
    maxKr = Max[maxKr, kr],
    {step, 0, steps - 1}
  ];

  metrics = frozenMetrics[Max[maxK, 10^-8], a, tau1, tau2];

  <|
    \"Case\" -> name,
    \"RMSError\" -> Sqrt[errorEnergy/finalTime],
    \"MaxError\" -> maxError,
    \"RMSControl\" -> Sqrt[controlEnergy/finalTime],
    \"MaxK\" -> maxK,
    \"MaxKr\" -> maxKr,
    \"FinalK\" -> k,
    \"FinalKr\" -> kr,
    \"PhaseMarginDegAtMaxK\" -> metrics[\"PhaseMarginDeg\"],
    \"GainMarginDbAtMaxK\" -> metrics[\"GainMarginDb\"],
    \"SmallGainIndexAtMaxK\" -> metrics[\"SmallGainIndex\"]
  |>
];

results = {
  runCase[\"robustified\", 8.0, 0.15],
  runCase[\"aggressive_plain\", 60.0, 0.0]
};

Dataset[results]
"], "Input"]
},
WindowSize -> {1100, 800},
StyleDefinitions -> "Default.nb"]

With the constants in the supplied programs, the robustified case typically reaches a much smaller peak feedback gain than the aggressive unmodified case. Representative values from the Python implementation are:

These numbers should not be universalized. They depend on the benchmark, initial conditions, projection limits, integration step, disturbance, and reference signal. Their role is to demonstrate a methodology:

1. robust modification limits parameter growth and control effort;
2. the adaptive trajectory then explores a smaller portion of controller space;
3. the smaller controller set retains more frozen frequency-domain reserve;
4. the nonlinear simulation confirms whether that reserve correlates with improved transient behavior.

Note that the aggressive case still has positive frozen margins in this example. This does not contradict its inferior behavior. Classical margins measure one aspect of frozen LTI robustness; they do not quantify all nonlinear transient effects, parameter-rate effects, saturation, or disturbance-driven adaptation.

1. Design a credible baseline loop first. Adaptation should augment, not conceal, a nominal controller with poor damping or inadequate actuator modeling.
2. Separate learnable parameters from dynamic uncertainty. Do not ask a finite-dimensional parameter estimator to reproduce unknown delay or arbitrarily high-order dynamics.
3. Project onto a robustness-compatible set. Choose parameter bounds by checking the entire frozen controller family, not by guessing wide numerical limits.
4. Do not tune adaptation gain from tracking error alone. Include control RMS, parameter-rate RMS, noise sensitivity, crossover movement, and saturation time.
5. Use leakage or a dead zone when irreducible error exists. Persistent disturbances and model mismatch otherwise create parameter drift.
6. Model actuator and sensor dynamics explicitly. Neglected poles, delay, filters, and rate limits often dominate the practical robustness of a fast adaptive law.
7. Validate uncertainty extremes. Sweep parametric corners, unmodeled-mode frequencies, damping ratios, delays, disturbance spectra, and sampling periods.

Problem 1: Multiplicative-Uncertainty Condition

A frozen nominal loop has complementary sensitivity \( T_0(s) \). The plant uncertainty is

\[ P=P_0(1+W_m\Delta_m),\qquad \|\Delta_m\|_\infty\leq 1. \]

Derive a sufficient robust-stability condition and explain its implication for an adaptive controller whose gain increases with time.

Solution:

The uncertain characteristic equation factors as

\[ 1+P_0C+P_0CW_m\Delta_m = (1+L_0)(1+T_0W_m\Delta_m). \]

Assuming nominal internal stability, the small-gain theorem gives

\[ \|T_0W_m\|_\infty<1. \]

If the adaptive gain raises crossover frequency, then \( |T_0| \) can remain large where the uncertainty weight is large. The sufficient robustness inequality may be violated. Therefore, projection bounds and adaptation-gain tuning should be checked against the full frozen controller family.

Problem 2: Frozen Margins Versus Time-Varying Stability

Suppose every frozen matrix \( A(\bar{\theta}) \) is Hurwitz for \( \bar{\theta}\in\Omega \). Is \( \dot{x}=A(\hat{\theta}(t))x \) necessarily uniformly exponentially stable?

Solution:

No. Pointwise Hurwitz stability does not by itself provide a common Lyapunov function or a sufficiently slow parameter variation. Stability can be concluded if stronger conditions hold, such as the existence of one matrix \( P=P^{\mathsf T}>0 \) satisfying

\[ A(\bar{\theta})^{\mathsf T}P+ PA(\bar{\theta})\leq -Q \]

for every \( \bar{\theta}\in\Omega \), with \( Q>0 \). Alternatively, slowly varying system theorems can be used when bounds on \( \dot{\hat{\theta}} \) are sufficiently small. Frozen gain and phase margins alone do not supply either condition.

Problem 3: Ultimate-Bound Radius

Suppose a Lyapunov analysis yields

\[ \dot{V}\leq-\alpha V+\beta,\qquad \alpha>0,\quad\beta\geq 0. \]

Derive a bound on \( V(t) \) and the asymptotic residual level.

Solution:

Multiply by \( e^{\alpha t} \) and integrate:

\[ V(t)\leq e^{-\alpha t}V(0)+ \frac{\beta}{\alpha}\left(1-e^{-\alpha t}\right). \]

Therefore,

\[ \limsup_{t\longrightarrow\infty}V(t)\leq\frac{\beta}{\alpha}. \]

If \( V\geq\lambda_{\min}(P)\|e\|^2 \), then

\[ \limsup_{t\longrightarrow\infty}\|e(t)\| \leq \sqrt{\frac{\beta} {\alpha\lambda_{\min}(P)}}. \]

Leakage typically increases the restoring term in parameter space but also introduces bias through the unknown ideal parameter. The design objective is an acceptable residual bound, not exact parameter convergence.

Problem 4: Exact Uncertainty of Neglected Actuator Poles

The nominal plant is \( P_0(s)=1/(s+a) \), while the true plant is

\[ P(s)= \frac{1} {(s+a)(1+\tau_1s)(1+\tau_2s)}. \]

Express the true plant as \( P=P_0(1+\Delta_m) \) and find the low-frequency approximation of \( \Delta_m \).

Solution:

\[ 1+\Delta_m(s)= \frac{1}{(1+\tau_1s)(1+\tau_2s)}, \]

\[ \Delta_m(s)= \frac{1}{(1+\tau_1s)(1+\tau_2s)}-1. \]

For small \( |s| \), use \( (1+\tau s)^{-1}\approx1-\tau s \):

\[ \Delta_m(s)\approx-(\tau_1+\tau_2)s. \]

Thus relative uncertainty is small at low frequency and grows with frequency, which explains why adaptation that pushes crossover upward can become fragile.

Problem 5: Projection-Set Design

A scalar adaptive feedback gain is projected onto \( 0\leq\hat{k}\leq k_{\max} \). Describe a systematic method for choosing \( k_{\max} \).

Solution:

1. Form \( L_f(s;k)=P(s)C(s;k) \) over a dense grid of \( k\in[0,k_{\max}] \).
2. Check internal stability for every frozen value.
3. Compute the worst-case phase margin, gain margin, \( M_s \), and uncertainty index \( \rho(k) \).
4. Reduce \( k_{\max} \) until all design thresholds are satisfied with reserve.
5. Verify that the resulting set still contains plausible ideal gains for the structured parametric uncertainty.
6. Perform nonlinear simulations including parameter motion, disturbances, actuator limits, delay, noise, and sampling.

Projection bounds are therefore control-design quantities, not merely numerical safeguards.

Robust-control notions provide a valuable diagnostic language for adaptive systems. Additive and multiplicative uncertainty models describe effects that are not captured by the finite-dimensional adaptive parameterization. Frozen-time sensitivity functions and classical margins show how the current controller interacts with plant dynamics, while the small-gain condition \( \|T_0W_m\|_\infty<1 \) quantifies one sufficient robustness reserve.

These tools must not be mistaken for a complete adaptive stability proof. Parameter motion makes the closed loop nonlinear and time varying. Lyapunov analysis, projection, leakage, normalization, parameter-rate bounds, actuator modeling, and nonlinear validation remain essential. Robust adaptive design is therefore a coordinated process: prove boundedness, restrict controller motion to a robustness-compatible set, and verify that the entire frozen family retains adequate frequency-domain reserve.

1. Rohrs, C.E., Valavani, L., Athans, M., & Stein, G. (1985). Robustness of continuous-time adaptive control algorithms in the presence of unmodeled dynamics. IEEE Transactions on Automatic Control, 30(9), 881–889.
2. Ioannou, P.A., & Kokotovic, P.V. (1984). Robust redesign of adaptive control. IEEE Transactions on Automatic Control, 29(3), 202–211.
3. Ioannou, P.A., & Kokotovic, P.V. (1984). Instability analysis and improvement of robustness of adaptive control. Automatica, 20(5), 583–594.
4. Kosut, R.L., & Johnson, C.R., Jr. (1984). An input-output view of robustness in adaptive control. Automatica, 20(5), 569–581.
5. Narendra, K.S., & Annaswamy, A.M. (1986). Robust adaptive control in the presence of bounded disturbances. IEEE Transactions on Automatic Control, 31(4), 306–315.
6. Ioannou, P.A., & Tsakalis, K.S. (1986). A robust direct adaptive controller. IEEE Transactions on Automatic Control, 31(11), 1033–1043.
7. Ortega, R., & Tang, Y. (1989). Robustness of adaptive controllers—A survey. Automatica, 25(5), 651–677.
8. Sastry, S.S. (1984). Model-reference adaptive control—Stability, parameter convergence, and robustness. IMA Journal of Mathematical Control and Information, 1(1), 27–66.