Chapter 5: MRAC for First-Order Linear Systems – Basics
Lesson 2: Reference Model and Error Dynamics Derivation
This lesson derives the desired first-order reference model, the ideal controller matching conditions, and the exact tracking-error dynamics for a scalar uncertain plant. The result isolates the effect of controller parameter mismatch and supplies the error model needed by the MIT-rule and Lyapunov designs in the next two lessons.
1. Learning Objectives and Lesson Boundary
After completing this lesson, the student should be able to:
- select a stable first-order reference model from desired time-domain behavior;
- derive the closed-loop plant obtained from a linearly parameterized controller;
- calculate the ideal controller gains that exactly match the plant to the reference model;
- derive the scalar and regressor forms of the tracking-error dynamics;
- distinguish tracking error, controller-parameter error, and closed-loop stability; and
- verify the derivation numerically without yet implementing an adaptation law.
The plant structure was introduced in Lesson 1. This lesson keeps the controller gains fixed while deriving the equations that an adaptive law must later regulate. The MIT rule is introduced in Lesson 3, and a Lyapunov-based update law is introduced in Lesson 4.
2. Plant, Command, and Controller Parameterization
Consider the first-order scalar plant
\[ \dot y(t)=-a\,y(t)+b\,u(t),\qquad a>0,\qquad b\neq 0, \]
where \(y(t)\) is the measured plant output, \(u(t)\) is the control input, and \(a,b\) are unknown constants. We assume that the sign of \(b\) is known. The command \(r(t)\) is bounded and piecewise continuous.
Use the two-parameter controller
\[ u(t)=k_y(t)y(t)+k_r(t)r(t). \]
For compactness, define the controller-parameter vector and regressor vector
\[ \boldsymbol{\theta}(t)= \begin{bmatrix}k_y(t)\\k_r(t)\end{bmatrix}, \qquad \boldsymbol{\phi}(t)= \begin{bmatrix}y(t)\\r(t)\end{bmatrix}, \qquad u(t)=\boldsymbol{\theta}^{T}(t)\boldsymbol{\phi}(t). \]
The gain \(k_y\) changes the effective closed-loop pole, whereas \(k_r\) changes the command-to-output gain. Two adjustable gains are therefore available to match the two coefficients of a first-order reference model.
3. Selecting the First-Order Reference Model
The desired response is generated by the stable reference model
\[ \dot y_m(t)=-a_m y_m(t)+b_m r(t), \qquad a_m>0. \]
Its transfer function from command to model output is
\[ M(s)=\frac{Y_m(s)}{R(s)}=\frac{b_m}{s+a_m}. \]
The model pole is \(-a_m\); hence the model time constant is \(\tau_m=1/a_m\). For a first-order system, the approximate two-percent settling time is
\[ T_s\approx \frac{4}{a_m}. \]
The model DC gain is \(b_m/a_m\). A common unit-gain choice is \(b_m=a_m\), which gives \(y_m(\infty)=r_0\) for a constant command \(r(t)=r_0\).
flowchart TD
R["Command r"] --> C["Controller u = k_y y + k_r r"]
C --> P["Unknown plant: y_dot = -a y + b u"]
P --> Y["Plant output y"]
Y --> C
R --> M["Reference model: y_m_dot = -a_m y_m + b_m r"]
M --> YM["Desired output y_m"]
Y --> E["Tracking error e = y - y_m"]
YM --> E
4. Closed-Loop Plant and Model-Matching Conditions
Substituting the controller into the plant gives
\[ \begin{aligned} \dot y &=-a y+b(k_y y+k_r r)\\ &=(-a+b k_y)y+b k_r r. \end{aligned} \]
Exact model matching requires the plant vector field to equal the reference-model vector field for every admissible value of \(y\) and \(r\). Therefore the coefficients must satisfy
\[ -a+b k_y^{*}=-a_m, \qquad b k_r^{*}=b_m. \]
Solving these equations yields the ideal controller gains
\[ \boxed{ k_y^{*}=\frac{a-a_m}{b}, \qquad k_r^{*}=\frac{b_m}{b} }. \]
These gains are called ideal because they are the values that would produce exact matching if the unknown plant parameters were known. Their existence is the first-order instance of the matching condition introduced earlier in the course.
Proposition 1 (exact matching).
For \(b\neq 0\), the closed-loop plant equals the reference model for all \(y,r\) if and only if \(k_y=k_y^{*}\) and \(k_r=k_r^{*}\).
Proof.
Sufficiency follows by direct substitution. For necessity, equality of the two affine functions \((-a+b k_y)y+b k_r r\) and \(-a_m y+b_m r\) for all independent \(y\) and \(r\) requires equality of their coefficients. Since \(b\neq 0\), the resulting two algebraic equations have the unique solution shown above.
5. Tracking Error and Parameter Error Definitions
Define the tracking error using the convention
\[ e(t)=y(t)-y_m(t). \]
Define controller-parameter errors relative to the ideal gains:
\[ \widetilde{k}_y(t)=k_y(t)-k_y^{*}, \qquad \widetilde{k}_r(t)=k_r(t)-k_r^{*}. \]
In vector notation,
\[ \widetilde{\boldsymbol{\theta} }(t) =\boldsymbol{\theta}(t)-\boldsymbol{\theta}^{*} = \begin{bmatrix} \widetilde{k}_y(t)\\ \widetilde{k}_r(t) \end{bmatrix}. \]
Tracking error and parameter error are different objects: \(e\) measures output mismatch, while \(\widetilde{\boldsymbol{\theta} }\) measures controller mismatch. Zero tracking error does not automatically imply zero parameter error.
6. Exact Error-Dynamics Derivation
Differentiate the tracking error:
\[ \dot e=\dot y-\dot y_m. \]
Substitute the closed-loop plant and reference-model equations:
\[ \dot e =(-a+b k_y)y+b k_r r+a_m y_m-b_m r. \]
Since \(y_m=y-e\),
\[ \begin{aligned} \dot e &=(-a+b k_y)y+b k_r r+a_m(y-e)-b_m r\\ &=-a_m e+(-a+b k_y+a_m)y+(b k_r-b_m)r. \end{aligned} \]
Add and subtract the ideal gains through \(k_y=k_y^{*}+\widetilde{k}_y\) and \(k_r=k_r^{*}+\widetilde{k}_r\). The matching conditions imply
\[ -a+b k_y+a_m=b\widetilde{k}_y, \qquad b k_r-b_m=b\widetilde{k}_r. \]
Hence the exact scalar error dynamics are
\[ \boxed{ \dot e(t) =-a_m e(t) +b\,\widetilde{k}_y(t)y(t) +b\,\widetilde{k}_r(t)r(t) }. \]
Equivalently, using the regressor vector,
\[ \boxed{ \dot e(t) =-a_m e(t) +b\,\widetilde{\boldsymbol{\theta} }^{T}(t) \boldsymbol{\phi}(t) }. \]
flowchart TD
A["Plant and controller equations"] --> B["Form e_dot = y_dot - y_m_dot"]
B --> C["Replace y_m by y - e"]
C --> D["Collect terms in e, y, and r"]
D --> F["Insert k_y = k_y_star + delta_k_y"]
F --> G["Insert k_r = k_r_star + delta_k_r"]
G --> H["Apply model-matching coefficient identities"]
H --> I["Obtain e_dot = -a_m e + b delta_theta_T phi"]
7. Structural Interpretation of the Error Equation
The error model separates into a stable nominal part and a mismatch forcing term:
\[ \underbrace{\dot e=-a_m e}_{\text{desired stable error dynamics} } \;+\; \underbrace{b\, \widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} } _{\text{forcing caused by controller mismatch} }. \]
Because \(a_m>0\), the homogeneous error system is exponentially stable. If the controller reaches the ideal gains, then
\[ \widetilde{\boldsymbol{\theta} }=\mathbf{0} \quad\text{implies}\quad \dot e=-a_m e, \]
and therefore
\[ e(t)=e(t_0)e^{-a_m(t-t_0)}. \]
Thus exact parameter matching guarantees exponential decay of tracking error even when the plant and model begin from different initial conditions.
When the parameter error is nonzero, the term \(\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}\) acts as an endogenous disturbance to the stable error filter \(1/(s+a_m)\). Future adaptive laws will adjust the gains so that this forcing term is controlled in a stability-compatible manner.
8. Integral Solution and an Input-to-Error Bound
The variation-of-constants formula gives the exact solution
\[ e(t)=e^{-a_m(t-t_0)}e(t_0) +b\int_{t_0}^{t} e^{-a_m(t-\tau)} \widetilde{\boldsymbol{\theta} }^{T}(\tau) \boldsymbol{\phi}(\tau)\,d\tau. \]
Suppose the mismatch product is bounded:
\[ \left| \widetilde{\boldsymbol{\theta} }^{T}(t) \boldsymbol{\phi}(t) \right|\le \bar d. \]
Taking absolute values and integrating the exponential kernel yields
\[ |e(t)| \le e^{-a_m(t-t_0)}|e(t_0)| +\frac{|b|\bar d}{a_m} \left(1-e^{-a_m(t-t_0)}\right). \]
Consequently,
\[ \limsup_{t→\infty}|e(t)| \le \frac{|b|\bar d}{a_m}. \]
This is not yet an adaptive-control convergence proof. It only shows how bounded mismatch forcing propagates through the stable reference error dynamics.
9. Closed-Loop Stability with Frozen Gains
Before adaptation is introduced, consider constant gains. The closed-loop plant is
\[ \dot y=-(a-b k_y)y+b k_r r. \]
Its pole is \(-(a-b k_y)\). Therefore frozen-gain internal stability requires
\[ a-b k_y>0. \]
Using \(k_y=k_y^{*}+\widetilde{k}_y\) and the matching relation, this condition becomes
\[ a_m-b\widetilde{k}_y>0. \]
A small tracking error at one instant does not guarantee this stability condition. Likewise, choosing \(k_r\) correctly cannot compensate for an unstable pole created by an unsuitable \(k_y\).
10. Transfer-Function and Steady-State Error Analysis
With constant gains and zero initial conditions,
\[ \frac{Y(s)}{R(s)} =\frac{b k_r}{s+a-b k_y}, \qquad \frac{Y_m(s)}{R(s)} =\frac{b_m}{s+a_m}. \]
Hence
\[ \frac{E(s)}{R(s)} = \frac{b k_r}{s+a-b k_y} - \frac{b_m}{s+a_m}. \]
Under the ideal gains, the two transfer functions are identical and \(E(s)/R(s)=0\). For a constant command \(r(t)=r_0\) and a stable frozen closed loop,
\[ e_{\mathrm{ss} } = \left( \frac{b k_r}{a-b k_y} - \frac{b_m}{a_m} \right)r_0. \]
Expressing the gains through their errors gives the exact relation
\[ e_{\mathrm{ss} } = \left( \frac{b_m+b\widetilde{k}_r} {a_m-b\widetilde{k}_y} - \frac{b_m}{a_m} \right)r_0. \]
For small parameter errors, first-order expansion around the ideal gains gives
\[ e_{\mathrm{ss} } \approx \left( \frac{b}{a_m}\widetilde{k}_r + \frac{b b_m}{a_m^2}\widetilde{k}_y \right)r_0. \]
Thus both feedback-gain error and command-gain error contribute to steady-state tracking error.
11. Sign Convention and Equivalent Forms
Some texts define the tracking error as \(e_a=y_m-y\). Since \(e_a=-e\), the corresponding error equation is
\[ \dot e_a =-a_m e_a -b\,\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}. \]
Both conventions are valid. However, the sign must remain consistent when deriving gradient or Lyapunov update laws. A copied adaptive law combined with the opposite error convention may reverse the stabilizing direction of parameter adjustment.
The plant is sometimes written as \(\tau\dot y+y=K u\). Dividing by \(\tau\) gives
\[ \dot y=-\frac{1}{\tau}y+\frac{K}{\tau}u, \qquad a=\frac{1}{\tau}, \qquad b=\frac{K}{\tau}. \]
Therefore the derivation is unchanged after identifying the equivalent \(a\) and \(b\).
12. Numerical Verification Strategy
The implementations below use \(a=1.4\), \(b=2.0\), \(a_m=b_m=2.5\). The ideal gains are
\[ k_y^{*}=-0.55, \qquad k_r^{*}=1.25. \]
To expose the forcing term, the simulations deliberately freeze the gains at \(k_y=k_y^{*}+0.18=-0.37\) and \(k_r=k_r^{*}-0.15=1.10\). At every simulated sample, the programs compare
\[ \dot y-\dot y_m \]
with the independently assembled expression
\[ -a_m e+b\widetilde{k}_y y+b\widetilde{k}_r r. \]
Their difference should be at floating-point roundoff level. This test verifies the algebraic identity; it is not a proof of parameter adaptation or parameter convergence.
Relevant software ecosystems include NumPy, SciPy, Matplotlib, and
python-control in Python; Eigen and Boost.Odeint in C++;
Apache Commons Math in Java; MATLAB, Control System Toolbox, and
Simulink; and NDSolve, StateSpaceModel, and
symbolic simplification in Wolfram Mathematica. The examples use
from-scratch RK4 integration where practical so that the numerical
mechanism remains visible.
13. Python Implementation
This implementation uses NumPy and Matplotlib and performs fixed-step RK4 integration from scratch. It saves the simulated signals and the error-identity residual.
Chapter5_Lesson2.py
"""
Chapter5_Lesson2.py
Reference-model and error-dynamics verification for a first-order MRAC setup.
The controller gains are deliberately frozen. No adaptive update law is used in
this lesson; the purpose is to verify the matching conditions and the exact
tracking-error differential equation.
Dependencies:
numpy
matplotlib
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import Callable
import matplotlib.pyplot as plt
import numpy as np
@dataclass(frozen=True)
class Parameters:
a: float = 1.4
b: float = 2.0
a_m: float = 2.5
b_m: float = 2.5
delta_k_y: float = 0.18
delta_k_r: float = -0.15
@property
def k_y_star(self) -> float:
return (self.a - self.a_m) / self.b
@property
def k_r_star(self) -> float:
return self.b_m / self.b
@property
def k_y(self) -> float:
return self.k_y_star + self.delta_k_y
@property
def k_r(self) -> float:
return self.k_r_star + self.delta_k_r
def reference(t: float) -> float:
"""Bounded piecewise-constant command."""
if t < 4.0:
return 1.0
if t < 8.0:
return -0.5
return 0.75
def dynamics(t: float, state: np.ndarray, p: Parameters) -> np.ndarray:
"""Return [plant derivative, reference-model derivative]."""
y, y_m = state
r = reference(t)
u = p.k_y * y + p.k_r * r
y_dot = -p.a * y + p.b * u
y_m_dot = -p.a_m * y_m + p.b_m * r
return np.array([y_dot, y_m_dot], dtype=float)
def rk4_step(
f: Callable[[float, np.ndarray, Parameters], np.ndarray],
t: float,
state: np.ndarray,
h: float,
p: Parameters,
) -> np.ndarray:
"""One classical fourth-order Runge-Kutta step."""
k1 = f(t, state, p)
k2 = f(t + 0.5 * h, state + 0.5 * h * k1, p)
k3 = f(t + 0.5 * h, state + 0.5 * h * k2, p)
k4 = f(t + h, state + h * k3, p)
return state + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
def simulate(p: Parameters, t_final: float = 12.0, h: float = 0.002):
n_steps = int(round(t_final / h))
time = np.linspace(0.0, t_final, n_steps + 1)
state = np.zeros((n_steps + 1, 2), dtype=float)
state[0] = np.array([0.4, -0.2], dtype=float)
for i in range(n_steps):
state[i + 1] = rk4_step(dynamics, time[i], state[i], h, p)
y = state[:, 0]
y_m = state[:, 1]
r = np.array([reference(t) for t in time])
e = y - y_m
plant_model_difference = (
-p.a * y + p.b * (p.k_y * y + p.k_r * r)
- (-p.a_m * y_m + p.b_m * r)
)
derived_error_rhs = (
-p.a_m * e
+ p.b * (p.delta_k_y * y + p.delta_k_r * r)
)
identity_residual = plant_model_difference - derived_error_rhs
return time, r, y, y_m, e, identity_residual
def main() -> None:
p = Parameters()
time, r, y, y_m, e, residual = simulate(p)
print(f"k_y* = {p.k_y_star:.6f}")
print(f"k_r* = {p.k_r_star:.6f}")
print(f"frozen k_y = {p.k_y:.6f}")
print(f"frozen k_r = {p.k_r:.6f}")
print(f"max error-dynamics identity residual = {np.max(np.abs(residual)):.3e}")
data = np.column_stack((time, r, y, y_m, e, residual))
np.savetxt(
"Chapter5_Lesson2_results.csv",
data,
delimiter=",",
header="time,r,y,y_m,e,identity_residual",
comments="",
)
plt.figure(figsize=(9, 5))
plt.plot(time, r, "--", label="reference r")
plt.plot(time, y, label="plant output y")
plt.plot(time, y_m, label="model output y_m")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.title("First-Order Plant and Reference Model")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig("Chapter5_Lesson2_outputs.png", dpi=180)
plt.figure(figsize=(9, 4))
plt.plot(time, e, label="tracking error e = y - y_m")
plt.xlabel("Time (s)")
plt.ylabel("Error")
plt.title("Tracking Error with Frozen Nonideal Gains")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig("Chapter5_Lesson2_error.png", dpi=180)
plt.show()
if __name__ == "__main__":
main()
14. C++ Implementation
This dependency-free C++17 implementation performs RK4 integration, checks the derived identity, and writes a CSV file. Eigen or Boost.Odeint could replace the small fixed-size routines in larger projects.
Chapter5_Lesson2.cpp
/*
Chapter5_Lesson2.cpp
Reference-model and error-dynamics verification for a first-order MRAC setup.
Build:
g++ -std=c++17 -O2 Chapter5_Lesson2.cpp -o Chapter5_Lesson2
Run:
./Chapter5_Lesson2
*/
#include <algorithm>
#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>
struct Parameters {
double a{1.4};
double b{2.0};
double a_m{2.5};
double b_m{2.5};
double delta_k_y{0.18};
double delta_k_r{-0.15};
[[nodiscard]] double k_y_star() const { return (a - a_m) / b; }
[[nodiscard]] double k_r_star() const { return b_m / b; }
[[nodiscard]] double k_y() const { return k_y_star() + delta_k_y; }
[[nodiscard]] double k_r() const { return k_r_star() + delta_k_r; }
};
using State = std::array<double, 2>;
double reference(const double t) {
if (t < 4.0) {
return 1.0;
}
if (t < 8.0) {
return -0.5;
}
return 0.75;
}
State dynamics(const double t, const State& x, const Parameters& p) {
const double y = x[0];
const double y_m = x[1];
const double r = reference(t);
const double u = p.k_y() * y + p.k_r() * r;
return {
-p.a * y + p.b * u,
-p.a_m * y_m + p.b_m * r
};
}
State add_scaled(const State& x, const State& k, const double scale) {
return {x[0] + scale * k[0], x[1] + scale * k[1]};
}
State rk4_step(
const double t,
const State& x,
const double h,
const Parameters& p
) {
const State k1 = dynamics(t, x, p);
const State k2 = dynamics(t + 0.5 * h, add_scaled(x, k1, 0.5 * h), p);
const State k3 = dynamics(t + 0.5 * h, add_scaled(x, k2, 0.5 * h), p);
const State k4 = dynamics(t + h, add_scaled(x, k3, h), p);
return {
x[0] + (h / 6.0) * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]),
x[1] + (h / 6.0) * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1])
};
}
int main() {
try {
const Parameters p;
constexpr double h = 0.002;
constexpr double t_final = 12.0;
const int n_steps = static_cast<int>(std::lround(t_final / h));
State x{0.4, -0.2};
double max_residual = 0.0;
std::ofstream csv("Chapter5_Lesson2_results_cpp.csv");
if (!csv) {
throw std::runtime_error("Cannot open the output CSV file.");
}
csv << "time,r,y,y_m,e,identity_residual\n";
csv << std::setprecision(12);
for (int i = 0; i <= n_steps; ++i) {
const double t = i * h;
const double r = reference(t);
const double y = x[0];
const double y_m = x[1];
const double e = y - y_m;
const double plant_model_difference =
(-p.a * y + p.b * (p.k_y() * y + p.k_r() * r))
- (-p.a_m * y_m + p.b_m * r);
const double derived_error_rhs =
-p.a_m * e
+ p.b * (p.delta_k_y * y + p.delta_k_r * r);
const double residual =
plant_model_difference - derived_error_rhs;
max_residual = std::max(max_residual, std::abs(residual));
csv << t << ',' << r << ',' << y << ',' << y_m << ','
<< e << ',' << residual << '\n';
if (i < n_steps) {
x = rk4_step(t, x, h, p);
}
}
std::cout << std::fixed << std::setprecision(6);
std::cout << "k_y* = " << p.k_y_star() << '\n';
std::cout << "k_r* = " << p.k_r_star() << '\n';
std::cout << "frozen k_y = " << p.k_y() << '\n';
std::cout << "frozen k_r = " << p.k_r() << '\n';
std::cout << std::scientific;
std::cout << "max error-dynamics identity residual = "
<< max_residual << '\n';
return 0;
} catch (const std::exception& ex) {
std::cerr << "Error: " << ex.what() << '\n';
return 1;
}
}
15. Java Implementation
This Java implementation uses only the standard library. Apache Commons Math provides reusable ODE solvers for larger adaptive-control simulations, while libraries such as XChart can be added for plotting.
Chapter5_Lesson2.java
/*
Chapter5_Lesson2.java
Reference-model and error-dynamics verification for a first-order MRAC setup.
Build:
javac Chapter5_Lesson2.java
Run:
java Chapter5_Lesson2
*/
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 Chapter5_Lesson2 {
private Chapter5_Lesson2() {
}
private static final class Parameters {
final double a = 1.4;
final double b = 2.0;
final double aM = 2.5;
final double bM = 2.5;
final double deltaKY = 0.18;
final double deltaKR = -0.15;
double kYStar() {
return (a - aM) / b;
}
double kRStar() {
return bM / b;
}
double kY() {
return kYStar() + deltaKY;
}
double kR() {
return kRStar() + deltaKR;
}
}
private static double reference(double t) {
if (t < 4.0) {
return 1.0;
}
if (t < 8.0) {
return -0.5;
}
return 0.75;
}
private static double[] dynamics(
double t,
double[] state,
Parameters p
) {
double y = state[0];
double yM = state[1];
double r = reference(t);
double u = p.kY() * y + p.kR() * r;
return new double[] {
-p.a * y + p.b * u,
-p.aM * yM + p.bM * r
};
}
private static double[] addScaled(
double[] state,
double[] slope,
double scale
) {
return new double[] {
state[0] + scale * slope[0],
state[1] + scale * slope[1]
};
}
private static double[] rk4Step(
double t,
double[] state,
double h,
Parameters p
) {
double[] k1 = dynamics(t, state, p);
double[] k2 = dynamics(
t + 0.5 * h,
addScaled(state, k1, 0.5 * h),
p
);
double[] k3 = dynamics(
t + 0.5 * h,
addScaled(state, k2, 0.5 * h),
p
);
double[] k4 = dynamics(
t + h,
addScaled(state, k3, h),
p
);
return new double[] {
state[0] + (h / 6.0)
* (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]),
state[1] + (h / 6.0)
* (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1])
};
}
public static void main(String[] args) {
Locale.setDefault(Locale.US);
Parameters p = new Parameters();
final double h = 0.002;
final double tFinal = 12.0;
final int nSteps = (int) Math.round(tFinal / h);
double[] state = {0.4, -0.2};
double maxResidual = 0.0;
Path output = Path.of("Chapter5_Lesson2_results_java.csv");
try (BufferedWriter writer = Files.newBufferedWriter(
output,
StandardCharsets.UTF_8
)) {
writer.write("time,r,y,y_m,e,identity_residual");
writer.newLine();
for (int i = 0; i <= nSteps; i++) {
double t = i * h;
double r = reference(t);
double y = state[0];
double yM = state[1];
double e = y - yM;
double plantModelDifference =
(-p.a * y + p.b * (p.kY() * y + p.kR() * r))
- (-p.aM * yM + p.bM * r);
double derivedErrorRhs =
-p.aM * e
+ p.b * (p.deltaKY * y + p.deltaKR * r);
double residual =
plantModelDifference - derivedErrorRhs;
maxResidual = Math.max(maxResidual, Math.abs(residual));
writer.write(String.format(
Locale.US,
"%.12f,%.12f,%.12f,%.12f,%.12f,%.12e",
t, r, y, yM, e, residual
));
writer.newLine();
if (i < nSteps) {
state = rk4Step(t, state, h, p);
}
}
} catch (IOException ex) {
System.err.println("Failed to write results: " + ex.getMessage());
System.exit(1);
}
System.out.printf(Locale.US, "k_y* = %.6f%n", p.kYStar());
System.out.printf(Locale.US, "k_r* = %.6f%n", p.kRStar());
System.out.printf(Locale.US, "frozen k_y = %.6f%n", p.kY());
System.out.printf(Locale.US, "frozen k_r = %.6f%n", p.kR());
System.out.printf(
Locale.US,
"max error-dynamics identity residual = %.3e%n",
maxResidual
);
}
}
16. MATLAB and Simulink Implementation
The script performs a from-scratch RK4 simulation, optionally checks the frozen-gain transfer function with Control System Toolbox, and programmatically creates a Simulink model when Simulink is available.
Chapter5_Lesson2.m
% Chapter5_Lesson2.m
% Reference-model and error-dynamics verification for a first-order MRAC setup.
%
% The controller gains are frozen. This lesson does not implement an adaptive
% update law. It verifies the matching conditions and exact error dynamics.
%
% Related MATLAB products:
% MATLAB, Control System Toolbox, Simulink
clear; clc; close all;
p.a = 1.4;
p.b = 2.0;
p.a_m = 2.5;
p.b_m = 2.5;
p.delta_k_y = 0.18;
p.delta_k_r = -0.15;
p.k_y_star = (p.a - p.a_m) / p.b;
p.k_r_star = p.b_m / p.b;
p.k_y = p.k_y_star + p.delta_k_y;
p.k_r = p.k_r_star + p.delta_k_r;
h = 0.002;
t = (0:h:12).';
n = numel(t);
x = zeros(n, 2);
x(1, :) = [0.4, -0.2];
for i = 1:n-1
x(i + 1, :) = rk4_step(t(i), x(i, :).', h, p).';
end
r = arrayfun(@reference_signal, t);
y = x(:, 1);
y_m = x(:, 2);
e = y - y_m;
plant_model_difference = ...
(-p.a .* y + p.b .* (p.k_y .* y + p.k_r .* r)) ...
- (-p.a_m .* y_m + p.b_m .* r);
derived_error_rhs = ...
-p.a_m .* e ...
+ p.b .* (p.delta_k_y .* y + p.delta_k_r .* r);
identity_residual = plant_model_difference - derived_error_rhs;
fprintf('k_y* = %.6f\n', p.k_y_star);
fprintf('k_r* = %.6f\n', p.k_r_star);
fprintf('frozen k_y = %.6f\n', p.k_y);
fprintf('frozen k_r = %.6f\n', p.k_r);
fprintf('max error-dynamics identity residual = %.3e\n', ...
max(abs(identity_residual)));
results = table(t, r, y, y_m, e, identity_residual);
writetable(results, 'Chapter5_Lesson2_results_matlab.csv');
figure;
plot(t, r, '--', t, y, t, y_m, 'LineWidth', 1.2);
grid on;
xlabel('Time (s)');
ylabel('Amplitude');
title('First-Order Plant and Reference Model');
legend('reference r', 'plant output y', 'model output y_m', ...
'Location', 'best');
figure;
plot(t, e, 'LineWidth', 1.2);
grid on;
xlabel('Time (s)');
ylabel('Error');
title('Tracking Error with Frozen Nonideal Gains');
legend('e = y - y_m', 'Location', 'best');
% Optional transfer-function check with Control System Toolbox.
if license('test', 'Control_Toolbox')
plant_closed_loop = tf(p.b * p.k_r, [1, p.a - p.b * p.k_y]);
reference_model = tf(p.b_m, [1, p.a_m]);
error_transfer = minreal(plant_closed_loop - reference_model);
disp('Frozen-gain transfer function R(s) -> E(s):');
disp(error_transfer);
end
% Optional programmatic Simulink construction.
if license('test', 'Simulink')
model = 'Chapter5_Lesson2_Simulink';
if bdIsLoaded(model)
close_system(model, 0);
end
if exist([model '.slx'], 'file')
delete([model '.slx']);
end
new_system(model);
open_system(model);
add_block('simulink/Sources/Step', [model '/Reference'], ...
'Time', '0', 'Before', '0', 'After', '1', ...
'Position', [30 80 60 110]);
add_block('simulink/Math Operations/Gain', [model '/k_r'], ...
'Gain', 'p.k_r', 'Position', [110 55 165 90]);
add_block('simulink/Math Operations/Gain', [model '/k_y'], ...
'Gain', 'p.k_y', 'Position', [110 145 165 180]);
add_block('simulink/Math Operations/Sum', [model '/ControlSum'], ...
'Inputs', '++', 'Position', [215 90 240 145]);
add_block('simulink/Continuous/Transfer Fcn', [model '/Plant'], ...
'Numerator', 'p.b', 'Denominator', '[1 p.a]', ...
'Position', [290 90 390 140]);
add_block('simulink/Continuous/Transfer Fcn', ...
[model '/ReferenceModel'], ...
'Numerator', 'p.b_m', 'Denominator', '[1 p.a_m]', ...
'Position', [290 220 390 270]);
add_block('simulink/Math Operations/Sum', [model '/Error'], ...
'Inputs', '+-', 'Position', [450 130 475 185]);
add_block('simulink/Sinks/Scope', [model '/Scope'], ...
'NumInputPorts', '3', 'Position', [545 105 580 205]);
add_line(model, 'Reference/1', 'k_r/1', 'autorouting', 'on');
add_line(model, 'Reference/1', 'ReferenceModel/1', 'autorouting', 'on');
add_line(model, 'k_r/1', 'ControlSum/1', 'autorouting', 'on');
add_line(model, 'Plant/1', 'k_y/1', 'autorouting', 'on');
add_line(model, 'k_y/1', 'ControlSum/2', 'autorouting', 'on');
add_line(model, 'ControlSum/1', 'Plant/1', 'autorouting', 'on');
add_line(model, 'Plant/1', 'Error/1', 'autorouting', 'on');
add_line(model, 'ReferenceModel/1', 'Error/2', 'autorouting', 'on');
add_line(model, 'Plant/1', 'Scope/1', 'autorouting', 'on');
add_line(model, 'ReferenceModel/1', 'Scope/2', 'autorouting', 'on');
add_line(model, 'Error/1', 'Scope/3', 'autorouting', 'on');
set_param(model, 'StopTime', '12', 'Solver', 'ode45');
save_system(model);
fprintf('Created %s.slx\n', model);
end
function r = reference_signal(t)
if t < 4
r = 1.0;
elseif t < 8
r = -0.5;
else
r = 0.75;
end
end
function dx = system_dynamics(t, x, p)
r = reference_signal(t);
y = x(1);
y_m = x(2);
u = p.k_y * y + p.k_r * r;
dx = [-p.a * y + p.b * u; ...
-p.a_m * y_m + p.b_m * r];
end
function x_next = rk4_step(t, x, h, p)
k1 = system_dynamics(t, x, p);
k2 = system_dynamics(t + h / 2, x + h * k1 / 2, p);
k3 = system_dynamics(t + h / 2, x + h * k2 / 2, p);
k4 = system_dynamics(t + h, x + h * k3, p);
x_next = x + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6;
end
17. Wolfram Mathematica Implementation
The notebook uses NDSolveValue for simulation, verifies the differential identity numerically, exports results, and performs a symbolic simplification of the error equation.
Chapter5_Lesson2.nb
Notebook[{
Cell["Chapter 5, Lesson 2: Reference Model and Error Dynamics Derivation", "Title"],
Cell["Frozen-gain verification; no adaptive update law is implemented.", "Text"],
Cell["ClearAll[\"Global`*\"];\na = 1.4; b = 2.0; am = 2.5; bm = 2.5;\ndeltaKy = 0.18; deltaKr = -0.15;\nkyStar = (a - am)/b;\nkrStar = bm/b;\nky = kyStar + deltaKy;\nkr = krStar + deltaKr;", "Input"],
Cell["r[t_?NumericQ] := Piecewise[{ {1.0, t < 4.0}, {-0.5, t < 8.0} }, 0.75];", "Input"],
Cell["solution = NDSolveValue[\n {\n y'[t] == -a y[t] + b (ky y[t] + kr r[t]),\n ym'[t] == -am ym[t] + bm r[t],\n y[0] == 0.4,\n ym[0] == -0.2\n },\n {y, ym},\n {t, 0, 12},\n MaxStepFraction -> 1/2000\n];\nySol[t_] := solution[[1]][t];\nymSol[t_] := solution[[2]][t];\neSol[t_] := ySol[t] - ymSol[t];", "Input"],
Cell["plantModelDifference[t_?NumericQ] :=\n (-a ySol[t] + b (ky ySol[t] + kr r[t])) -\n (-am ymSol[t] + bm r[t]);\nderivedErrorRHS[t_?NumericQ] :=\n -am eSol[t] + b (deltaKy ySol[t] + deltaKr r[t]);\nresidual[t_?NumericQ] :=\n plantModelDifference[t] - derivedErrorRHS[t];", "Input"],
Cell["sampleTimes = Subdivide[0, 12, 6000];\nmaxResidual = Max[Abs[residual /@ sampleTimes]];\nPrint[\"ky* = \", N[kyStar, 8]];\nPrint[\"kr* = \", N[krStar, 8]];\nPrint[\"frozen ky = \", N[ky, 8]];\nPrint[\"frozen kr = \", N[kr, 8]];\nPrint[\"max error-dynamics identity residual = \", ScientificForm[maxResidual]];", "Input"],
Cell["outputPlot = Plot[\n Evaluate[{r[t], ySol[t], ymSol[t]}],\n {t, 0, 12},\n PlotLegends -> {\"reference r\", \"plant output y\", \"model output ym\"},\n AxesLabel -> {\"Time (s)\", \"Amplitude\"},\n PlotLabel -> \"First-Order Plant and Reference Model\",\n PlotRange -> All,\n ImageSize -> Large\n];\nerrorPlot = Plot[\n eSol[t],\n {t, 0, 12},\n PlotLegends -> {\"e = y - ym\"},\n AxesLabel -> {\"Time (s)\", \"Error\"},\n PlotLabel -> \"Tracking Error with Frozen Nonideal Gains\",\n PlotRange -> All,\n ImageSize -> Large\n];\nColumn[{outputPlot, errorPlot}]", "Input"],
Cell["Export[\n \"Chapter5_Lesson2_results_mathematica.csv\",\n Prepend[\n Table[{tt, r[tt], ySol[tt], ymSol[tt], eSol[tt], residual[tt]},\n {tt, sampleTimes}],\n {\"time\", \"r\", \"y\", \"y_m\", \"e\", \"identity_residual\"}\n ]\n];", "Input"],
Cell["Symbolic derivation", "Section"],
Cell["errorEquation = FullSimplify[\n (-a y + b (ky y + kr rr)) - (-am ym + bm rr) /. {\n y -> e + ym,\n ky -> kyStar + deltaKy,\n kr -> krStar + deltaKr,\n kyStar -> (a - am)/b,\n krStar -> bm/b\n }\n];\nCollect[errorEquation, {e, deltaKy, deltaKr}]", "Input"]
},
WindowTitle -> "Chapter5_Lesson2",
StyleDefinitions -> "Default.nb"
]
18. Common Derivation and Interpretation Errors
- Changing the error convention midway. Defining \(e=y-y_m\) and later using a formula derived for \(e=y_m-y\) changes the sign of the mismatch term.
- Replacing \(y_m\) incorrectly. From \(e=y-y_m\), the correct identity is \(y_m=y-e\).
- Matching only the pole. Choosing \(k_y\) to match the pole does not generally match the command gain; \(k_r\) must also satisfy its matching equation.
- Assuming arbitrary gains preserve stability. For frozen gains, \(a-b k_y>0\) is required.
- Confusing tracking with identification. Tracking may be satisfactory even when the controller parameters do not equal the ideal values, particularly for insufficiently informative commands.
- Introducing an update law too early. This lesson establishes the error model. The selection and proof of an adaptive law require additional arguments developed in subsequent lessons.
19. Problems and Solutions
Problem 1 (Ideal matching gains): Consider \(\dot y=-1.4y+2u\) and the reference model \(\dot y_m=-2.5y_m+2.5r\). Find \(k_y^{*}\) and \(k_r^{*}\) for \(u=k_y y+k_r r\).
Solution:
\[ k_y^{*} =\frac{a-a_m}{b} =\frac{1.4-2.5}{2} =-0.55, \qquad k_r^{*} =\frac{b_m}{b} =\frac{2.5}{2} =1.25. \]
Substitution gives \(\dot y=-2.5y+2.5r\), exactly equal to the reference model vector field.
Problem 2 (Derivation with the opposite error convention): Define \(e_a=y_m-y\). Derive its dynamics using the same parameter-error convention \(\widetilde{\boldsymbol{\theta} } =\boldsymbol{\theta}-\boldsymbol{\theta}^{*}\).
Solution:
Since \(e_a=-e\), differentiation gives \(\dot e_a=-\dot e\). Using the previously derived equation,
\[ \dot e_a =-\left(-a_m e+ b\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}\right) =-a_m e_a -b\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}. \]
The nominal stable term retains the form \(-a_m e_a\), while the parameter-mismatch term changes sign.
Problem 3 (Frozen-gain stability region): For \(a=1.4\) and \(b=2\), determine the set of constant \(k_y\) values that stabilize the closed-loop plant. Also express the condition in terms of \(\widetilde{k}_y\) when \(a_m=2.5\).
Solution:
The stability condition is
\[ a-b k_y>0 \quad\text{is equivalent to}\quad 1.4-2k_y>0 \quad\text{is equivalent to}\quad k_y<0.7. \]
Since \(k_y=k_y^{*}+\widetilde{k}_y\) and \(k_y^{*}=-0.55\),
\[ -0.55+\widetilde{k}_y<0.7 \quad\text{is equivalent to}\quad \widetilde{k}_y<1.25. \]
This agrees with \(a_m-b\widetilde{k}_y>0\), namely \(2.5-2\widetilde{k}_y>0\).
Problem 4 (Steady-state error under mismatch): Use the numerical values in Section 12 and a unit step command. Calculate the exact frozen-gain steady-state tracking error.
Solution:
The frozen gains are \(k_y=-0.37\) and \(k_r=1.10\). Therefore
\[ \frac{y_{\mathrm{ss} } }{r_0} =\frac{b k_r}{a-b k_y} =\frac{2(1.10)}{1.4-2(-0.37)} =\frac{2.2}{2.14} \approx 1.028037. \]
Since \(b_m/a_m=1\),
\[ e_{\mathrm{ss} } \approx (1.028037-1)r_0 =0.028037 \]
for \(r_0=1\).
Problem 5 (Error bound): Suppose \(a_m=2\), \(|b|=3\), \(|e(0)|=0.5\), and \(|\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}| \le 0.1\). Give an upper bound for \(|e(t)|\) and its asymptotic limit superior.
Solution:
\[ |e(t)| \le 0.5e^{-2t} +\frac{3(0.1)}{2}(1-e^{-2t}) = 0.5e^{-2t}+0.15(1-e^{-2t}). \]
Therefore
\[ \limsup_{t→\infty}|e(t)|\le 0.15. \]
Problem 6 (Why zero tracking error need not identify the gains): Let \(r(t)=0\), \(y(0)=y_m(0)=0\). Show that \(e(t)=0\) for arbitrary constant \(k_r\), and explain the implication.
Solution:
With zero command and zero initial conditions, both model and plant remain at zero as long as the plant closed loop is well posed: \(y(t)=y_m(t)=0\). The regressor is \(\boldsymbol{\phi}=[0,0]^T\), so
\[ \dot e=-a_m e+ b\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} =-a_m e, \qquad e(0)=0, \]
and hence \(e(t)=0\). The value of \(k_r\) cannot be inferred from this experiment because the command channel is never excited. This anticipates the distinction between tracking convergence and parameter convergence.
20. Summary
For the uncertain first-order plant \(\dot y=-ay+bu\), the controller \(u=k_y y+k_r r\) can exactly reproduce the reference model \(\dot y_m=-a_m y_m+b_m r\) when
\[ k_y^{*}=\frac{a-a_m}{b}, \qquad k_r^{*}=\frac{b_m}{b}. \]
With \(e=y-y_m\) and parameter errors measured from these ideal gains, the exact tracking-error model is
\[ \dot e=-a_m e+ b\widetilde{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}, \qquad \boldsymbol{\phi}= \begin{bmatrix}y\\r\end{bmatrix}. \]
This equation exposes a stable nominal error filter driven by controller mismatch. It is the central starting point for the MIT-rule update law in Lesson 3 and the Lyapunov-based adaptive design in Lesson 4.
21. References
- Parks, P.C. (1966). Liapunov redesign of model reference adaptive control systems. IEEE Transactions on Automatic Control, 11(3), 362–367. doi:10.1109/TAC.1966.1098361.
- Monopoli, R.V. (1974). Model reference adaptive control with an augmented error signal. IEEE Transactions on Automatic Control, 19(5), 474–484.
- Narendra, K.S., & Valavani, L.S. (1978). Stable adaptive controller design—Direct control. IEEE Transactions on Automatic Control, 23(4), 570–583. doi:10.1109/TAC.1978.1101823.
- Feuer, A., & Morse, A.S. (1978). Adaptive control of single-input, single-output linear systems. IEEE Transactions on Automatic Control, 23(4), 557–569.
- Morse, A.S. (1980). Global stability of parameter-adaptive control systems. IEEE Transactions on Automatic Control, 25(3), 433–439. doi:10.1109/TAC.1980.1102364.
- Narendra, K.S., Lin, Y.-H., & Valavani, L.S. (1980). Stable adaptive controller design, Part II: Proof of stability. IEEE Transactions on Automatic Control, 25(3), 440–448. doi:10.1109/TAC.1980.1102362.
- Egardt, B. (1980). Stability analysis of continuous-time adaptive control systems. SIAM Journal on Control and Optimization, 18(5), 540–558. doi:10.1137/0318040.
- Gibson, T.E., Annaswamy, A.M., & Lavretsky, E. (2013). On adaptive control with closed-loop reference models: Transients, oscillations, and peaking. IEEE Access, 1, 703–717. doi:10.1109/ACCESS.2013.2284005.
- Serrani, A. (2013). An output regulation perspective on the model reference adaptive control problem. International Journal of Adaptive Control and Signal Processing, 27(1–2), 22–34. doi:10.1002/acs.2340.
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.
