In Lessons 1–3 of this chapter, composite adaptation combined two kinds of information: a tracking-error channel, which corrects parameters because the controlled state is not following its reference, and a prediction-error channel, which corrects parameters because the uncertain model does not reproduce measured data. DREM is best viewed as a signal-processing and regression-restructuring tool for the second channel. It is not, by itself, a complete feedback controller.

The starting point is a measurable regression that is linear in an unknown constant parameter vector:

\[ y(t)=\boldsymbol{\phi}^{\top}(t)\boldsymbol{\theta}, \qquad y(t)\in\mathbb{R},\quad \boldsymbol{\phi}(t)\in\mathbb{R}^{q},\quad \boldsymbol{\theta}\in\mathbb{R}^{q}. \]

DREM transforms this one scalar equation with a vector regressor into q scalar equations, each containing only one unknown parameter. This decoupling is the central structural result of the method.

Learning objectives

After this lesson, a student should be able to:

• construct an extended regression with stable dynamic operators;
• derive determinant-adjugate mixing for a general parameter dimension;
• prove elementwise monotonicity of the ideal DREM estimation errors;
• distinguish vector persistent excitation from the DREM scalar condition;
• embed a DREM prediction term into a composite adaptive law;
• recognize filter-transient, noise, conditioning, and discretization issues.

flowchart TD
  A["Measured plant signals"] --> B["Linear regression: y = phi^T theta"]
  B --> C["Dynamic extension: original plus filtered regressions"]
  C --> D["Square extended regressor matrix"]
  D --> E["Mixing: determinant and adjugate"]
  E --> F["q scalar regressions"]
  F --> G["Scalar parameter update laws"]
  G --> H["Prediction term for composite adaptation"]
        

For the regression above, a classical gradient estimator is

\[ \dot{\hat{\boldsymbol{\theta}}} =\boldsymbol{\Gamma}\boldsymbol{\phi} \left(y-\boldsymbol{\phi}^{\top}\hat{\boldsymbol{\theta}}\right), \qquad \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{\top}>0. \]

Define the parameter error using the convention \( \tilde{\boldsymbol{\theta}} =\hat{\boldsymbol{\theta}}-\boldsymbol{\theta} \). Since \( y-\boldsymbol{\phi}^{\top}\hat{\boldsymbol{\theta}} =-\boldsymbol{\phi}^{\top}\tilde{\boldsymbol{\theta}} \), its dynamics are

\[ \dot{\tilde{\boldsymbol{\theta}}} =-\boldsymbol{\Gamma}\boldsymbol{\phi}\boldsymbol{\phi}^{\top} \tilde{\boldsymbol{\theta}}. \]

With \( V_{\theta}=\tfrac{1}{2} \tilde{\boldsymbol{\theta}}^{\top}\boldsymbol{\Gamma}^{-1} \tilde{\boldsymbol{\theta}} \),

\[ \dot V_{\theta} =-\left(\boldsymbol{\phi}^{\top} \tilde{\boldsymbol{\theta}}\right)^2\leq 0. \]

Thus a weighted norm of the vector error is nonincreasing. However, the individual coordinates are coupled by \( \boldsymbol{\phi}\boldsymbol{\phi}^{\top} \). A component such as \( |\tilde\theta_1| \) can increase temporarily even while the total Lyapunov function decreases. Moreover, standard exponential convergence is normally associated with a vector persistent-excitation inequality:

\[ \exists T>0,\ \alpha>0:\qquad \int_{t}^{t+T} \boldsymbol{\phi}(\tau)\boldsymbol{\phi}^{\top}(\tau)\,d\tau \geq \alpha\mathbf{I}_{q} \quad \forall t\geq 0. \]

DREM changes the geometry of the estimation problem. It does not merely retune \( \boldsymbol{\Gamma} \); it constructs a new regression whose parameter channels are scalar and decoupled.

Select \( q-1 \) linear, bounded-input bounded-output stable operators \( \mathcal{H}_1,\ldots,\mathcal{H}_{q-1} \). Typical choices are stable first-order filters or time-delay operators. Apply every operator to both the scalar output and every component of the regressor:

\[ y_i(t)=\mathcal{H}_i[y](t), \qquad \boldsymbol{\phi}_i(t)=\mathcal{H}_i[\boldsymbol{\phi}](t), \qquad i=1,\ldots,q-1. \]

Because \( \boldsymbol{\theta} \) is constant, an ideal operator with compatible initial conditions preserves the regression:

\[ y_i(t)=\boldsymbol{\phi}_i^{\top}(t)\boldsymbol{\theta}. \]

With ordinary filter initialization, an exponentially decaying transient is generally present:

\[ y_i(t)=\boldsymbol{\phi}_i^{\top}(t)\boldsymbol{\theta} +\varepsilon_i(t), \qquad |\varepsilon_i(t)|\leq c_i e^{-a_i t}. \]

Stack the original and generated equations:

\[ \mathbf{Y}_{e}(t)= \begin{bmatrix} y(t)\\ y_1(t)\\ \vdots\\ y_{q-1}(t) \end{bmatrix}, \qquad \mathbf{M}_{e}(t)= \begin{bmatrix} \boldsymbol{\phi}^{\top}(t)\\ \boldsymbol{\phi}_1^{\top}(t)\\ \vdots\\ \boldsymbol{\phi}_{q-1}^{\top}(t) \end{bmatrix}. \]

The extended vector regression becomes

\[ \mathbf{Y}_{e}(t)=\mathbf{M}_{e}(t)\boldsymbol{\theta} +\boldsymbol{\varepsilon}_{e}(t). \]

First-order filter realization

A common operator is \( \mathcal{H}(s)=\lambda/(s+\lambda) \), with \( \lambda>0 \). It is implemented without signal differentiation:

\[ \dot z_y=-\lambda z_y+\lambda y, \qquad \dot{\mathbf{z}}_{\phi}=-\lambda\mathbf{z}_{\phi} +\lambda\boldsymbol{\phi}, \]

