Chapter 24: Implementation Issues and Practical Pitfalls
Lesson 4: Parameter Drift, Reset Strategies, and Monitoring
This lesson studies why adaptive parameters may move far from physically meaningful values even when the tracking error initially appears small. We develop drift indicators, leakage and projection safeguards, hybrid reset logic, hysteresis, dwell-time conditions, and bumpless recovery. The analysis connects continuous-time Lyapunov inequalities to discrete-time monitoring algorithms suitable for real implementations.
1. Learning Objectives and Practical Context
After completing this lesson, students should be able to:
- distinguish parameter adaptation from undesirable parameter drift;
- identify drift caused by bias, noise, unmodeled dynamics, saturation, delay, and loss of excitation;
- derive ultimate-boundedness inequalities for leakage-modified adaptive laws;
- construct parameter-boundary, residual, update-rate, excitation, and saturation monitors;
- compare adaptation freeze, soft reset, hard reset, covariance reset, and fallback control;
- design hysteresis and dwell-time logic that avoids rapid reset cycling;
- implement monitored adaptation in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.
A changing estimate is not automatically a fault. When the true plant changes, a well-designed adaptive controller should change its parameters. The term drift is reserved here for sustained parameter motion caused primarily by noninformative or corrupted adaptation signals, rather than by identifiable plant variation.
flowchart TD
R["Reference r"] --> M["Reference model"]
R --> C["Adaptive controller"]
Y["Measured output"] --> C
C --> SAT["Actuator saturation"]
SAT --> P["Plant"]
P --> Y
Y --> E["Tracking and prediction residuals"]
M --> E
E --> A["Adaptive law"]
Y --> A
A --> TH["Parameter estimate"]
TH --> C
E --> MON["Monitor: residual, \nrate, bounds, \nexcitation, saturation"]
TH --> MON
SAT --> MON
MON --> SUP["Supervisor: continue, \nfreeze, soft reset, \nhard reset, fallback"]
SUP --> A
SUP --> C
2. Baseline MRAC Model and the Meaning of Drift
Consider the first-order plant and reference model
\[ \dot{x}=a_p x+b_p u+d(t), \qquad \dot{x}_m=a_m x_m+b_m r, \]
where \( a_m<0 \), the sign of \( b_p \) is known, and \( d(t) \) represents disturbances, bias-equivalent errors, and the effect of neglected dynamics. Use the direct adaptive control law
\[ u=\hat{k}_x x+\hat{k}_r r =\hat{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}, \qquad \boldsymbol{\phi}= \begin{bmatrix}x\\r\end{bmatrix}. \]
If ideal constant gains exist, they satisfy
\[ a_p+b_p k_x^{*}=a_m,\qquad b_p k_r^{*}=b_m. \]
Define tracking and parameter errors by \( e=x-x_m \) and \( \tilde{\boldsymbol{\theta}} =\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*} \). The error dynamics are
\[ \dot{e} =a_m e+b_p\tilde{\boldsymbol{\theta}}^{T}\boldsymbol{\phi}+d(t). \]
Parameter drift occurs when \( \|\hat{\boldsymbol{\theta}}(t)\| \) grows, approaches an unsafe boundary, or moves persistently without evidence that the motion represents a real and identifiable plant change. Typical causes are:
- measurement bias: a small nonzero mean error is integrated by the update law;
- unmodeled dynamics: adaptation attempts to reproduce dynamics that the chosen parameterization cannot represent;
- weak excitation: many parameter vectors explain the same observed trajectory, so estimates can move in poorly observable directions;
- actuator saturation: the adaptive law assumes a commanded input was applied even when the actuator clipped it;
- sampling, delay, and numerical error: the measured error is no longer synchronized with the regressor used by the update law;
- excessive adaptation gain: noise is converted into rapid parameter motion.
3. Why a Small Persistent Error Can Produce Unbounded Drift
Let \( P>0 \) solve the scalar Lyapunov equation
\[ 2a_mP=-q,\qquad q>0, \quad\Longrightarrow\quad P=\frac{q}{-2a_m}. \]
With the standard gradient law
\[ \dot{\hat{\boldsymbol{\theta}}} =-\boldsymbol{\Gamma}\boldsymbol{\phi}\,eP\, \operatorname{sgn}(b_p), \qquad \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{T}>0, \]
choose
\[ V=Pe^2+|b_p| \tilde{\boldsymbol{\theta}}^{T} \boldsymbol{\Gamma}^{-1} \tilde{\boldsymbol{\theta}}. \]
The cross term cancels and gives
\[ \dot{V}=-qe^2+2Pe\,d(t). \]
In the ideal case \( d(t)=0 \), \( \dot V=-qe^2\leq 0 \). With a persistent disturbance, however, the derivative is not negative semidefinite for all states. More importantly, it contains no strictly negative term in \( \tilde{\boldsymbol{\theta}} \). Therefore a small biased adaptation signal may accumulate over a long period.
A scalar illustration makes the mechanism explicit. Suppose the regressor and error settle temporarily to nonzero constants \( \phi_0 \) and \( e_0 \). Then
\[ \dot{\hat{\theta}}=-\gamma\phi_0e_0, \qquad \hat{\theta}(t) =\hat{\theta}(0)-\gamma\phi_0e_0t. \]
Even a very small product \( \phi_0e_0 \) causes linear-in-time drift. This is why monitoring only the instantaneous tracking error is insufficient.
4. Leakage, Dead Zones, Normalization, and Projection
4.1 Leakage or sigma modification
Introduce a validated center \( \boldsymbol{\theta}_0 \) and use
\[ \dot{\hat{\boldsymbol{\theta}}} =-\boldsymbol{\Gamma}\boldsymbol{\phi}\,eP\, \operatorname{sgn}(b_p) -\sigma\boldsymbol{\Gamma} \left(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_0\right), \qquad \sigma>0. \]
Since \( \hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_0 =\tilde{\boldsymbol{\theta}} +\boldsymbol{\theta}^{*}-\boldsymbol{\theta}_0 \), Young's inequality gives
\[ -2\sigma|b_p| \tilde{\boldsymbol{\theta}}^{T} \left(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_0\right) \leq -\sigma|b_p|\|\tilde{\boldsymbol{\theta}}\|^2 +\sigma|b_p| \|\boldsymbol{\theta}^{*}-\boldsymbol{\theta}_0\|^2. \]
For \( |d(t)|\leq\bar d \),
\[ 2P|e|\,|d| \leq \frac{q}{2}e^2+\frac{2P^2}{q}\bar d^2. \]
Therefore
\[ \dot V \leq -\frac{q}{2}e^2 -\sigma|b_p|\|\tilde{\boldsymbol{\theta}}\|^2 +C, \]
\[ C= \frac{2P^2}{q}\bar d^2 +\sigma|b_p| \|\boldsymbol{\theta}^{*}-\boldsymbol{\theta}_0\|^2. \]
Let \( \alpha=\min\!\left\{ q/(2P),\,\sigma\lambda_{\min}(\boldsymbol{\Gamma}) \right\} \). Then
\[ \dot V\leq-\alpha V+C \quad\Longrightarrow\quad V(t)\leq e^{-\alpha t}V(0) +\frac{C}{\alpha}\left(1-e^{-\alpha t}\right). \]
Thus the error and parameter error are uniformly ultimately bounded. Leakage trades exact asymptotic parameter recovery for robustness.
4.2 Dead-zone adaptation
A dead zone prevents small noise-dominated errors from driving the estimator:
\[ \dot{\hat{\boldsymbol{\theta}}} = \begin{cases} \mathbf{0}, & |s(t)|\leq\delta,\\ -\boldsymbol{\Gamma}\boldsymbol{\phi}s(t), & |s(t)|>\delta, \end{cases} \]
where \( s(t) \) is the adaptation error signal and \( \delta \) is selected above the expected noise and modeling-error floor.
4.3 Normalization
To prevent large regressors from creating excessively large parameter increments, use
\[ \dot{\hat{\boldsymbol{\theta}}} = -\frac{\boldsymbol{\Gamma}\boldsymbol{\phi}s} {\delta_n+\boldsymbol{\phi}^{T}\boldsymbol{\phi}}, \qquad \delta_n>0. \]
4.4 Projection
If engineering knowledge provides a compact admissible set \( \Omega \), use
\[ \dot{\hat{\boldsymbol{\theta}}} =\operatorname{Proj}_{\Omega} \left(\hat{\boldsymbol{\theta}},\mathbf{y}\right), \]
where a Lyapunov-compatible projection satisfies, for \( \boldsymbol{\theta}^{*}\in\Omega \),
\[ \tilde{\boldsymbol{\theta}}^{T} \boldsymbol{\Gamma}^{-1} \left[ \operatorname{Proj}_{\Omega} \left(\hat{\boldsymbol{\theta}},\mathbf{y}\right) -\mathbf{y} \right] \leq 0. \]
Simple componentwise clipping is useful in software, but a formal projection operator is preferred when preservation of a Lyapunov proof is required.
5. Monitoring Signals for Drift Detection
No single monitor is reliable in all operating conditions. A practical supervisor combines several indicators over finite windows.
5.1 Exponentially weighted tracking residual
\[ E_k=\beta_eE_{k-1}+(1-\beta_e)e_k^2, \qquad e_{\mathrm{rms},k}=\sqrt{E_k}, \qquad 0<\beta_e<1. \]
5.2 Parameter update-rate monitor
\[ R_k=\beta_rR_{k-1} +(1-\beta_r) \left\| \frac{\hat{\boldsymbol{\theta}}_k -\hat{\boldsymbol{\theta}}_{k-1}}{T_s} \right\|^2. \]
5.3 Boundary-proximity index
For componentwise bounds \( \theta_{i,\min}\leq\hat\theta_i \leq\theta_{i,\max} \), define center and half-span
\[ c_i=\frac{\theta_{i,\max}+\theta_{i,\min}}{2}, \qquad h_i=\frac{\theta_{i,\max}-\theta_{i,\min}}{2}, \]
\[ B_{\theta,k} = \sqrt{ \frac{1}{p} \sum_{i=1}^{p} \left( \frac{\hat\theta_{i,k}-c_i}{h_i} \right)^2 }. \]
A value near one means that the estimate is approaching the boundary in a normalized multi-parameter sense.
5.4 Finite-window excitation monitor
Use the exponentially weighted information matrix
\[ \mathbf{G}_k =\beta_g\mathbf{G}_{k-1} +(1-\beta_g) \boldsymbol{\phi}_k\boldsymbol{\phi}_k^{T}. \]
The quantity \( \lambda_{\min}(\mathbf{G}_k) \) indicates whether all parameter directions are being excited over the selected memory horizon. A small eigenvalue does not prove drift, but it warns that parameter interpretation is unreliable.
5.5 Saturation and implementation-condition monitors
\[ S_k = \beta_s S_{k-1} + (1-\beta_s) \mathbf{1}_{\left\{ \left|u_{\mathrm{cmd},k}-u_{\mathrm{act},k}\right| > \varepsilon_u \right\}}. \]
Similar binary or continuous monitors should cover timing overruns, invalid sensor data, excessive delay, arithmetic overflow, and failed plausibility checks.
5.6 Composite score
A dimensionless engineering score can be formed as
\[ J_k= w_e\frac{\sqrt{E_k}}{E_0} +w_r\frac{\sqrt{R_k}}{R_0} +w_bB_{\theta,k} +w_sS_k +w_g \max\left( 0, \frac{\lambda_0-\lambda_{\min}(\mathbf{G}_k)} {\lambda_0} \right), \]
with nonnegative weights. This score is a supervisory diagnostic, not a Lyapunov function. Thresholds should be derived from validated nominal tests and worst-case simulations, not selected only to make one demonstration look successful.
6. Reset Strategies and Their Mathematical Effects
6.1 Adaptation freeze
Set \( \dot{\hat{\boldsymbol{\theta}}}=\mathbf{0} \) while retaining the current controller. This is minimally disruptive, but it does not repair a parameter vector that is already unsafe.
6.2 Soft or partial reset
\[ \hat{\boldsymbol{\theta}}^{+} =(1-\rho)\hat{\boldsymbol{\theta}}^{-} +\rho\boldsymbol{\theta}_s, \qquad 0<\rho<1, \]
where \( \boldsymbol{\theta}_s \) is a validated safe estimate. With \( \tilde{\boldsymbol{\theta}}_s =\boldsymbol{\theta}_s-\boldsymbol{\theta}^{*} \),
\[ \tilde{\boldsymbol{\theta}}^{+} =(1-\rho)\tilde{\boldsymbol{\theta}}^{-} +\rho\tilde{\boldsymbol{\theta}}_s. \]
For any symmetric \( \mathbf{W}>0 \),
\[ \left\| (1-\rho)\mathbf{a}+\rho\mathbf{b} \right\|_{\mathbf{W}}^2 = (1-\rho)\|\mathbf{a}\|_{\mathbf{W}}^2 +\rho\|\mathbf{b}\|_{\mathbf{W}}^2 -\rho(1-\rho) \|\mathbf{a}-\mathbf{b}\|_{\mathbf{W}}^2. \]
Hence a soft reset reduces parameter error when the safe estimate is sufficiently closer to the ideal parameter than the current estimate. Since the ideal parameter is unknown, this is an engineering assumption that must be supported by validation.
6.3 Hard reset
\[ \hat{\boldsymbol{\theta}}^{+} =\boldsymbol{\theta}_s. \]
Hard reset gives the fastest return to a known configuration but may produce a control discontinuity. The controller output should therefore be rate-limited, blended, or initialized for bumpless transfer.
6.4 RLS covariance reset
For recursive least squares, the covariance matrix controls adaptation aggressiveness. A covariance reset may use
\[ \mathbf{P}^{+} =(1-\alpha)\mathbf{P}^{-} +\alpha\mathbf{P}_0, \qquad 0\leq\alpha\leq1, \]
where \( \mathbf{P}_0>0 \). Convexity guarantees \( \mathbf{P}^{+}>0 \). Resetting covariance without resetting the parameter estimate restores responsiveness while avoiding an immediate controller-gain jump.
6.5 Fallback controller
A validated fixed-gain controller can be selected when monitors indicate that the assumptions required for adaptation are not credible. Switching must respect the stability and actuator constraints of both controllers.
7. Hybrid Supervisor: Hysteresis, Persistence, and Dwell Time
A threshold crossing caused by one noisy sample should not trigger a reset. Let \( J_{\mathrm{on}}>J_{\mathrm{off}} \). An alarm is declared only when \( J_k>J_{\mathrm{on}} \) persists for \( T_p \). It is cleared only after \( J_k<J_{\mathrm{off}} \). A minimum dwell time \( T_d \) between resets prevents rapid cycling.
\[ N_p=\left\lceil\frac{T_p}{T_s}\right\rceil, \qquad N_d=\left\lceil\frac{T_d}{T_s}\right\rceil. \]
flowchart TD
N["NORMAL: adapt and monitor"] --> C{"Score above ON threshold \nfor persistence time?"}
C -->|"no"| N
C -->|"yes and dwell time satisfied"| R["RESET: soft or hard parameter action"]
R --> H["HOLD: freeze adaptation and blend control"]
H --> V{"Signals valid and score \nbelow OFF threshold?"}
V -->|"no"| F["FALLBACK: fixed validated controller"]
V -->|"yes"| N
F --> Q{"Recovery test passed and \ndwell time satisfied?"}
Q -->|"no"| F
Q -->|"yes"| H
The reset event is a jump in a hybrid system. A complete stability argument must analyze both continuous flows and jumps:
\[ \dot V\leq-\alpha V+C \quad\text{during flows}, \qquad V^{+}-V^{-}\leq\Delta_J \quad\text{at jumps}. \]
If \( \Delta_J\leq0 \), resets do not increase the Lyapunov function. If a bounded positive jump is possible, dwell-time and recovery conditions must ensure that flow intervals dissipate more energy than jumps add.
8. Tuning and Validation Procedure
- Establish physically justified parameter bounds and a validated safe parameter vector.
- Run nominal tests that cover expected references, loads, noise, sampling jitter, and actuator use.
- Select monitor time constants longer than dominant noise periods but shorter than the maximum acceptable drift-detection delay.
- Normalize every metric so that score weights have interpretable relative influence.
- Test bias, dropout, delay, saturation, weak excitation, plant change, and unmodeled high-frequency dynamics separately and jointly.
- Verify that a real plant change is not mistaken for drift before the adaptive loop has sufficient time and excitation to respond.
- Validate bumpless transfer, control-rate limits, reset counters, event logs, and fallback recovery.
- Store the reason, thresholds, parameter vector before and after the reset, and relevant signal windows for post-test analysis.
Monitor thresholds are part of the control-system design. They should be configuration-controlled and reviewed with the same discipline as controller gains and actuator limits.
9. Python Implementation
The implementation uses only the Python standard library. For larger systems, NumPy provides vectorized linear algebra, SciPy provides numerical tools, and the Python Control Systems Library can support plant and reference-model analysis. The algorithm below remains explicit so every monitoring and reset step is visible.
Chapter24_Lesson4.py
# Chapter24_Lesson4.py
"""
Monitored MRAC demonstration for:
Chapter 24, Lesson 4 - Parameter Drift, Reset Strategies, and Monitoring.
The example uses a first-order plant, a reference model, normalized gradient
adaptation, leakage, projection, exponentially weighted monitoring metrics,
hysteresis, dwell time, soft parameter resets, and temporary adaptation freeze.
No third-party packages are required. The script writes a CSV log that can be
plotted with a spreadsheet, pandas, MATLAB, or Mathematica.
"""
from __future__ import annotations
import csv
import math
import random
from dataclasses import dataclass
from pathlib import Path
@dataclass
class MonitorState:
error_sq: float = 0.0
rate_sq: float = 0.0
saturation: float = 0.0
g00: float = 0.0
g01: float = 0.0
g11: float = 0.0
exceed_count: int = 0
alarm: bool = False
freeze_count: int = 0
last_reset_time: float = -1.0e9
def clamp(value: float, lower: float, upper: float) -> float:
return max(lower, min(upper, value))
def min_eigenvalue_2x2(g00: float, g01: float, g11: float) -> float:
discriminant = math.sqrt(max(0.0, (g00 - g11) ** 2 + 4.0 * g01 * g01))
return max(0.0, 0.5 * (g00 + g11 - discriminant))
def run_simulation(output_csv: Path = Path("Chapter24_Lesson4_results.csv")) -> list[dict[str, float]]:
dt = 0.01
final_time = 120.0
steps = int(final_time / dt)
# Plant and reference-model parameters.
a_plant_nominal = -0.8
a_plant_changed = -1.05
b_plant = 1.1
a_model = -1.5
b_model = 1.5
# Ideal gains for the nominal plant:
# a_p + b_p*k_x = a_m and b_p*k_r = b_m.
theta_safe = [
(a_model - a_plant_nominal) / b_plant,
b_model / b_plant,
]
theta = theta_safe.copy()
lower = [-3.0, 0.2]
upper = [1.0, 3.0]
gamma = [10.0, 8.0]
leakage = 0.002
dead_zone = 0.005
control_limit = 4.0
# Monitoring and reset settings.
score_on = 0.65
score_off = 0.36
persistence_time = 0.70
dwell_time = 8.0
freeze_time = 2.0
reset_fraction = 0.55
beta_fast = math.exp(-dt / 1.0)
beta_gramian = math.exp(-dt / 3.0)
x = 0.0
x_model = 0.0
monitor = MonitorState()
rows: list[dict[str, float]] = []
reset_events: list[tuple[float, float]] = []
random.seed(4)
for k in range(steps):
time = k * dt
# A genuine plant change occurs after 75 s.
a_plant = a_plant_nominal if time < 75.0 else a_plant_changed
# Rich reference first, then a weakly exciting command after 85 s.
if time < 85.0:
reference = math.sin(0.45 * time) + 0.4 * math.sin(1.2 * time)
else:
reference = 0.15 * math.sin(0.2 * time)
# A deliberately large sensor bias demonstrates drift and reset behavior.
sensor_bias = 0.50 if 35.0 <= time < 65.0 else 0.0
measurement_noise = 0.025 * random.gauss(0.0, 1.0)
x_model += dt * (a_model * x_model + b_model * reference)
measured_x = x + sensor_bias + measurement_noise
unsaturated_control = theta[0] * measured_x + theta[1] * reference
control = clamp(unsaturated_control, -control_limit, control_limit)
plant_disturbance = 0.12 * math.sin(2.3 * time)
x += dt * (a_plant * x + b_plant * control + plant_disturbance)
tracking_error = measured_x - x_model
phi0 = measured_x
phi1 = reference
normalization = 1.0 + phi0 * phi0 + phi1 * phi1
if monitor.freeze_count > 0:
theta_dot = [0.0, 0.0]
monitor.freeze_count -= 1
else:
theta_dot = [
-gamma[0] * phi0 * tracking_error / normalization
- leakage * (theta[0] - theta_safe[0]),
-gamma[1] * phi1 * tracking_error / normalization
- leakage * (theta[1] - theta_safe[1]),
]
if abs(tracking_error) <= dead_zone:
theta_dot = [
-leakage * (theta[0] - theta_safe[0]),
-leakage * (theta[1] - theta_safe[1]),
]
# Euler update followed by box projection.
for i in range(2):
theta[i] = clamp(theta[i] + dt * theta_dot[i], lower[i], upper[i])
# Exponentially weighted monitoring signals.
monitor.error_sq = (
beta_fast * monitor.error_sq
+ (1.0 - beta_fast) * tracking_error * tracking_error
)
monitor.rate_sq = (
beta_fast * monitor.rate_sq
+ (1.0 - beta_fast) * (theta_dot[0] ** 2 + theta_dot[1] ** 2)
)
saturated = 1.0 if abs(unsaturated_control - control) > 1.0e-12 else 0.0
monitor.saturation = (
beta_fast * monitor.saturation + (1.0 - beta_fast) * saturated
)
monitor.g00 = beta_gramian * monitor.g00 + (1.0 - beta_gramian) * phi0 * phi0
monitor.g01 = beta_gramian * monitor.g01 + (1.0 - beta_gramian) * phi0 * phi1
monitor.g11 = beta_gramian * monitor.g11 + (1.0 - beta_gramian) * phi1 * phi1
lambda_min = min_eigenvalue_2x2(
monitor.g00, monitor.g01, monitor.g11
)
center = [(lower[i] + upper[i]) / 2.0 for i in range(2)]
half_span = [(upper[i] - lower[i]) / 2.0 for i in range(2)]
parameter_ratio = math.sqrt(
sum(
((theta[i] - center[i]) / half_span[i]) ** 2
for i in range(2)
)
/ 2.0
)
error_metric = math.sqrt(monitor.error_sq) / 0.20
rate_metric = math.sqrt(monitor.rate_sq) / 1.50
excitation_deficit = max(0.0, (0.04 - lambda_min) / 0.04)
score = (
0.40 * error_metric
+ 0.25 * parameter_ratio
+ 0.20 * rate_metric
+ 0.05 * monitor.saturation
+ 0.10 * excitation_deficit
)
# A near-bound estimate is treated as a high-severity condition.
if parameter_ratio > 0.90:
score = max(score, 1.10)
if score > score_on:
monitor.exceed_count += 1
elif score < 0.75 * score_on:
monitor.exceed_count = max(0, monitor.exceed_count - 2)
persistent_alarm = monitor.exceed_count * dt >= persistence_time
dwell_satisfied = time - monitor.last_reset_time >= dwell_time
if (not monitor.alarm) and persistent_alarm and dwell_satisfied:
monitor.alarm = True
# Soft reset: move only partway toward the validated safe estimate.
for i in range(2):
theta[i] = (
(1.0 - reset_fraction) * theta[i]
+ reset_fraction * theta_safe[i]
)
monitor.freeze_count = int(freeze_time / dt)
monitor.last_reset_time = time
monitor.exceed_count = 0
reset_events.append((time, score))
print(
f"Reset at t={time:6.2f} s, score={score:.3f}, "
f"theta=[{theta[0]:.4f}, {theta[1]:.4f}]"
)
if monitor.alarm and score < score_off:
monitor.alarm = False
rows.append(
{
"time": time,
"x": x,
"x_model": x_model,
"reference": reference,
"tracking_error": tracking_error,
"control": control,
"kx_hat": theta[0],
"kr_hat": theta[1],
"monitor_score": score,
"error_rms": math.sqrt(monitor.error_sq),
"parameter_ratio": parameter_ratio,
"lambda_min_gramian": lambda_min,
"saturation_ewma": monitor.saturation,
"alarm": 1.0 if monitor.alarm else 0.0,
"frozen": 1.0 if monitor.freeze_count > 0 else 0.0,
}
)
with output_csv.open("w", newline="", encoding="utf-8") as csv_file:
writer = csv.DictWriter(csv_file, fieldnames=list(rows[0].keys()))
writer.writeheader()
writer.writerows(rows)
print(f"Wrote {len(rows)} samples to {output_csv.resolve()}")
print(f"Number of reset events: {len(reset_events)}")
print(f"Final estimate: kx={theta[0]:.5f}, kr={theta[1]:.5f}")
return rows
if __name__ == "__main__":
run_simulation()
10. C++ Implementation
The C++17 version uses the standard library. In higher-dimensional controllers, Eigen is a common choice for matrices, eigenvalues, and factorizations. Real-time deployments should replace dynamic logging and console output with bounded, nonblocking telemetry.
Chapter24_Lesson4.cpp
// Chapter24_Lesson4.cpp
// Monitored first-order MRAC with projection, drift monitoring, soft reset,
// hysteresis, dwell time, and temporary adaptation freeze.
//
// Build:
// g++ -std=c++17 -O2 -Wall -Wextra -pedantic Chapter24_Lesson4.cpp -o Chapter24_Lesson4
// Run:
// ./Chapter24_Lesson4
#include <algorithm>
#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <random>
#include <stdexcept>
#include <string>
#include <vector>
struct MonitorState {
double error_sq = 0.0;
double rate_sq = 0.0;
double saturation = 0.0;
double g00 = 0.0;
double g01 = 0.0;
double g11 = 0.0;
int exceed_count = 0;
bool alarm = false;
int freeze_count = 0;
double last_reset_time = -1.0e9;
};
double clamp_value(double value, double lower, double upper) {
return std::max(lower, std::min(upper, value));
}
double min_eigenvalue_2x2(double g00, double g01, double g11) {
const double discriminant =
std::sqrt(std::max(0.0, (g00 - g11) * (g00 - g11) + 4.0 * g01 * g01));
return std::max(0.0, 0.5 * (g00 + g11 - discriminant));
}
int main() {
try {
constexpr double dt = 0.01;
constexpr double final_time = 120.0;
const int steps = static_cast<int>(final_time / dt);
constexpr double a_plant_nominal = -0.8;
constexpr double a_plant_changed = -1.05;
constexpr double b_plant = 1.1;
constexpr double a_model = -1.5;
constexpr double b_model = 1.5;
const std::array<double, 2> theta_safe{
(a_model - a_plant_nominal) / b_plant,
b_model / b_plant
};
std::array<double, 2> theta = theta_safe;
const std::array<double, 2> lower{-3.0, 0.2};
const std::array<double, 2> upper{1.0, 3.0};
const std::array<double, 2> gamma{10.0, 8.0};
constexpr double leakage = 0.002;
constexpr double dead_zone = 0.005;
constexpr double control_limit = 4.0;
constexpr double score_on = 0.65;
constexpr double score_off = 0.36;
constexpr double persistence_time = 0.70;
constexpr double dwell_time = 8.0;
constexpr double freeze_time = 2.0;
constexpr double reset_fraction = 0.55;
const double beta_fast = std::exp(-dt / 1.0);
const double beta_gramian = std::exp(-dt / 3.0);
double x = 0.0;
double x_model = 0.0;
MonitorState monitor;
int reset_count = 0;
std::mt19937 generator(4);
std::normal_distribution<double> normal(0.0, 1.0);
std::ofstream csv("Chapter24_Lesson4_results.csv");
if (!csv) {
throw std::runtime_error("Cannot open Chapter24_Lesson4_results.csv");
}
csv << "time,x,x_model,reference,tracking_error,control,"
"kx_hat,kr_hat,monitor_score,error_rms,parameter_ratio,"
"lambda_min_gramian,saturation_ewma,alarm,frozen\n";
csv << std::setprecision(12);
for (int k = 0; k < steps; ++k) {
const double time = k * dt;
const double a_plant =
(time < 75.0) ? a_plant_nominal : a_plant_changed;
const double reference =
(time < 85.0)
? std::sin(0.45 * time) + 0.4 * std::sin(1.2 * time)
: 0.15 * std::sin(0.2 * time);
const double sensor_bias =
(time >= 35.0 && time < 65.0) ? 0.50 : 0.0;
const double measurement_noise = 0.025 * normal(generator);
x_model += dt * (a_model * x_model + b_model * reference);
const double measured_x = x + sensor_bias + measurement_noise;
const double unsaturated_control =
theta[0] * measured_x + theta[1] * reference;
const double control =
clamp_value(unsaturated_control, -control_limit, control_limit);
const double plant_disturbance = 0.12 * std::sin(2.3 * time);
x += dt * (a_plant * x + b_plant * control + plant_disturbance);
const double tracking_error = measured_x - x_model;
const double phi0 = measured_x;
const double phi1 = reference;
const double normalization = 1.0 + phi0 * phi0 + phi1 * phi1;
std::array<double, 2> theta_dot{0.0, 0.0};
if (monitor.freeze_count > 0) {
--monitor.freeze_count;
} else {
theta_dot[0] =
-gamma[0] * phi0 * tracking_error / normalization
- leakage * (theta[0] - theta_safe[0]);
theta_dot[1] =
-gamma[1] * phi1 * tracking_error / normalization
- leakage * (theta[1] - theta_safe[1]);
if (std::abs(tracking_error) <= dead_zone) {
theta_dot[0] = -leakage * (theta[0] - theta_safe[0]);
theta_dot[1] = -leakage * (theta[1] - theta_safe[1]);
}
}
for (std::size_t i = 0; i < theta.size(); ++i) {
theta[i] = clamp_value(
theta[i] + dt * theta_dot[i], lower[i], upper[i]);
}
monitor.error_sq =
beta_fast * monitor.error_sq
+ (1.0 - beta_fast) * tracking_error * tracking_error;
monitor.rate_sq =
beta_fast * monitor.rate_sq
+ (1.0 - beta_fast)
* (theta_dot[0] * theta_dot[0]
+ theta_dot[1] * theta_dot[1]);
const double saturated =
(std::abs(unsaturated_control - control) > 1.0e-12) ? 1.0 : 0.0;
monitor.saturation =
beta_fast * monitor.saturation + (1.0 - beta_fast) * saturated;
monitor.g00 =
beta_gramian * monitor.g00 + (1.0 - beta_gramian) * phi0 * phi0;
monitor.g01 =
beta_gramian * monitor.g01 + (1.0 - beta_gramian) * phi0 * phi1;
monitor.g11 =
beta_gramian * monitor.g11 + (1.0 - beta_gramian) * phi1 * phi1;
const double lambda_min =
min_eigenvalue_2x2(monitor.g00, monitor.g01, monitor.g11);
double parameter_ratio_sq = 0.0;
for (std::size_t i = 0; i < theta.size(); ++i) {
const double center = 0.5 * (lower[i] + upper[i]);
const double half_span = 0.5 * (upper[i] - lower[i]);
const double normalized = (theta[i] - center) / half_span;
parameter_ratio_sq += normalized * normalized;
}
const double parameter_ratio =
std::sqrt(parameter_ratio_sq / 2.0);
const double error_metric = std::sqrt(monitor.error_sq) / 0.20;
const double rate_metric = std::sqrt(monitor.rate_sq) / 1.50;
const double excitation_deficit =
std::max(0.0, (0.04 - lambda_min) / 0.04);
double score =
0.40 * error_metric
+ 0.25 * parameter_ratio
+ 0.20 * rate_metric
+ 0.05 * monitor.saturation
+ 0.10 * excitation_deficit;
if (parameter_ratio > 0.90) {
score = std::max(score, 1.10);
}
if (score > score_on) {
++monitor.exceed_count;
} else if (score < 0.75 * score_on) {
monitor.exceed_count = std::max(0, monitor.exceed_count - 2);
}
const bool persistent_alarm =
monitor.exceed_count * dt >= persistence_time;
const bool dwell_satisfied =
time - monitor.last_reset_time >= dwell_time;
if (!monitor.alarm && persistent_alarm && dwell_satisfied) {
monitor.alarm = true;
for (std::size_t i = 0; i < theta.size(); ++i) {
theta[i] =
(1.0 - reset_fraction) * theta[i]
+ reset_fraction * theta_safe[i];
}
monitor.freeze_count = static_cast<int>(freeze_time / dt);
monitor.last_reset_time = time;
monitor.exceed_count = 0;
++reset_count;
std::cout << std::fixed << std::setprecision(3)
<< "Reset at t=" << time
<< " s, score=" << score
<< ", theta=[" << theta[0]
<< ", " << theta[1] << "]\n";
}
if (monitor.alarm && score < score_off) {
monitor.alarm = false;
}
csv << time << ','
<< x << ','
<< x_model << ','
<< reference << ','
<< tracking_error << ','
<< control << ','
<< theta[0] << ','
<< theta[1] << ','
<< score << ','
<< std::sqrt(monitor.error_sq) << ','
<< parameter_ratio << ','
<< lambda_min << ','
<< monitor.saturation << ','
<< (monitor.alarm ? 1 : 0) << ','
<< (monitor.freeze_count > 0 ? 1 : 0) << '\n';
}
std::cout << "Wrote " << steps
<< " samples to Chapter24_Lesson4_results.csv\n";
std::cout << "Number of reset events: " << reset_count << '\n';
std::cout << std::setprecision(6)
<< "Final estimate: kx=" << theta[0]
<< ", kr=" << theta[1] << '\n';
return 0;
} catch (const std::exception& exception) {
std::cerr << "Error: " << exception.what() << '\n';
return 1;
}
}
11. Java Implementation
The Java version uses JDK classes only. EJML or Apache Commons Math can be used when the parameter vector and information matrix require general matrix operations.
Chapter24_Lesson4.java
// Chapter24_Lesson4.java
// Monitored first-order MRAC with projection, drift monitoring, soft reset,
// hysteresis, dwell time, and temporary adaptation freeze.
//
// Build:
// javac Chapter24_Lesson4.java
// Run:
// java Chapter24_Lesson4
import java.io.BufferedWriter;
import java.io.IOException;
import java.nio.charset.StandardCharsets;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.Locale;
import java.util.Random;
public final class Chapter24_Lesson4 {
private Chapter24_Lesson4() {
}
private static final class MonitorState {
double errorSq;
double rateSq;
double saturation;
double g00;
double g01;
double g11;
int exceedCount;
boolean alarm;
int freezeCount;
double lastResetTime = -1.0e9;
}
private static double clamp(double value, double lower, double upper) {
return Math.max(lower, Math.min(upper, value));
}
private static double minEigenvalue2x2(double g00, double g01, double g11) {
double discriminant = Math.sqrt(
Math.max(0.0, Math.pow(g00 - g11, 2.0) + 4.0 * g01 * g01)
);
return Math.max(0.0, 0.5 * (g00 + g11 - discriminant));
}
public static void main(String[] args) {
Locale.setDefault(Locale.US);
final double dt = 0.01;
final double finalTime = 120.0;
final int steps = (int) (finalTime / dt);
final double aPlantNominal = -0.8;
final double aPlantChanged = -1.05;
final double bPlant = 1.1;
final double aModel = -1.5;
final double bModel = 1.5;
final double[] thetaSafe = {
(aModel - aPlantNominal) / bPlant,
bModel / bPlant
};
final double[] theta = thetaSafe.clone();
final double[] lower = {-3.0, 0.2};
final double[] upper = {1.0, 3.0};
final double[] gamma = {10.0, 8.0};
final double leakage = 0.002;
final double deadZone = 0.005;
final double controlLimit = 4.0;
final double scoreOn = 0.65;
final double scoreOff = 0.36;
final double persistenceTime = 0.70;
final double dwellTime = 8.0;
final double freezeTime = 2.0;
final double resetFraction = 0.55;
final double betaFast = Math.exp(-dt / 1.0);
final double betaGramian = Math.exp(-dt / 3.0);
double x = 0.0;
double xModel = 0.0;
int resetCount = 0;
MonitorState monitor = new MonitorState();
Random random = new Random(4L);
Path output = Path.of("Chapter24_Lesson4_results.csv");
try (BufferedWriter writer = Files.newBufferedWriter(
output, StandardCharsets.UTF_8)) {
writer.write(
"time,x,x_model,reference,tracking_error,control,"
+ "kx_hat,kr_hat,monitor_score,error_rms,parameter_ratio,"
+ "lambda_min_gramian,saturation_ewma,alarm,frozen"
);
writer.newLine();
for (int k = 0; k < steps; ++k) {
double time = k * dt;
double aPlant = time < 75.0
? aPlantNominal
: aPlantChanged;
double reference = time < 85.0
? Math.sin(0.45 * time) + 0.4 * Math.sin(1.2 * time)
: 0.15 * Math.sin(0.2 * time);
double sensorBias =
(time >= 35.0 && time < 65.0) ? 0.50 : 0.0;
double measurementNoise = 0.025 * random.nextGaussian();
xModel += dt * (aModel * xModel + bModel * reference);
double measuredX = x + sensorBias + measurementNoise;
double unsaturatedControl =
theta[0] * measuredX + theta[1] * reference;
double control = clamp(
unsaturatedControl, -controlLimit, controlLimit
);
double plantDisturbance = 0.12 * Math.sin(2.3 * time);
x += dt * (aPlant * x + bPlant * control + plantDisturbance);
double trackingError = measuredX - xModel;
double phi0 = measuredX;
double phi1 = reference;
double normalization =
1.0 + phi0 * phi0 + phi1 * phi1;
double[] thetaDot = {0.0, 0.0};
if (monitor.freezeCount > 0) {
monitor.freezeCount -= 1;
} else {
thetaDot[0] =
-gamma[0] * phi0 * trackingError / normalization
- leakage * (theta[0] - thetaSafe[0]);
thetaDot[1] =
-gamma[1] * phi1 * trackingError / normalization
- leakage * (theta[1] - thetaSafe[1]);
if (Math.abs(trackingError) <= deadZone) {
thetaDot[0] =
-leakage * (theta[0] - thetaSafe[0]);
thetaDot[1] =
-leakage * (theta[1] - thetaSafe[1]);
}
}
for (int i = 0; i < theta.length; ++i) {
theta[i] = clamp(
theta[i] + dt * thetaDot[i],
lower[i],
upper[i]
);
}
monitor.errorSq =
betaFast * monitor.errorSq
+ (1.0 - betaFast) * trackingError * trackingError;
monitor.rateSq =
betaFast * monitor.rateSq
+ (1.0 - betaFast)
* (thetaDot[0] * thetaDot[0]
+ thetaDot[1] * thetaDot[1]);
double saturated =
Math.abs(unsaturatedControl - control) > 1.0e-12
? 1.0
: 0.0;
monitor.saturation =
betaFast * monitor.saturation
+ (1.0 - betaFast) * saturated;
monitor.g00 =
betaGramian * monitor.g00
+ (1.0 - betaGramian) * phi0 * phi0;
monitor.g01 =
betaGramian * monitor.g01
+ (1.0 - betaGramian) * phi0 * phi1;
monitor.g11 =
betaGramian * monitor.g11
+ (1.0 - betaGramian) * phi1 * phi1;
double lambdaMin = minEigenvalue2x2(
monitor.g00, monitor.g01, monitor.g11
);
double parameterRatioSq = 0.0;
for (int i = 0; i < theta.length; ++i) {
double center = 0.5 * (lower[i] + upper[i]);
double halfSpan = 0.5 * (upper[i] - lower[i]);
double normalized =
(theta[i] - center) / halfSpan;
parameterRatioSq += normalized * normalized;
}
double parameterRatio =
Math.sqrt(parameterRatioSq / 2.0);
double errorMetric =
Math.sqrt(monitor.errorSq) / 0.20;
double rateMetric =
Math.sqrt(monitor.rateSq) / 1.50;
double excitationDeficit =
Math.max(0.0, (0.04 - lambdaMin) / 0.04);
double score =
0.40 * errorMetric
+ 0.25 * parameterRatio
+ 0.20 * rateMetric
+ 0.05 * monitor.saturation
+ 0.10 * excitationDeficit;
if (parameterRatio > 0.90) {
score = Math.max(score, 1.10);
}
if (score > scoreOn) {
monitor.exceedCount += 1;
} else if (score < 0.75 * scoreOn) {
monitor.exceedCount =
Math.max(0, monitor.exceedCount - 2);
}
boolean persistentAlarm =
monitor.exceedCount * dt >= persistenceTime;
boolean dwellSatisfied =
time - monitor.lastResetTime >= dwellTime;
if (!monitor.alarm
&& persistentAlarm
&& dwellSatisfied) {
monitor.alarm = true;
for (int i = 0; i < theta.length; ++i) {
theta[i] =
(1.0 - resetFraction) * theta[i]
+ resetFraction * thetaSafe[i];
}
monitor.freezeCount = (int) (freezeTime / dt);
monitor.lastResetTime = time;
monitor.exceedCount = 0;
resetCount += 1;
System.out.printf(
"Reset at t=%6.2f s, score=%.3f, "
+ "theta=[%.4f, %.4f]%n",
time, score, theta[0], theta[1]
);
}
if (monitor.alarm && score < scoreOff) {
monitor.alarm = false;
}
writer.write(String.format(
Locale.US,
"%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,"
+ "%.12g,%.12g,%.12g,%.12g,%.12g,%.12g,"
+ "%.12g,%d,%d",
time,
x,
xModel,
reference,
trackingError,
control,
theta[0],
theta[1],
score,
Math.sqrt(monitor.errorSq),
parameterRatio,
lambdaMin,
monitor.saturation,
monitor.alarm ? 1 : 0,
monitor.freezeCount > 0 ? 1 : 0
));
writer.newLine();
}
System.out.println(
"Wrote " + steps
+ " samples to " + output.toAbsolutePath()
);
System.out.println(
"Number of reset events: " + resetCount
);
System.out.printf(
"Final estimate: kx=%.6f, kr=%.6f%n",
theta[0],
theta[1]
);
} catch (IOException exception) {
System.err.println(
"I/O error: " + exception.getMessage()
);
System.exit(1);
}
}
}
12. MATLAB and Simulink Implementation
The script uses core MATLAB operations and creates diagnostic plots. In Simulink, the continuous plant and reference model can be built with Integrator, Sum, Gain, and Saturation blocks. A MATLAB Function block can contain the projection, monitor, hysteresis, reset, and freeze logic. Control System Toolbox supports linear analysis; Simulink and Stateflow are useful for production-grade supervisor logic.
Chapter24_Lesson4.m
% Chapter24_Lesson4.m
% Monitored first-order MRAC with projection, drift monitoring, soft reset,
% hysteresis, dwell time, and temporary adaptation freeze.
%
% The script uses core MATLAB syntax only. It also maps directly to Simulink:
% 1. Integrator blocks for x and x_model.
% 2. Gain/Sum blocks for plant, model, and u = kx_hat*x + kr_hat*r.
% 3. Saturation block for the actuator.
% 4. MATLAB Function block for adaptation, monitoring, projection, and reset.
% 5. To Workspace blocks for the signals stored below.
%
% Run:
% Chapter24_Lesson4
clear;
clc;
rng(4, 'twister');
dt = 0.01;
finalTime = 120.0;
steps = floor(finalTime / dt);
aPlantNominal = -0.8;
aPlantChanged = -1.05;
bPlant = 1.1;
aModel = -1.5;
bModel = 1.5;
thetaSafe = [
(aModel - aPlantNominal) / bPlant;
bModel / bPlant
];
theta = thetaSafe;
lower = [-3.0; 0.2];
upper = [1.0; 3.0];
gamma = [10.0; 8.0];
leakage = 0.002;
deadZone = 0.005;
controlLimit = 4.0;
scoreOn = 0.65;
scoreOff = 0.36;
persistenceTime = 0.70;
dwellTime = 8.0;
freezeTime = 2.0;
resetFraction = 0.55;
betaFast = exp(-dt / 1.0);
betaGramian = exp(-dt / 3.0);
x = 0.0;
xModel = 0.0;
errorSq = 0.0;
rateSq = 0.0;
saturationEwma = 0.0;
g00 = 0.0;
g01 = 0.0;
g11 = 0.0;
exceedCount = 0;
alarm = false;
freezeCount = 0;
lastResetTime = -1.0e9;
resetCount = 0;
% Columns:
% time, x, x_model, reference, tracking_error, control, kx_hat, kr_hat,
% monitor_score, error_rms, parameter_ratio, lambda_min_gramian,
% saturation_ewma, alarm, frozen
logData = zeros(steps, 15);
for k = 1:steps
time = (k - 1) * dt;
if time < 75.0
aPlant = aPlantNominal;
else
aPlant = aPlantChanged;
end
if time < 85.0
reference = sin(0.45 * time) + 0.4 * sin(1.2 * time);
else
reference = 0.15 * sin(0.2 * time);
end
if time >= 35.0 && time < 65.0
sensorBias = 0.50;
else
sensorBias = 0.0;
end
measurementNoise = 0.025 * randn();
xModel = xModel + dt * (aModel * xModel + bModel * reference);
measuredX = x + sensorBias + measurementNoise;
unsaturatedControl = theta(1) * measuredX + theta(2) * reference;
control = min(controlLimit, max(-controlLimit, unsaturatedControl));
plantDisturbance = 0.12 * sin(2.3 * time);
x = x + dt * (aPlant * x + bPlant * control + plantDisturbance);
trackingError = measuredX - xModel;
phi = [measuredX; reference];
normalization = 1.0 + phi' * phi;
if freezeCount > 0
thetaDot = [0.0; 0.0];
freezeCount = freezeCount - 1;
else
thetaDot = -gamma .* phi * trackingError / normalization ...
- leakage * (theta - thetaSafe);
if abs(trackingError) <= deadZone
thetaDot = -leakage * (theta - thetaSafe);
end
end
theta = theta + dt * thetaDot;
theta = min(upper, max(lower, theta));
errorSq = betaFast * errorSq ...
+ (1.0 - betaFast) * trackingError^2;
rateSq = betaFast * rateSq ...
+ (1.0 - betaFast) * (thetaDot' * thetaDot);
saturated = double(abs(unsaturatedControl - control) > 1.0e-12);
saturationEwma = betaFast * saturationEwma ...
+ (1.0 - betaFast) * saturated;
g00 = betaGramian * g00 + (1.0 - betaGramian) * phi(1)^2;
g01 = betaGramian * g01 + (1.0 - betaGramian) * phi(1) * phi(2);
g11 = betaGramian * g11 + (1.0 - betaGramian) * phi(2)^2;
discriminant = sqrt(max(0.0, (g00 - g11)^2 + 4.0 * g01^2));
lambdaMin = max(0.0, 0.5 * (g00 + g11 - discriminant));
center = 0.5 * (lower + upper);
halfSpan = 0.5 * (upper - lower);
parameterRatio = sqrt(mean(((theta - center) ./ halfSpan).^2));
errorMetric = sqrt(errorSq) / 0.20;
rateMetric = sqrt(rateSq) / 1.50;
excitationDeficit = max(0.0, (0.04 - lambdaMin) / 0.04);
score = 0.40 * errorMetric ...
+ 0.25 * parameterRatio ...
+ 0.20 * rateMetric ...
+ 0.05 * saturationEwma ...
+ 0.10 * excitationDeficit;
if parameterRatio > 0.90
score = max(score, 1.10);
end
if score > scoreOn
exceedCount = exceedCount + 1;
elseif score < 0.75 * scoreOn
exceedCount = max(0, exceedCount - 2);
end
persistentAlarm = exceedCount * dt >= persistenceTime;
dwellSatisfied = time - lastResetTime >= dwellTime;
if ~alarm && persistentAlarm && dwellSatisfied
alarm = true;
theta = (1.0 - resetFraction) * theta ...
+ resetFraction * thetaSafe;
freezeCount = floor(freezeTime / dt);
lastResetTime = time;
exceedCount = 0;
resetCount = resetCount + 1;
fprintf( ...
'Reset at t=%6.2f s, score=%.3f, theta=[%.4f, %.4f]\n', ...
time, score, theta(1), theta(2) ...
);
end
if alarm && score < scoreOff
alarm = false;
end
logData(k, :) = [
time, ...
x, ...
xModel, ...
reference, ...
trackingError, ...
control, ...
theta(1), ...
theta(2), ...
score, ...
sqrt(errorSq), ...
parameterRatio, ...
lambdaMin, ...
saturationEwma, ...
double(alarm), ...
double(freezeCount > 0)
];
end
outputFile = 'Chapter24_Lesson4_results.csv';
fileId = fopen(outputFile, 'w');
if fileId < 0
error('Cannot open %s for writing.', outputFile);
end
fprintf(fileId, [ ...
'time,x,x_model,reference,tracking_error,control,' ...
'kx_hat,kr_hat,monitor_score,error_rms,parameter_ratio,' ...
'lambda_min_gramian,saturation_ewma,alarm,frozen\n' ...
]);
fclose(fileId);
dlmwrite( ...
outputFile, ...
logData, ...
'-append', ...
'delimiter', ',', ...
'precision', '%.12g' ...
);
fprintf('Wrote %d samples to %s\n', steps, outputFile);
fprintf('Number of reset events: %d\n', resetCount);
fprintf('Final estimate: kx=%.6f, kr=%.6f\n', theta(1), theta(2));
% Optional diagnostic plots.
figure('Name', 'Chapter 24 Lesson 4: Monitored MRAC');
subplot(3, 1, 1);
plot(logData(:, 1), logData(:, 2), ...
logData(:, 1), logData(:, 3), '--');
grid on;
ylabel('State');
legend('x', 'x_m');
subplot(3, 1, 2);
plot(logData(:, 1), logData(:, 7), ...
logData(:, 1), logData(:, 8));
grid on;
ylabel('Adaptive gains');
legend('k_x hat', 'k_r hat');
subplot(3, 1, 3);
plot(logData(:, 1), logData(:, 9));
hold on;
yline(scoreOn, '--');
grid on;
xlabel('Time (s)');
ylabel('Monitor score');
legend('score', 'reset threshold');
13. Wolfram Mathematica Implementation
Wolfram Language provides built-in matrix, symbolic, statistical, and
plotting operations. The downloadable
Chapter24_Lesson4.nb places the following program in an
executable notebook input cell.
Chapter24_Lesson4.nb
(* Chapter24_Lesson4.nb
Monitored first-order MRAC with projection, drift monitoring, soft reset,
hysteresis, dwell time, and temporary adaptation freeze. *)
ClearAll["Global`*"];
SeedRandom[4];
dt = 0.01;
finalTime = 120.0;
steps = Floor[finalTime/dt];
aPlantNominal = -0.8;
aPlantChanged = -1.05;
bPlant = 1.1;
aModel = -1.5;
bModel = 1.5;
thetaSafe = {
(aModel - aPlantNominal)/bPlant,
bModel/bPlant
};
theta = thetaSafe;
lower = {-3.0, 0.2};
upper = {1.0, 3.0};
gamma = {10.0, 8.0};
leakage = 0.002;
deadZone = 0.005;
controlLimit = 4.0;
scoreOn = 0.65;
scoreOff = 0.36;
persistenceTime = 0.70;
dwellTime = 8.0;
freezeTime = 2.0;
resetFraction = 0.55;
betaFast = Exp[-dt/1.0];
betaGramian = Exp[-dt/3.0];
x = 0.0;
xModel = 0.0;
errorSq = 0.0;
rateSq = 0.0;
saturationEwma = 0.0;
g00 = 0.0;
g01 = 0.0;
g11 = 0.0;
exceedCount = 0;
alarm = False;
freezeCount = 0;
lastResetTime = -10.^9;
resetCount = 0;
ClampValue[value_, low_, high_] := Max[low, Min[high, value]];
logData = ConstantArray[0.0, {steps, 15}];
For[k = 1, k <= steps, k++,
time = (k - 1) dt;
aPlant = If[time < 75.0, aPlantNominal, aPlantChanged];
reference = If[
time < 85.0,
Sin[0.45 time] + 0.4 Sin[1.2 time],
0.15 Sin[0.2 time]
];
sensorBias = If[35.0 <= time < 65.0, 0.50, 0.0];
measurementNoise = 0.025 RandomVariate[NormalDistribution[0.0, 1.0]];
xModel = xModel + dt (aModel xModel + bModel reference);
measuredX = x + sensorBias + measurementNoise;
unsaturatedControl = theta[[1]] measuredX + theta[[2]] reference;
control = ClampValue[unsaturatedControl, -controlLimit, controlLimit];
plantDisturbance = 0.12 Sin[2.3 time];
x = x + dt (aPlant x + bPlant control + plantDisturbance);
trackingError = measuredX - xModel;
phi = {measuredX, reference};
normalization = 1.0 + phi.phi;
If[freezeCount > 0,
thetaDot = {0.0, 0.0};
freezeCount--,
thetaDot =
-gamma phi trackingError/normalization
- leakage (theta - thetaSafe);
If[Abs[trackingError] <= deadZone,
thetaDot = -leakage (theta - thetaSafe)
];
];
theta = theta + dt thetaDot;
theta = MapThread[ClampValue, {theta, lower, upper}];
errorSq =
betaFast errorSq + (1.0 - betaFast) trackingError^2;
rateSq =
betaFast rateSq + (1.0 - betaFast) (thetaDot.thetaDot);
saturated = Boole[Abs[unsaturatedControl - control] > 10.^-12];
saturationEwma =
betaFast saturationEwma + (1.0 - betaFast) saturated;
g00 = betaGramian g00 + (1.0 - betaGramian) phi[[1]]^2;
g01 = betaGramian g01 + (1.0 - betaGramian) phi[[1]] phi[[2]];
g11 = betaGramian g11 + (1.0 - betaGramian) phi[[2]]^2;
discriminant = Sqrt[Max[0.0, (g00 - g11)^2 + 4.0 g01^2]];
lambdaMin = Max[0.0, 0.5 (g00 + g11 - discriminant)];
center = 0.5 (lower + upper);
halfSpan = 0.5 (upper - lower);
parameterRatio = Sqrt[Mean[((theta - center)/halfSpan)^2]];
errorMetric = Sqrt[errorSq]/0.20;
rateMetric = Sqrt[rateSq]/1.50;
excitationDeficit = Max[0.0, (0.04 - lambdaMin)/0.04];
score =
0.40 errorMetric
+ 0.25 parameterRatio
+ 0.20 rateMetric
+ 0.05 saturationEwma
+ 0.10 excitationDeficit;
If[parameterRatio > 0.90, score = Max[score, 1.10]];
If[score > scoreOn,
exceedCount++,
If[score < 0.75 scoreOn,
exceedCount = Max[0, exceedCount - 2]
]
];
persistentAlarm = exceedCount dt >= persistenceTime;
dwellSatisfied = time - lastResetTime >= dwellTime;
If[! alarm && persistentAlarm && dwellSatisfied,
alarm = True;
theta =
(1.0 - resetFraction) theta
+ resetFraction thetaSafe;
freezeCount = Floor[freezeTime/dt];
lastResetTime = time;
exceedCount = 0;
resetCount++;
Print[
Row[{
"Reset at t=", NumberForm[time, {6, 2}],
" s, score=", NumberForm[score, {4, 3}],
", theta=", NumberForm[theta, {6, 4}]
}]
];
];
If[alarm && score < scoreOff, alarm = False];
logData[[k]] = {
time,
x,
xModel,
reference,
trackingError,
control,
theta[[1]],
theta[[2]],
score,
Sqrt[errorSq],
parameterRatio,
lambdaMin,
saturationEwma,
Boole[alarm],
Boole[freezeCount > 0]
};
];
header = {
"time", "x", "x_model", "reference", "tracking_error", "control",
"kx_hat", "kr_hat", "monitor_score", "error_rms", "parameter_ratio",
"lambda_min_gramian", "saturation_ewma", "alarm", "frozen"
};
Export[
"Chapter24_Lesson4_results.csv",
Prepend[logData, header],
"CSV"
];
Print["Wrote ", steps, " samples to Chapter24_Lesson4_results.csv"];
Print["Number of reset events: ", resetCount];
Print[
"Final estimate: kx=", NumberForm[theta[[1]], {8, 6}],
", kr=", NumberForm[theta[[2]], {8, 6}]
];
statePlot = ListLinePlot[
{
logData[[All, {1, 2}]],
logData[[All, {1, 3}]]
},
PlotLegends -> {"x", "x_m"},
Frame -> True,
FrameLabel -> {"Time (s)", "State"},
PlotRange -> All,
ImageSize -> Large
];
gainPlot = ListLinePlot[
{
logData[[All, {1, 7}]],
logData[[All, {1, 8}]]
},
PlotLegends -> {"k_x hat", "k_r hat"},
Frame -> True,
FrameLabel -> {"Time (s)", "Adaptive gains"},
PlotRange -> All,
ImageSize -> Large
];
scorePlot = ListLinePlot[
{
logData[[All, {1, 9}]],
Table[{t, scoreOn}, {t, 0.0, finalTime - dt, dt}]
},
PlotLegends -> {"monitor score", "reset threshold"},
Frame -> True,
FrameLabel -> {"Time (s)", "Score"},
PlotRange -> All,
ImageSize -> Large
];
Column[{statePlot, gainPlot, scorePlot}]
14. Problems and Solutions
Problem 1: Linear Drift Under Constant Bias
A scalar adaptive law is \( \dot{\hat\theta}=-\gamma\phi e \). During a low-excitation operating interval, suppose \( \phi(t)=0.4 \), \( e(t)=0.01 \), \( \gamma=20 \), and \( \hat\theta(0)=1 \). Find \( \hat\theta(100) \).
Solution:
\[ \dot{\hat\theta} =-(20)(0.4)(0.01)=-0.08. \]
\[ \hat\theta(t)=1-0.08t, \qquad \hat\theta(100)=1-8=-7. \]
The tracking error is only 0.01, but its nonzero mean produces a large parameter change over a long interval.
Problem 2: Ultimate Bound with Leakage
Suppose a Lyapunov analysis gives
\[ \dot V\leq-0.4V+0.02, \qquad V(0)=1. \]
Derive an upper bound on \( V(t) \) and its ultimate value.
Solution:
By the comparison lemma,
\[ V(t) \leq e^{-0.4t}V(0) +\frac{0.02}{0.4} \left(1-e^{-0.4t}\right). \]
\[ V(t) \leq 0.05+0.95e^{-0.4t}, \qquad \limsup_{t\to\infty}V(t)\leq0.05. \]
Problem 3: Soft-Reset Identity
Prove that for \( \mathbf{z}^{+}=(1-\rho)\mathbf{a}+\rho\mathbf{b} \) and \( \mathbf{W}>0 \),
\[ \|\mathbf{z}^{+}\|_{\mathbf{W}}^2 = (1-\rho)\|\mathbf{a}\|_{\mathbf{W}}^2 +\rho\|\mathbf{b}\|_{\mathbf{W}}^2 -\rho(1-\rho) \|\mathbf{a}-\mathbf{b}\|_{\mathbf{W}}^2. \]
Solution:
Expand the left side:
\[ \begin{aligned} \|\mathbf{z}^{+}\|_{\mathbf{W}}^2 &= (1-\rho)^2\mathbf{a}^{T}\mathbf{W}\mathbf{a} +2\rho(1-\rho)\mathbf{a}^{T}\mathbf{W}\mathbf{b} +\rho^2\mathbf{b}^{T}\mathbf{W}\mathbf{b}. \end{aligned} \]
Also,
\[ \|\mathbf{a}-\mathbf{b}\|_{\mathbf{W}}^2 = \mathbf{a}^{T}\mathbf{W}\mathbf{a} -2\mathbf{a}^{T}\mathbf{W}\mathbf{b} +\mathbf{b}^{T}\mathbf{W}\mathbf{b}. \]
Substitution and collection of coefficients gives the required identity.
Problem 4: Covariance Reset Preserves Positive Definiteness
Let \( \mathbf{P}^{-}>0 \) and \( \mathbf{P}_0>0 \). Show that
\[ \mathbf{P}^{+} =(1-\alpha)\mathbf{P}^{-} +\alpha\mathbf{P}_0, \qquad 0\leq\alpha\leq1, \]
is positive definite.
Solution:
For every nonzero vector \( \mathbf{v} \),
\[ \mathbf{v}^{T}\mathbf{P}^{+}\mathbf{v} = (1-\alpha) \mathbf{v}^{T}\mathbf{P}^{-}\mathbf{v} +\alpha \mathbf{v}^{T}\mathbf{P}_0\mathbf{v}. \]
Both quadratic forms are positive. Their convex combination is therefore positive, so \( \mathbf{P}^{+}>0 \).
Problem 5: Persistence and Dwell-Time Samples
A controller samples at \( T_s=5\text{ ms} \). The alarm score must remain above threshold for \( T_p=0.75\text{ s} \), and resets must be separated by at least \( T_d=10\text{ s} \). Find the required sample counts.
Solution:
\[ N_p= \left\lceil\frac{0.75}{0.005}\right\rceil =150, \qquad N_d= \left\lceil\frac{10}{0.005}\right\rceil =2000. \]
The implementation should require 150 persistent alarm samples and reject a new reset until at least 2000 samples have elapsed after the preceding reset.
Problem 6: Distinguishing Adaptation from Drift
A parameter estimate changes rapidly after a load change. The tracking residual initially increases, the information-matrix minimum eigenvalue is large, saturation remains zero, and the residual decreases as the estimate settles inside its admissible set. Should the supervisor reset the estimate?
Solution:
Not solely because the parameter moved rapidly. The regressor is informative, the actuator is not saturated, the estimate remains admissible, and the residual improves. These observations are consistent with useful adaptation to a real plant change. A reset would discard information and could worsen tracking. The monitor should distinguish a transient adaptation event from persistent, unexplained motion with poor residual improvement.
15. Summary
Parameter drift is a long-duration robustness problem, not merely a large instantaneous error. Ordinary gradient adaptation can integrate small biased signals because its Lyapunov derivative need not contain a strictly negative parameter-error term. Leakage creates an ultimate bound, dead zones suppress noise-driven updates, normalization limits regressor-induced update size, and projection enforces admissible parameter regions. A practical implementation should monitor tracking residuals, update rate, parameter-bound proximity, finite-window excitation, saturation, timing, and sensor validity. Reset logic must use persistence, hysteresis, dwell time, adaptation freeze, and bumpless transfer. A reset is not automatically stabilizing; its jump map and fallback logic require explicit analysis and validation.
16. References
- Peterson, B.B., & Narendra, K.S. (1982). Bounded error adaptive control. IEEE Transactions on Automatic Control, 27(6), 1161–1168.
- Kreisselmeier, G., & Narendra, K.S. (1982). Stable model reference adaptive control in the presence of bounded disturbances. IEEE Transactions on Automatic Control, 27(6), 1169–1175.
- Åström, K.J. (1983). Theory and applications of adaptive control—A survey. Automatica, 19(5), 471–486.
- Ioannou, P.A., & Kokotović, P.V. (1984). Instability analysis and improvement of robustness of adaptive control. Automatica, 20(5), 583–594.
- 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.
- Ioannou, P.A., & Tsakalis, K.S. (1986). A robust direct adaptive controller. IEEE Transactions on Automatic Control, 31(11), 1033–1043.
- Narendra, K.S., & Annaswamy, A.M. (1987). A new adaptive law for robust adaptation without persistent excitation. IEEE Transactions on Automatic Control, 32(2), 134–145.
- Morse, A.S., Mayne, D.Q., & Goodwin, G.C. (1992). Applications of hysteresis switching in parameter adaptive control. IEEE Transactions on Automatic Control, 37(9), 1343–1354.
- Morse, A.S. (1996). Supervisory control of families of linear set-point controllers—Part I: Exact matching. IEEE Transactions on Automatic Control, 41(10), 1413–1431.
- Narendra, K.S., & Balakrishnan, J. (1997). Adaptive control using multiple models. IEEE Transactions on Automatic Control, 42(2), 171–187.
- Hespanha, J.P., Liberzon, D., & Morse, A.S. (2003). Overcoming the limitations of adaptive control by means of logic-based switching. Systems & Control Letters, 49(1), 49–65.
Help keep these engineering tutorials free and growing
If these lessons, examples, and project pages help you, a small donation supports the continued creation and improvement of free control, robotics, software, and engineering education resources.
Created and maintained by Abolfazl Mohammadijoo.