where \( y_f=z_y \) and \( \boldsymbol{\phi}_f=\mathbf{z}_{\phi} \). For a q-parameter regression, enough distinct operators must be used to form a square \( q\times q \) matrix.

Define the scalar extended regressor

\[ \Delta(t)=\det\!\left(\mathbf{M}_{e}(t)\right). \]

Premultiply the extended regression by the adjugate matrix:

\[ \boldsymbol{\mathcal{Y}}(t) =\operatorname{adj}\!\left(\mathbf{M}_{e}(t)\right) \mathbf{Y}_{e}(t). \]

The algebraic identity \( \operatorname{adj}(\mathbf{M})\mathbf{M} =\det(\mathbf{M})\mathbf{I} \) gives

\[ \boldsymbol{\mathcal{Y}}(t) =\Delta(t)\boldsymbol{\theta} +\operatorname{adj}\!\left(\mathbf{M}_{e}(t)\right) \boldsymbol{\varepsilon}_{e}(t). \]

Ignoring the decaying filter transient for the moment, the vector equation is equivalent to q scalar regressions:

\[ \mathcal{Y}_i(t)=\Delta(t)\theta_i, \qquad i=1,\ldots,q. \]

This is the mixing step. Every parameter now appears in its own scalar equation, while all excitation information is concentrated into the one scalar signal \( \Delta(t) \).

No matrix inverse is required

Although the determinant and adjugate resemble the inverse formula \( \mathbf{M}^{-1}=\operatorname{adj}(\mathbf{M})/\det(\mathbf{M}) \), the standard DREM gradient law does not divide by \( \Delta \). Therefore isolated determinant zeros do not create algebraic singularities. They merely suspend or slow parameter adaptation at those instants.

For each scalar regression, use

\[ \dot{\hat\theta}_i =\gamma_i\Delta \left(\mathcal{Y}_i-\Delta\hat\theta_i\right), \qquad \gamma_i>0. \]

With \( \tilde\theta_i=\hat\theta_i-\theta_i \) and the ideal relation \( \mathcal{Y}_i=\Delta\theta_i \),

\[ \dot{\tilde\theta}_i =-\gamma_i\Delta^2(t)\tilde\theta_i. \]

Closed-form solution

Separating variables gives

\[ \frac{d\tilde\theta_i}{\tilde\theta_i} =-\gamma_i\Delta^2(t)dt. \]

Hence, for any initial time \( t_0 \),

\[ \tilde\theta_i(t) =\exp\!\left[-\gamma_i \int_{t_0}^{t}\Delta^2(\tau)d\tau\right] \tilde\theta_i(t_0). \]

This identity immediately establishes three properties.

Property 1 — sign preservation

\[ \operatorname{sign}\!\left(\tilde\theta_i(t)\right) =\operatorname{sign}\!\left(\tilde\theta_i(t_0)\right) \]

unless the initial error is zero. The ideal scalar error approaches zero without crossing it.

Property 2 — elementwise monotonicity

\[ |\tilde\theta_i(t)| \leq |\tilde\theta_i(t_0)|, \qquad t\geq t_0. \]

This is stronger than merely proving that a norm of the full parameter error vector is nonincreasing.

Property 3 — asymptotic convergence under a scalar condition

\[ \int_{t_0}^{\infty}\Delta^2(\tau)d\tau=\infty \quad \Longrightarrow \quad \lim_{t\rightarrow\infty}\tilde\theta_i(t)=0. \]

Equivalently, the basic ideal DREM estimator converges when \( \Delta\notin\mathcal{L}_2 \). This scalar non-square-integrability condition is the signature convergence test of the basic construction.

Excitation statements must be made carefully. DREM does not imply that every weakly varying input identifies every parameter. The chosen operators and the original regressor must produce a determinant carrying sufficient energy.

6.1 Asymptotic convergence of basic DREM

\[ \Delta\notin\mathcal{L}_2 \quad \Longleftrightarrow \quad \int_{0}^{\infty}\Delta^2(\tau)d\tau=\infty. \]

This condition is scalar and may hold even when the original vector regressor is not persistently exciting in the classical uniform sense. However, it is a condition on the constructed determinant, not an automatic consequence of filtering.

6.2 Exponential convergence

A sufficient scalar persistent-excitation condition is

\[ \exists T>0,\ \alpha_{\Delta}>0:\qquad \int_{t}^{t+T}\Delta^2(\tau)d\tau \geq \alpha_{\Delta} \quad \forall t\geq 0. \]

Under this condition, the integral in the error solution grows at least linearly with time, and each scalar error converges exponentially.

6.3 Interval excitation and finite-time variants

Some later DREM variants add memory, resetting, or modified mixing and can obtain finite-time recovery after excitation over a finite interval. That result should not be attributed automatically to the elementary gradient estimator derived in Section 5. For the basic estimator, a finite burst with finite \( \int_0^{\infty}\Delta^2 d\tau \) reduces the error by a finite factor but does not force exact asymptotic convergence.

6.4 A simple failure case

If every row of \( \mathbf{M}_e(t) \) remains linearly dependent, then \( \Delta(t)\equiv 0 \). The mixed regressions contain no adaptation signal and \( \dot{\hat\theta}_i=0 \). Filter diversity is therefore essential.

Consider a tracking-error model of the familiar linearly parameterized form

\[ \dot{\mathbf{e}} =\mathbf{A}_{m}\mathbf{e} +\mathbf{B}\boldsymbol{\psi}^{\top}(t) \tilde{\boldsymbol{\theta}}, \qquad \mathbf{A}_{m}\ \text{Hurwitz}. \]

Let \( \mathbf{P}=\mathbf{P}^{\top}>0 \) solve

\[ \mathbf{A}_{m}^{\top}\mathbf{P} +\mathbf{P}\mathbf{A}_{m}=-\mathbf{Q}, \qquad \mathbf{Q}=\mathbf{Q}^{\top}>0. \]

A composite parameter update can combine the tracking term and a DREM prediction term:

\[ \dot{\hat{\boldsymbol{\theta}}} =-\boldsymbol{\Gamma}\boldsymbol{\psi} \mathbf{e}^{\top}\mathbf{P}\mathbf{B} +\boldsymbol{\Gamma}\mathbf{K}_{D}\Delta \left(\boldsymbol{\mathcal{Y}} -\Delta\hat{\boldsymbol{\theta}}\right), \]

where \( \mathbf{K}_{D}=\operatorname{diag}(k_{D1},\ldots,k_{Dq}) \) and \( k_{Di}\geq 0 \). Use the composite Lyapunov function

\[ V=\mathbf{e}^{\top}\mathbf{P}\mathbf{e} +\tilde{\boldsymbol{\theta}}^{\top} \boldsymbol{\Gamma}^{-1} \tilde{\boldsymbol{\theta}}. \]

In the ideal regression case, the tracking update cancels the state-parameter cross term and the DREM term contributes a negative semidefinite parameter-error term:

\[ \dot V =-\mathbf{e}^{\top}\mathbf{Q}\mathbf{e} -2\Delta^2 \tilde{\boldsymbol{\theta}}^{\top} \mathbf{K}_{D}\tilde{\boldsymbol{\theta}} \leq 0. \]

Conceptually, tracking error protects closed-loop behavior, while the mixed prediction regressions inject direct parameter information. The channels reinforce each other but play different roles.

flowchart TD
  R["Reference model"] --> E["Tracking error"]
  P["Uncertain plant"] --> E
  E --> T["Tracking-error adaptation term"]
  P --> L["Measured-data regression"]
  L --> X["Extension filters"]
  X --> M["Determinant-adjugate mixing"]
  M --> D["DREM prediction term"]
  T --> U["Composite parameter update"]
  D --> U
  U --> C["Adaptive controller parameters"]
  C --> P
        

Let

\[ y=\phi_1\theta_1+\phi_2\theta_2. \]

Apply one stable filter to obtain \( y_f=\phi_{1f}\theta_1+\phi_{2f}\theta_2 \). The extended regression is

\[ \underbrace{ \begin{bmatrix}y\\y_f\end{bmatrix}}_{\mathbf{Y}_e} = \underbrace{ \begin{bmatrix} \phi_1 & \phi_2\\ \phi_{1f} & \phi_{2f} \end{bmatrix}}_{\mathbf{M}_e} \begin{bmatrix}\theta_1\\\theta_2\end{bmatrix}. \]

The determinant is

\[ \Delta=\phi_1\phi_{2f}-\phi_2\phi_{1f}. \]

For a two-by-two matrix,

\[ \operatorname{adj}(\mathbf{M}_e) = \begin{bmatrix} \phi_{2f} & -\phi_2\\ -\phi_{1f} & \phi_1 \end{bmatrix}. \]

Therefore,

\[ \begin{aligned} \mathcal{Y}_1&=\phi_{2f}y-\phi_2y_f=\Delta\theta_1,\\ \mathcal{Y}_2&=-\phi_{1f}y+\phi_1y_f=\Delta\theta_2. \end{aligned} \]

The two update laws are

\[ \begin{aligned} \dot{\hat\theta}_1 &=\gamma_1\Delta \left(\mathcal{Y}_1-\Delta\hat\theta_1\right),\\ \dot{\hat\theta}_2 &=\gamma_2\Delta \left(\mathcal{Y}_2-\Delta\hat\theta_2\right). \end{aligned} \]

Observe that neither equation contains the other parameter estimate. Their only shared time-varying quantity is \( \Delta \).

Geometric interpretation

The two rows of \( \mathbf{M}_e \) are two regressor directions in the parameter plane. The magnitude \( |\Delta| \) equals the area of the parallelogram generated by these rows. If the rows are nearly parallel, the area is small and adaptation is slow. If they remain exactly parallel, the area is zero and the two parameters cannot be separated by this extension.

9.1 Distinct dynamics are required

Repeating identical filters creates duplicate rows and may force the determinant to zero. Common choices include different first-order poles, different delays, or combinations of stable filters and delays.

\[ \mathcal{H}_i(s)=\frac{\lambda_i}{s+\lambda_i}, \qquad \lambda_i>0, \qquad \lambda_i\neq\lambda_j\ \text{for selected channels}. \]

9.2 Fast versus slow filters

Very fast filters make filtered regressors closely resemble the original regressor, which can reduce row diversity. Very slow filters create stronger phase separation but increase transient duration and lag. Pole selection is therefore an excitation-conditioning trade-off rather than a purely bandwidth-based choice.

9.3 Determinant scaling

For large parameter dimension, determinants and cofactors may become very large or very small. A normalized scalar update can be used:

\[ \dot{\hat\theta}_i =\gamma_i\frac{\Delta}{1+\nu\Delta^2} \left(\mathcal{Y}_i-\Delta\hat\theta_i\right), \qquad \nu>0. \]

In the ideal case, this gives

\[ \dot{\tilde\theta}_i =-\gamma_i\frac{\Delta^2}{1+\nu\Delta^2} \tilde\theta_i, \]

preserving the scalar sign and monotonicity properties while limiting effective adaptation gain for very large \( |\Delta| \).

9.4 Monitoring an excitation certificate

A practical implementation can accumulate

\[ S_{\Delta}(t)=\int_{0}^{t}\Delta^2(\tau)d\tau, \qquad \dot S_{\Delta}=\Delta^2, \]

which measures the total determinant energy observed so far. Growth of \( S_{\Delta} \) is informative, but it is not by itself a proof of future excitation.

10.1 Filter initialization error

With nonzero transient terms, the mixed scalar regression is

\[ \mathcal{Y}_i =\Delta\theta_i+d_i(t), \qquad \mathbf{d}(t)=\operatorname{adj}(\mathbf{M}_e) \boldsymbol{\varepsilon}_e(t). \]

The parameter error then satisfies

\[ \dot{\tilde\theta}_i =-\gamma_i\Delta^2\tilde\theta_i +\gamma_i\Delta d_i(t). \]

The perfect sign-preservation statement belongs to the unperturbed model. A decaying or noisy \( d_i \) can cause small crossings or residual fluctuations.

10.2 Measurement noise

Filtering attenuates some high-frequency noise, but the adjugate contains cofactors formed from products of regressor entries. Consequently, mixing can amplify noise when the extended matrix is poorly scaled. Implementations should log \( \Delta \), cofactor magnitudes, residuals, and parameter-update rates.

10.3 Leakage, dead zone, and projection

The robust modifications introduced earlier in the course remain applicable. One possible scalar robust law is

\[ \dot{\hat\theta}_i =\gamma_i\Delta \operatorname{dz}_{\epsilon_i} \left(\mathcal{Y}_i-\Delta\hat\theta_i\right) -\sigma_i\hat\theta_i, \]

followed by projection when known parameter bounds exist. Leakage and a dead zone improve robustness but generally replace exact convergence with an ultimate bound.

10.4 Discrete-time implementation

With sampling period \( h \), a forward-Euler scalar update is

\[ \hat\theta_i[k+1] =\hat\theta_i[k] +h\gamma_i\Delta[k] \left(\mathcal{Y}_i[k]-\Delta[k]\hat\theta_i[k]\right). \]

For a frozen value of \( \Delta[k] \), the scalar error multiplier is \( 1-h\gamma_i\Delta^2[k] \). A conservative numerical requirement is

\[ 0\leq h\gamma_i\Delta^2[k]\leq 1 \]

when discrete elementwise non-overshoot is desired. The wider linear stability interval is \( 0<h\gamma_i\Delta^2[k]<2 \), but values above one permit sign alternation of the discrete error.

The downloadable programs use

\[ \boldsymbol{\theta} =\begin{bmatrix}2\\-1.5\end{bmatrix}, \qquad \phi_1(t)=1+0.5\sin(0.7t), \]

\[ \phi_2(t)=0.7\cos(1.3t)+0.2\sin(2.4t), \qquad y(t)=\boldsymbol{\phi}^{\top}(t)\boldsymbol{\theta}. \]

One first-order extension filter is used because the parameter dimension is two:

\[ \mathcal{H}(s)=\frac{1.5}{s+1.5}. \]

The adaptation gains are \( \gamma_1=\gamma_2=30 \), the simulation interval is \( 0\leq t\leq25 \) seconds, and the Euler step is \( h=0.002 \) seconds. Each program exports the original and filtered regressors, \( \Delta \), mixed outputs, parameter estimates, errors, and accumulated determinant energy to CSV.

Algorithm

1. sample \( \boldsymbol{\phi}(t) \) and compute \( y(t) \);
2. update filters for \( y_f \) and \( \boldsymbol{\phi}_f \);
3. compute \( \Delta=\phi_1\phi_{2f}-\phi_2\phi_{1f} \);
4. compute \( \mathcal{Y}_1,\mathcal{Y}_2 \) with the adjugate formulas;
5. integrate the two scalar update laws;
6. log errors and \( \int\Delta^2 \).

Chapter21_Lesson4.py

"""Chapter 21, Lesson 4: Dynamic Regressor Extension and Mixing.

A dependency-free NumPy demonstration of a two-parameter DREM estimator.
The original scalar regression is y(t) = phi(t)^T theta. A stable first-order
filter creates a second regression. The 2 x 2 extended regressor is then mixed
with its adjugate to produce two scalar regressions.
"""

from __future__ import annotations

import csv
import math
from dataclasses import dataclass
from pathlib import Path

import numpy as np


@dataclass(frozen=True)
class Configuration:
    dt: float = 0.002
    final_time: float = 25.0
    filter_rate: float = 1.5
    gamma_1: float = 30.0
    gamma_2: float = 30.0
    output_csv: str = "Chapter21_Lesson4_results.csv"


def regressor(t: float) -> np.ndarray:
    """Known bounded regressor used by the numerical experiment."""
    return np.array(
        [
            1.0 + 0.5 * math.sin(0.7 * t),
            0.7 * math.cos(1.3 * t) + 0.2 * math.sin(2.4 * t),
        ],
        dtype=float,
    )


def simulate(config: Configuration) -> tuple[np.ndarray, np.ndarray]:
    theta_true = np.array([2.0, -1.5], dtype=float)
    theta_hat = np.zeros(2, dtype=float)
    gamma = np.array([config.gamma_1, config.gamma_2], dtype=float)

    filtered_y = 0.0
    filtered_phi = np.zeros(2, dtype=float)
    integral_delta_squared = 0.0

    rows: list[list[float]] = []
    number_of_steps = int(round(config.final_time / config.dt))

    for step in range(number_of_steps + 1):
        time = step * config.dt
        phi = regressor(time)
        y = float(phi @ theta_true)

        # Stable strictly proper extension operator H(s) = lambda / (s + lambda).
        filtered_y += config.dt * (
            -config.filter_rate * filtered_y + config.filter_rate * y
        )
        filtered_phi += config.dt * (
            -config.filter_rate * filtered_phi + config.filter_rate * phi
        )

        # Phi_e = [[phi_1, phi_2], [phi_1f, phi_2f]].
        delta = phi[0] * filtered_phi[1] - phi[1] * filtered_phi[0]

        # mathcalY = adj(Phi_e) Y_e for a 2 x 2 matrix.
        mixed_1 = filtered_phi[1] * y - phi[1] * filtered_y
        mixed_2 = -filtered_phi[0] * y + phi[0] * filtered_y
        mixed = np.array([mixed_1, mixed_2], dtype=float)

        # Two decoupled scalar gradient estimators.
        theta_hat += (
            config.dt * gamma * delta * (mixed - delta * theta_hat)
        )
        integral_delta_squared += config.dt * delta * delta

        rows.append(
            [
                time,
                phi[0],
                phi[1],
                y,
                filtered_phi[0],
                filtered_phi[1],
                filtered_y,
                delta,
                mixed_1,
                mixed_2,
                theta_hat[0],
                theta_hat[1],
                theta_true[0] - theta_hat[0],
                theta_true[1] - theta_hat[1],
                integral_delta_squared,
            ]
        )

    output_path = Path(config.output_csv)
    with output_path.open("w", newline="", encoding="utf-8") as csv_file:
        writer = csv.writer(csv_file)
        writer.writerow(
            [
                "time",
                "phi1",
                "phi2",
                "y",
                "phi1_filtered",
                "phi2_filtered",
                "y_filtered",
                "delta",
                "mixed1",
                "mixed2",
                "theta1_hat",
                "theta2_hat",
                "theta1_error",
                "theta2_error",
                "integral_delta_squared",
            ]
        )
        writer.writerows(rows)

    return theta_true, theta_hat


def main() -> None:
    config = Configuration()
    theta_true, theta_hat = simulate(config)
    error = theta_true - theta_hat

    print("DREM two-parameter simulation")
    print(f"true theta      = {theta_true}")
    print(f"estimated theta = {theta_hat}")
    print(f"final error     = {error}")
    print(f"results written to {config.output_csv}")


if __name__ == "__main__":
    main()

Chapter21_Lesson4.cpp

// Chapter 21, Lesson 4: Dynamic Regressor Extension and Mixing.
// Build: g++ -std=c++17 -O2 Chapter21_Lesson4.cpp -o Chapter21_Lesson4

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

namespace {
constexpr double kDt = 0.002;
constexpr double kFinalTime = 25.0;
constexpr double kFilterRate = 1.5;
constexpr std::array<double, 2> kThetaTrue{2.0, -1.5};
constexpr std::array<double, 2> kGamma{30.0, 30.0};

std::array<double, 2> regressor(const double t) {
    return {
        1.0 + 0.5 * std::sin(0.7 * t),
        0.7 * std::cos(1.3 * t) + 0.2 * std::sin(2.4 * t)
    };
}

double dot(const std::array<double, 2>& a,
           const std::array<double, 2>& b) {
    return a[0] * b[0] + a[1] * b[1];
}
}  // namespace

int main() {
    std::array<double, 2> thetaHat{0.0, 0.0};
    std::array<double, 2> filteredPhi{0.0, 0.0};
    double filteredY = 0.0;
    double integralDeltaSquared = 0.0;

    std::ofstream csv("Chapter21_Lesson4_results.csv");
    if (!csv) {
        throw std::runtime_error("Cannot open Chapter21_Lesson4_results.csv");
    }
    csv << std::setprecision(15);
    csv << "time,phi1,phi2,y,phi1_filtered,phi2_filtered,y_filtered,"
           "delta,mixed1,mixed2,theta1_hat,theta2_hat,theta1_error,"
           "theta2_error,integral_delta_squared\n";

    const int steps = static_cast<int>(std::llround(kFinalTime / kDt));
    for (int step = 0; step <= steps; ++step) {
        const double t = step * kDt;
        const auto phi = regressor(t);
        const double y = dot(phi, kThetaTrue);

        filteredY += kDt * (-kFilterRate * filteredY + kFilterRate * y);
        for (std::size_t i = 0; i < filteredPhi.size(); ++i) {
            filteredPhi[i] +=
                kDt * (-kFilterRate * filteredPhi[i] + kFilterRate * phi[i]);
        }

        const double delta =
            phi[0] * filteredPhi[1] - phi[1] * filteredPhi[0];
        const std::array<double, 2> mixed{
            filteredPhi[1] * y - phi[1] * filteredY,
            -filteredPhi[0] * y + phi[0] * filteredY
        };

        for (std::size_t i = 0; i < thetaHat.size(); ++i) {
            thetaHat[i] +=
                kDt * kGamma[i] * delta * (mixed[i] - delta * thetaHat[i]);
        }
        integralDeltaSquared += kDt * delta * delta;

        csv << t << ',' << phi[0] << ',' << phi[1] << ',' << y << ','
            << filteredPhi[0] << ',' << filteredPhi[1] << ',' << filteredY
            << ',' << delta << ',' << mixed[0] << ',' << mixed[1] << ','
            << thetaHat[0] << ',' << thetaHat[1] << ','
            << kThetaTrue[0] - thetaHat[0] << ','
            << kThetaTrue[1] - thetaHat[1] << ','
            << integralDeltaSquared << '\n';
    }

    std::cout << std::setprecision(10)
              << "DREM two-parameter simulation\n"
              << "true theta      = [" << kThetaTrue[0] << ", "
              << kThetaTrue[1] << "]\n"
              << "estimated theta = [" << thetaHat[0] << ", "
              << thetaHat[1] << "]\n"
              << "final error     = [" << kThetaTrue[0] - thetaHat[0]
              << ", " << kThetaTrue[1] - thetaHat[1] << "]\n"
              << "results written to Chapter21_Lesson4_results.csv\n";

    return 0;
}

Chapter21_Lesson4.java

// Chapter 21, Lesson 4: Dynamic Regressor Extension and Mixing.
// Build: javac Chapter21_Lesson4.java
// Run:   java Chapter21_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;

public final class Chapter21_Lesson4 {
    private static final double DT = 0.002;
    private static final double FINAL_TIME = 25.0;
    private static final double FILTER_RATE = 1.5;
    private static final double[] THETA_TRUE = {2.0, -1.5};
    private static final double[] GAMMA = {30.0, 30.0};

    private Chapter21_Lesson4() {
    }

    private static double[] regressor(double time) {
        return new double[] {
            1.0 + 0.5 * Math.sin(0.7 * time),
            0.7 * Math.cos(1.3 * time) + 0.2 * Math.sin(2.4 * time)
        };
    }

    private static double dot(double[] a, double[] b) {
        return a[0] * b[0] + a[1] * b[1];
    }

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

        double[] thetaHat = {0.0, 0.0};
        double[] filteredPhi = {0.0, 0.0};
        double filteredY = 0.0;
        double integralDeltaSquared = 0.0;

        Path output = Path.of("Chapter21_Lesson4_results.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(
                output, StandardCharsets.UTF_8)) {
            writer.write(
                "time,phi1,phi2,y,phi1_filtered,phi2_filtered,y_filtered,"
                + "delta,mixed1,mixed2,theta1_hat,theta2_hat,theta1_error,"
                + "theta2_error,integral_delta_squared\n"
            );

            int steps = (int) Math.round(FINAL_TIME / DT);
            for (int step = 0; step <= steps; step++) {
                double time = step * DT;
                double[] phi = regressor(time);
                double y = dot(phi, THETA_TRUE);

                filteredY += DT * (-FILTER_RATE * filteredY + FILTER_RATE * y);
                for (int i = 0; i < filteredPhi.length; i++) {
                    filteredPhi[i] += DT * (
                        -FILTER_RATE * filteredPhi[i] + FILTER_RATE * phi[i]
                    );
                }

                double delta =
                    phi[0] * filteredPhi[1] - phi[1] * filteredPhi[0];
                double[] mixed = {
                    filteredPhi[1] * y - phi[1] * filteredY,
                    -filteredPhi[0] * y + phi[0] * filteredY
                };

                for (int i = 0; i < thetaHat.length; i++) {
                    thetaHat[i] += DT * GAMMA[i] * delta
                        * (mixed[i] - delta * thetaHat[i]);
                }
                integralDeltaSquared += DT * delta * delta;

                writer.write(String.format(
                    Locale.ROOT,
                    "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,"
                    + "%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f%n",
                    time, phi[0], phi[1], y, filteredPhi[0], filteredPhi[1],
                    filteredY, delta, mixed[0], mixed[1], thetaHat[0],
                    thetaHat[1], THETA_TRUE[0] - thetaHat[0],
                    THETA_TRUE[1] - thetaHat[1], integralDeltaSquared
                ));
            }
        }

        System.out.printf(Locale.ROOT,
            "DREM two-parameter simulation%n"
            + "true theta      = [%.10f, %.10f]%n"
            + "estimated theta = [%.10f, %.10f]%n"
            + "final error     = [%.10e, %.10e]%n"
            + "results written to %s%n",
            THETA_TRUE[0], THETA_TRUE[1], thetaHat[0], thetaHat[1],
            THETA_TRUE[0] - thetaHat[0], THETA_TRUE[1] - thetaHat[1], output
        );
    }
}

Chapter21_Lesson4.m

% Chapter 21, Lesson 4: Dynamic Regressor Extension and Mixing
% Two-parameter DREM simulation using a first-order stable filter.

clear; clc;

dt = 0.002;
T = 25.0;
lambda = 1.5;
gamma = [30.0; 30.0];
theta_true = [2.0; -1.5];

time = (0:dt:T).';
number_of_samples = numel(time);
theta_hat = zeros(2, 1);
filtered_phi = zeros(2, 1);
filtered_y = 0.0;
integral_delta_squared = 0.0;

history = zeros(number_of_samples, 15);

for k = 1:number_of_samples
    t = time(k);
    phi = [1.0 + 0.5*sin(0.7*t); ...
           0.7*cos(1.3*t) + 0.2*sin(2.4*t)];
    y = phi.' * theta_true;

    % H(s) = lambda / (s + lambda), implemented by forward Euler.
    filtered_y = filtered_y + dt*(-lambda*filtered_y + lambda*y);
    filtered_phi = filtered_phi + dt*(-lambda*filtered_phi + lambda*phi);

    % Extended regressor Phi_e = [phi.'; filtered_phi.'].
    delta = phi(1)*filtered_phi(2) - phi(2)*filtered_phi(1);

    % mathcalY = adj(Phi_e) * Y_e.
    mixed = [filtered_phi(2)*y - phi(2)*filtered_y; ...
             -filtered_phi(1)*y + phi(1)*filtered_y];

    % Decoupled scalar DREM estimators.
    theta_hat = theta_hat + dt*gamma.*delta.*(mixed - delta*theta_hat);
    integral_delta_squared = integral_delta_squared + dt*delta^2;

    history(k, :) = [t, phi.', y, filtered_phi.', filtered_y, delta, ...
        mixed.', theta_hat.', (theta_true - theta_hat).', ...
        integral_delta_squared];
end

variable_names = {'time','phi1','phi2','y','phi1_filtered', ...
    'phi2_filtered','y_filtered','delta','mixed1','mixed2', ...
    'theta1_hat','theta2_hat','theta1_error','theta2_error', ...
    'integral_delta_squared'};
results = array2table(history, 'VariableNames', variable_names);
writetable(results, 'Chapter21_Lesson4_results.csv');

fprintf('DREM two-parameter simulation\n');
fprintf('true theta      = [%.10f, %.10f]\n', theta_true(1), theta_true(2));
fprintf('estimated theta = [%.10f, %.10f]\n', theta_hat(1), theta_hat(2));
fprintf('final error     = [%.3e, %.3e]\n', ...
    theta_true(1)-theta_hat(1), theta_true(2)-theta_hat(2));

figure('Name', 'Chapter21 Lesson4 DREM');
plot(time, history(:, 11), 'LineWidth', 1.2); hold on;
plot(time, history(:, 12), 'LineWidth', 1.2);
yline(theta_true(1), '--');
yline(theta_true(2), '--');
grid on;
xlabel('Time (s)');
ylabel('Parameter estimate');
legend('theta1 hat','theta2 hat','theta1 true','theta2 true', ...
    'Location','best');
title('DREM parameter estimates');

% Simulink implementation note:
% Use one State-Space or Transfer Fcn block per filtered signal with
% transfer function lambda/(s+lambda). Feed the original and filtered
% regressions into a MATLAB Function block that computes delta, mixed1,
% mixed2, and the two scalar update laws. Integrator blocks realize the
% parameter estimates.

Simulink block structure

In Simulink, implement one transfer-function block \( 1.5/(s+1.5) \) for the regression output and one for each regressor component. A MATLAB Function block can calculate the determinant and mixed outputs. Two Integrator blocks implement the scalar parameter laws. Add Scope or To Workspace blocks for \( \Delta \), parameter estimates, and residuals. Use saturation or projection around the integrator states if parameter bounds are known.

Chapter21_Lesson4.nb

Notebook[{
 Cell["Chapter 21, Lesson 4: Dynamic Regressor Extension and Mixing", "Title"],
 Cell["Two-parameter DREM simulation using a first-order stable filter.", "Text"],
 Cell[BoxData[
  RowBox[{"ClearAll", "[", "\"Global`*\"", "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"dt", "=", "0.002"}], ";",
   RowBox[{"tFinal", "=", "25.0"}], ";",
   RowBox[{"lambda", "=", "1.5"}], ";",
   RowBox[{"gamma", "=", RowBox[{"{", RowBox[{"30.0", ",", "30.0"}], "}"}]}], ";",
   RowBox[{"thetaTrue", "=", RowBox[{"{", RowBox[{"2.0", ",", RowBox[{"-", "1.5"}]}], "}"}]}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"phi", "[", "t_", "]", ":=", RowBox[{"{",
    RowBox[{
     RowBox[{"1.0", "+", RowBox[{"0.5", " ", RowBox[{"Sin", "[", RowBox[{"0.7", " ", "t"}], "]"}]}]}], ",",
     RowBox[{RowBox[{"0.7", " ", RowBox[{"Cos", "[", RowBox[{"1.3", " ", "t"}], "]"}]}], "+", RowBox[{"0.2", " ", RowBox[{"Sin", "[", RowBox[{"2.4", " ", "t"}], "]"}]}]}]
    }], "}"}]}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"thetaHat", "=", RowBox[{"{", RowBox[{"0.0", ",", "0.0"}], "}"}]}], ";",
   RowBox[{"filteredPhi", "=", RowBox[{"{", RowBox[{"0.0", ",", "0.0"}], "}"}]}], ";",
   RowBox[{"filteredY", "=", "0.0"}], ";",
   RowBox[{"integralDeltaSquared", "=", "0.0"}], ";",
   RowBox[{"history", "=", RowBox[{"Reap", "[",
    RowBox[{"Do", "[",
     RowBox[{
      RowBox[{
       RowBox[{"p", "=", RowBox[{"phi", "[", "t", "]"}]}], ";",
       RowBox[{"y", "=", RowBox[{"p", ".", "thetaTrue"}]}], ";",
       RowBox[{"filteredY", "=", RowBox[{"filteredY", "+", RowBox[{"dt", " ", RowBox[{"(", RowBox[{RowBox[{"-", "lambda", " ", "filteredY"}], "+", RowBox[{"lambda", " ", "y"}]}], ")"}]}]}]}], ";",
       RowBox[{"filteredPhi", "=", RowBox[{"filteredPhi", "+", RowBox[{"dt", " ", RowBox[{"(", RowBox[{RowBox[{"-", "lambda", " ", "filteredPhi"}], "+", RowBox[{"lambda", " ", "p"}]}], ")"}]}]}]}], ";",
       RowBox[{"delta", "=", RowBox[{RowBox[{RowBox[{"p", "[", RowBox[{"[", "1", "]"}], "]"}], " ", RowBox[{"filteredPhi", "[", RowBox[{"[", "2", "]"}], "]"}]}], "-", RowBox[{RowBox[{"p", "[", RowBox[{"[", "2", "]"}], "]"}], " ", RowBox[{"filteredPhi", "[", RowBox[{"[", "1", "]"}], "]"}]}]}]}], ";",
       RowBox[{"mixed", "=", RowBox[{"{",
        RowBox[{
         RowBox[{RowBox[{RowBox[{"filteredPhi", "[", RowBox[{"[", "2", "]"}], "]"}], " ", "y"}], "-", RowBox[{RowBox[{"p", "[", RowBox[{"[", "2", "]"}], "]"}], " ", "filteredY"}]}], ",",
         RowBox[{RowBox[{"-", RowBox[{"filteredPhi", "[", RowBox[{"[", "1", "]"}], "]"}]}], " ", "y", "+", RowBox[{RowBox[{"p", "[", RowBox[{"[", "1", "]"}], "]"}], " ", "filteredY"}]}]
        }], "}"}]}], ";",
       RowBox[{"thetaHat", "=", RowBox[{"thetaHat", "+", RowBox[{"dt", " ", "gamma", " ", "delta", " ", RowBox[{"(", RowBox[{"mixed", "-", RowBox[{"delta", " ", "thetaHat"}]}], ")"}]}]}]}], ";",
       RowBox[{"integralDeltaSquared", "=", RowBox[{"integralDeltaSquared", "+", RowBox[{"dt", " ", RowBox[{"delta", "^", "2"}]}]}]}], ";",
       RowBox[{"Sow", "[", RowBox[{"Join", "[", RowBox[{RowBox[{"{", "t", "}"}], ",", "p", ",", RowBox[{"{", RowBox[{"y", ",", "delta"}], "}"}], ",", "mixed", ",", "thetaHat", ",", RowBox[{"thetaTrue", "-", "thetaHat"}], ",", RowBox[{"{", "integralDeltaSquared", "}"}]}], "]"}], "]"}]
      }], ",", RowBox[{"{", RowBox[{"t", ",", "0.0", ",", "tFinal", ",", "dt"}], "}"}]}], "]"}], "]"}], "[", RowBox[{"[", RowBox[{"2", ",", "1"}], "]"}], "]"}]}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{"Print", "[", RowBox[{"\"True theta: \"", ",", "thetaTrue", ",", "\"   Estimated theta: \"", ",", "thetaHat"}], "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{"ListLinePlot", "[",
   RowBox[{
    RowBox[{"{",
     RowBox[{
      RowBox[{"history", "[", RowBox[{"[", RowBox[{"All", ",", RowBox[{"{", RowBox[{"1", ",", "8"}], "}"}]}], "]"}], "]"}], ",",
      RowBox[{"history", "[", RowBox[{"[", RowBox[{"All", ",", RowBox[{"{", RowBox[{"1", ",", "9"}], "}"}]}], "]"}], "]"}]
     }], "}"}], ",",
    RowBox[{"PlotLegends", "->", RowBox[{"{", RowBox[{"\"theta1 hat\"", ",", "\"theta2 hat\""}], "}"}]}], ",",
    RowBox[{"AxesLabel", "->", RowBox[{"{", RowBox[{"\"time\"", ",", "\"estimate\""}], "}"}]}], ",",
    RowBox[{"PlotRange", "->", "All"}]
   }], "]"}]], "Input"]
}, WindowTitle -> "Chapter21_Lesson4", StyleDefinitions -> "Default.nb"]

Problem 1 (Two-parameter mixing): Starting from \( y=\phi_1\theta_1+\phi_2\theta_2 \) and \( y_f=\phi_{1f}\theta_1+\phi_{2f}\theta_2 \), derive the scalar regressions without using a matrix inverse.

Solution: Form

\[ \mathbf{M}_e= \begin{bmatrix}\phi_1&\phi_2\\\phi_{1f}&\phi_{2f}\end{bmatrix}, \qquad \Delta=\phi_1\phi_{2f}-\phi_2\phi_{1f}. \]

Premultiplication by the adjugate gives

\[ \begin{bmatrix} \phi_{2f}&-\phi_2\\-\phi_{1f}&\phi_1 \end{bmatrix} \begin{bmatrix}y\\y_f\end{bmatrix} =\Delta\begin{bmatrix}\theta_1\\\theta_2\end{bmatrix}. \]

Therefore,

\[ \mathcal{Y}_1=\phi_{2f}y-\phi_2y_f=\Delta\theta_1, \qquad \mathcal{Y}_2=-\phi_{1f}y+\phi_1y_f=\Delta\theta_2. \]

Problem 2 (Exact scalar error solution): Solve \( \dot{\tilde\theta}_i=-\gamma_i\Delta^2(t)\tilde\theta_i \) and state the weakest direct integral condition that forces asymptotic convergence.

Solution: Integration gives

\[ \tilde\theta_i(t)= \exp\!\left[-\gamma_i\int_{0}^{t}\Delta^2(\tau)d\tau\right] \tilde\theta_i(0). \]

The exponential factor tends to zero exactly when the accumulated integral diverges. Thus the direct condition is \( \int_0^{\infty}\Delta^2d\tau=\infty \), or \( \Delta\notin\mathcal{L}_2 \).

Problem 3 (Finite excitation burst): Suppose \( \Delta(t)=0 \) for all \( t>t_1 \), but \( E=\int_0^{t_1}\Delta^2(\tau)d\tau>0 \). Find the final error of the basic DREM gradient estimator.

Solution:

\[ \tilde\theta_i(t) =e^{-\gamma_i E}\tilde\theta_i(0), \qquad t\geq t_1. \]

The excitation burst contracts the error but, for finite \( E \), does not make it identically zero. Exact finite-time recovery requires a modified memory-based or finite-time DREM construction, not the elementary law alone.

Problem 4 (Composite Lyapunov derivative): For the error model and composite law in Section 7, derive \( \dot V \) in the ideal mixed-regression case.

Solution: Differentiating the state term yields

\[ \frac{d}{dt}\left(\mathbf{e}^{\top}\mathbf{P}\mathbf{e}\right) =-\mathbf{e}^{\top}\mathbf{Q}\mathbf{e} +2\mathbf{e}^{\top}\mathbf{P}\mathbf{B} \boldsymbol{\psi}^{\top}\tilde{\boldsymbol{\theta}}. \]

The tracking part of \( \dot{\hat{\boldsymbol{\theta}}} \) cancels the cross term. Since \( \boldsymbol{\mathcal{Y}}-\Delta\hat{\boldsymbol{\theta}} =-\Delta\tilde{\boldsymbol{\theta}} \), the prediction part contributes \( -2\Delta^2\tilde{\boldsymbol{\theta}}^{\top} \mathbf{K}_{D}\tilde{\boldsymbol{\theta}} \). Hence

\[ \dot V=-\mathbf{e}^{\top}\mathbf{Q}\mathbf{e} -2\Delta^2\tilde{\boldsymbol{\theta}}^{\top} \mathbf{K}_{D}\tilde{\boldsymbol{\theta}}\leq0. \]

Problem 5 (Discrete non-overshoot condition): For the Euler update, derive a condition that keeps the scalar error multiplier nonnegative and no greater than one.

Solution: The ideal discrete error is

\[ \tilde\theta_i[k+1] =\left(1-h\gamma_i\Delta^2[k]\right)\tilde\theta_i[k]. \]

Requiring the multiplier to belong to \( [0,1] \) gives \( 0\leq h\gamma_i\Delta^2[k]\leq1 \). If only asymptotic stability of the frozen scalar recursion is required, the wider interval is \( 0<h\gamma_i\Delta^2[k]<2 \).

Problem 6 (Filter design diagnosis): A three-parameter regression uses the original row and two identical first-order filters with the same pole and initial conditions. What happens to the extended determinant?

Solution: The two filtered rows are identical. Therefore the square extended matrix has two equal rows and its determinant is identically zero. At least two dynamically distinct extension operators are needed in addition to the original regression.

DREM begins with a regression linear in unknown parameters. Stable dynamic operators generate additional regressions, which are stacked into a square extended regressor matrix. Multiplication by its adjugate creates scalar equations \( \mathcal{Y}_i=\Delta\theta_i \). The resulting scalar gradient laws have explicit error solutions, elementwise monotonicity in the ideal case, and asymptotic convergence when \( \Delta\notin\mathcal{L}_2 \). In composite adaptive control, the DREM term strengthens model-based parameter correction while the tracking-error term continues to enforce closed-loop behavior.

The method does not eliminate the need for informative data. Success depends on the original regression, extension operators, determinant energy, numerical conditioning, and treatment of filter transients and noise. These distinctions are essential before applying DREM to the mechanical-system example in Lesson 5.