Chapter 7: MRAC for General SISO Systems in Canonical Forms
Lesson 2: Control Law Structure and Parameterization
This lesson constructs the adjustable state-feedback and command-feedforward controller used in canonical-form model reference adaptive control. The main objective is to express the control input as a scalar product between an adjustable parameter vector and a measurable regressor, derive the exact model-matching conditions, and obtain the error equation that will support the Lyapunov design in Lesson 3.
1. Learning Objectives and Scope
After completing this lesson, students should be able to:
- write an \(n\)-th order SISO plant and reference model in controllable canonical form;
- construct the controller \(u=\hat{\mathbf{k}}_x^T\mathbf{x}+\hat{k}_r r\);
- stack the controller coefficients into a linearly parameterized form \(u=\hat{\boldsymbol{\theta}}^T\boldsymbol{\omega}\);
- derive the matching equations and ideal controller parameters;
- distinguish existence of ideal parameters from knowledge of their values;
- derive the state tracking-error model in terms of parameter error.
The adjustable-parameter update law is intentionally not designed here. Lesson 2 establishes the algebraic structure; Lesson 3 will select the update law by a Lyapunov argument.
2. Canonical-Form Plant and Reference Model
Consider a controllable SISO plant whose full state \(\mathbf{x}\in\mathbb{R}^n\) is measured:
\[ \dot{\mathbf{x}}=\mathbf{A}_p\mathbf{x}+\mathbf{b}_p u, \qquad y=\mathbf{c}^T\mathbf{x}. \]
In controllable canonical form, using the coefficient ordering \(a_0,a_1,\ldots,a_{n-1}\),
\[ \mathbf{A}_p= \begin{bmatrix} 0 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & 1\\ -a_0 & -a_1 & -a_2 & \cdots & -a_{n-1} \end{bmatrix}, \qquad \mathbf{b}_p=b\,\mathbf{e}_n, \]
\[ \mathbf{e}_n= \begin{bmatrix}0&\cdots&0&1\end{bmatrix}^T, \qquad b\neq 0. \]
The desired dynamics are specified by the stable reference model
\[ \dot{\mathbf{x}}_m=\mathbf{A}_m\mathbf{x}_m+\mathbf{b}_m r, \qquad \mathbf{b}_m=b_m\mathbf{e}_n, \]
\[ \mathbf{A}_m= \begin{bmatrix} 0 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & 1\\ -a_{m0} & -a_{m1} & -a_{m2} & \cdots & -a_{m,n-1} \end{bmatrix}. \]
The polynomial \(s^n+a_{m,n-1}s^{n-1}+\cdots+a_{m1}s+a_{m0}\) must be Hurwitz:
\[ \operatorname{Re}\!\left(\lambda_i(\mathbf{A}_m)\right)<0, \qquad i=1,\ldots,n. \]
A common command-scaling choice is \(b_m=a_{m0}\), which gives unit steady-state gain from a constant command to the first canonical state.
3. Adjustable Control-Law Structure
The natural direct-MRAC controller for a full-state canonical plant is
\[ u(t)=\hat{\mathbf{k}}_x^T(t)\mathbf{x}(t)+\hat{k}_r(t)r(t), \]
where \(\hat{\mathbf{k}}_x\in\mathbb{R}^n\) is the adjustable state-feedback vector and \(\hat{k}_r\in\mathbb{R}\) is the adjustable command-feedforward gain. The hats indicate controller variables generated online; they are not assumed to equal plant coefficients.
The controller has two distinct jobs:
- state feedback changes the plant characteristic polynomial so that it matches the reference-model polynomial;
- command feedforward matches the reference input channel.
flowchart TD
R["Reference command r"] --> W["Regressor omega = [x; r]"]
X["Measured plant state x"] --> W
TH["Adjustable parameters \ntheta_hat = [kx_hat; kr_hat]"] --> U["Control u = theta_hat^T omega"]
W --> U
U --> P["Canonical-form SISO plant"]
P --> X
R --> M["Reference model"]
M --> XM["Reference state xm"]
X --> E["Tracking error e = x - xm"]
XM --> E
4. Exact Model-Matching Equations
First suppose that constant ideal controller parameters \(\mathbf{k}_x^*\) and \(k_r^*\) exist. Substitution of \(u=\mathbf{k}_x^{*T}\mathbf{x}+k_r^*r\) into the plant gives
\[ \dot{\mathbf{x}} = \left(\mathbf{A}_p+\mathbf{b}_p\mathbf{k}_x^{*T}\right)\mathbf{x} +\mathbf{b}_p k_r^*r. \]
Exact replication of the reference-model state equation requires
\[ \boxed{ \mathbf{A}_p+\mathbf{b}_p\mathbf{k}_x^{*T}=\mathbf{A}_m }, \qquad \boxed{ \mathbf{b}_p k_r^*=\mathbf{b}_m }. \]
These are the matching conditions. Since both canonical matrices differ only in their final rows, the first equation reduces to \(n\) scalar equations:
\[ -a_i+b k_{x,i}^*=-a_{m,i}, \qquad i=0,\ldots,n-1. \]
Therefore,
\[ \boxed{ k_{x,i}^*=\frac{a_i-a_{m,i}}{b} }, \qquad \boxed{ k_r^*=\frac{b_m}{b} }. \]
Proposition: Existence and Uniqueness in the Canonical Case
For the plant and model defined above, if \(b\neq0\), there is a unique pair \((\mathbf{k}_x^*,k_r^*)\) satisfying the matching equations.
Proof.
The rank-one correction \(\mathbf{b}_p\mathbf{k}_x^{*T} =b\mathbf{e}_n\mathbf{k}_x^{*T}\) modifies only the final row of \(\mathbf{A}_p\). Equality with \(\mathbf{A}_m\) produces the diagonal scalar system
\[ b\mathbf{k}_x^* = \begin{bmatrix} a_0-a_{m0}\\ a_1-a_{m1}\\ \vdots\\ a_{n-1}-a_{m,n-1} \end{bmatrix}. \]
Division by nonzero \(b\) gives one and only one \(\mathbf{k}_x^*\). Likewise, \(bk_r^*=b_m\) gives one and only one \(k_r^*\). Therefore the ideal controller parameter pair exists and is unique. \(\square\)
In adaptive control, the ideal values are generally unknown because \(a_i\) and \(b\) are uncertain. Their existence is a structural assumption; online adaptation attempts to produce controller parameters that recover the associated behavior.
5. Linear Parameterization of the Controller
Define the adjustable parameter vector and measurable regressor as
\[ \hat{\boldsymbol{\theta}} = \begin{bmatrix} \hat{\mathbf{k}}_x\\ \hat{k}_r \end{bmatrix} \in\mathbb{R}^{n+1}, \qquad \boldsymbol{\omega} = \begin{bmatrix} \mathbf{x}\\ r \end{bmatrix} \in\mathbb{R}^{n+1}. \]
The control law becomes
\[ \boxed{u=\hat{\boldsymbol{\theta}}^T\boldsymbol{\omega}}. \]
This expression is linear in the adjustable parameters. That property is central: differentiation with respect to the controller parameters is immediate,
\[ \frac{\partial u}{\partial\hat{\boldsymbol{\theta}}} =\boldsymbol{\omega}. \]
It is important not to confuse three different statements:
- the plant is linear in the state and input;
- the controller is linear in the measured signals;
- the controller is linear in its adjustable parameters.
The third statement is the one needed for standard direct MRAC. More generally, a nonlinear basis vector \(\boldsymbol{\omega}(\mathbf{x},r)\) could still yield a controller that is linear in parameters.
Ideal and Parameter-Error Vectors
\[ \boldsymbol{\theta}^* = \begin{bmatrix} \mathbf{k}_x^*\\ k_r^* \end{bmatrix}, \qquad \tilde{\boldsymbol{\theta}} = \hat{\boldsymbol{\theta}}-\boldsymbol{\theta}^*. \]
Consequently,
\[ u = \boldsymbol{\theta}^{*T}\boldsymbol{\omega} + \tilde{\boldsymbol{\theta}}^T\boldsymbol{\omega}. \]
The first term is the unknown ideal control action, and the second term is the control perturbation caused by parameter mismatch.
6. Derivation of the Parameterized Error Model
Define the state tracking error
\[ \mathbf{e}=\mathbf{x}-\mathbf{x}_m. \]
Using the adjustable controller and adding and subtracting the ideal parameter contribution,
\[ \begin{aligned} \dot{\mathbf{x}} &= \mathbf{A}_p\mathbf{x} +\mathbf{b}_p \hat{\boldsymbol{\theta}}^T\boldsymbol{\omega}\\ &= \mathbf{A}_p\mathbf{x} +\mathbf{b}_p \boldsymbol{\theta}^{*T}\boldsymbol{\omega} +\mathbf{b}_p \tilde{\boldsymbol{\theta}}^T\boldsymbol{\omega}. \end{aligned} \]
The matching equations imply
\[ \mathbf{A}_p\mathbf{x} +\mathbf{b}_p \boldsymbol{\theta}^{*T}\boldsymbol{\omega} = \mathbf{A}_m\mathbf{x}+\mathbf{b}_m r. \]
Hence
\[ \dot{\mathbf{x}} = \mathbf{A}_m\mathbf{x}+\mathbf{b}_m r +\mathbf{b}_p \tilde{\boldsymbol{\theta}}^T\boldsymbol{\omega}. \]
Subtracting the reference-model equation gives the central error model:
\[ \boxed{ \dot{\mathbf{e}} = \mathbf{A}_m\mathbf{e} + \mathbf{b}_p \tilde{\boldsymbol{\theta}}^T\boldsymbol{\omega} }. \]
If the parameter error is zero, the forced term vanishes and the error evolves according to the stable matrix \(\mathbf{A}_m\). For nonzero parameter error, all mismatch enters through the scalar quantity \(\tilde{\boldsymbol{\theta}}^T\boldsymbol{\omega}\) and the known input direction \(\mathbf{b}_p\).
Equivalent Scalar Differential Equation
For the first canonical-state error \(e_1=x_1-x_{m1}\), the vector equation is equivalent to
\[ e_1^{(n)}+ a_{m,n-1}e_1^{(n-1)}+\cdots+a_{m1}\dot{e}_1+a_{m0}e_1 = b\,\tilde{\boldsymbol{\theta}}^T\boldsymbol{\omega}. \]
This representation shows that the stable reference polynomial filters the parameter-mismatch signal before it appears in the tracking error.
7. Sign Conventions and Equivalent Parameterizations
Some texts define negative feedback explicitly:
\[ u=-\hat{\mathbf{q}}^T\mathbf{x}+\hat{q}_r r. \]
This is equivalent to the convention used here under
\[ \hat{\mathbf{k}}_x=-\hat{\mathbf{q}}, \qquad \hat{k}_r=\hat{q}_r. \]
The matching formulas must always be derived from the chosen plant and controller signs rather than memorized. A reliable check is the matrix residual
\[ \mathbf{R}_A= \mathbf{A}_p+\mathbf{b}_p\mathbf{k}_x^{*T}-\mathbf{A}_m, \qquad \mathbf{r}_b=\mathbf{b}_p k_r^*-\mathbf{b}_m. \]
Correct matching requires
\[ \|\mathbf{R}_A\|=0, \qquad \|\mathbf{r}_b\|=0. \]
8. Worked Third-Order Example
Consider
\[ \mathbf{A}_p= \begin{bmatrix} 0&1&0\\ 0&0&1\\ -1.5&-2&-1 \end{bmatrix}, \qquad \mathbf{b}_p= \begin{bmatrix}0\\0\\1.4\end{bmatrix}. \]
Choose the reference model
\[ \mathbf{A}_m= \begin{bmatrix} 0&1&0\\ 0&0&1\\ -6&-11&-6 \end{bmatrix}, \qquad \mathbf{b}_m= \begin{bmatrix}0\\0\\6\end{bmatrix}. \]
The model polynomial is \(s^3+6s^2+11s+6=(s+1)(s+2)(s+3)\). Thus the desired poles are \(-1,-2,-3\).
The ideal feedback parameters are
\[ \mathbf{k}_x^* = \frac{1}{1.4} \begin{bmatrix} 1.5-6\\ 2-11\\ 1-6 \end{bmatrix} = \begin{bmatrix} -3.2142857\\ -6.4285714\\ -3.5714286 \end{bmatrix}, \]
and
\[ k_r^*=\frac{6}{1.4}=4.2857143. \]
Therefore \(\boldsymbol{\theta}^* =[-3.2142857,-6.4285714,-3.5714286,4.2857143]^T\) and \(\boldsymbol{\omega}=[x_1,x_2,x_3,r]^T\).
Substitution into the final plant row gives
\[ \begin{aligned} -1.5+1.4(-3.2142857)&=-6,\\ -2+1.4(-6.4285714)&=-11,\\ -1+1.4(-3.5714286)&=-6,\\ 1.4(4.2857143)&=6. \end{aligned} \]
Thus, under the ideal fixed controller and equal initial plant/model states, the plant and reference-model trajectories are identical.
9. Structural Assumptions and Limitations
The simple parameterization depends on the following assumptions:
- full state availability: every canonical state used in \(\boldsymbol{\omega}\) is measured or otherwise supplied;
- known model order: plant and model have the same state dimension;
- matched input direction: uncertainty enters through the same final canonical channel;
- nonzero high-frequency gain: \(b\neq0\);
- constant ideal parameters: the unknown plant coefficients are treated as constant over the stability analysis;
- reference-model stability: \(\mathbf{A}_m\) is Hurwitz.
The exact magnitude of \(b\) need not be known by a direct adaptive implementation, but the sign of \(b\) will enter the standard Lyapunov update law. Chapter 7, Lesson 4 treats the sign-definite high-frequency-gain issue explicitly.
If the plant is not in canonical coordinates but is controllable, a similarity transformation may produce the required form. If only output measurements are available, a dynamic output-feedback parameterization is needed; that topic is deferred to Chapter 15.
10. Verification Workflow
flowchart TD
A["Write plant Ap, bp in controllable canonical form"] --> B["Choose Hurwitz model Am and command gain bm"]
B --> C["Assume u = kx_hat^T x + kr_hat r"]
C --> D["Solve ideal matching equations for kx_star and kr_star"]
D --> E["Stack theta_hat and omega so u = theta_hat^T omega"]
E --> F["Check matrix and input-channel residuals"]
F --> G["Compare ideal and deliberately mismatched fixed controllers"]
G --> H["Use parameterized error model in Lesson 3"]
The accompanying programs implement this workflow. They do not use an adaptive law. Instead, they compare the ideal parameter vector with a deliberately perturbed fixed vector, making the role of exact matching visible before adaptation is introduced.
11. Python Implementation
This implementation uses NumPy for vector operations and a from-scratch fourth-order Runge-Kutta integrator. Matplotlib is used only for plotting.
Chapter7_Lesson2.py
"""
Chapter7_Lesson2.py
Structural verification of the MRAC control law parameterization for a
third-order SISO plant in controllable canonical form.
This lesson intentionally uses fixed parameter vectors. The adaptive update
law is developed in Chapter 7, Lesson 3.
"""
from __future__ import annotations
import numpy as np
import matplotlib.pyplot as plt
PLANT_A = np.array([1.5, 2.0, 1.0], dtype=float) # [a0, a1, a2]
MODEL_A = np.array([6.0, 11.0, 6.0], dtype=float) # [am0, am1, am2]
B = 1.4
BM = 6.0
def ideal_parameters(
plant_a: np.ndarray,
model_a: np.ndarray,
b: float,
bm: float,
) -> np.ndarray:
"""Return theta* = [kx*, kr*] from the canonical-form matching equations."""
if abs(b) < 1.0e-12:
raise ValueError("The plant high-frequency gain b must be nonzero.")
kx_star = (plant_a - model_a) / b
kr_star = bm / b
return np.concatenate((kx_star, np.array([kr_star])))
def reference(t: float) -> float:
"""Bounded command used for structural verification."""
return 1.0 + 0.25 * np.sin(0.6 * t)
def derivative(
t: float,
z: np.ndarray,
theta_ideal: np.ndarray,
theta_fixed: np.ndarray,
) -> np.ndarray:
"""Combined dynamics: ideal plant, mismatched plant, and reference model."""
x_ideal = z[0:3]
x_fixed = z[3:6]
x_model = z[6:9]
r = reference(t)
omega_ideal = np.concatenate((x_ideal, np.array([r])))
omega_fixed = np.concatenate((x_fixed, np.array([r])))
u_ideal = float(theta_ideal @ omega_ideal)
u_fixed = float(theta_fixed @ omega_fixed)
dx_ideal = np.array(
[
x_ideal[1],
x_ideal[2],
-PLANT_A @ x_ideal + B * u_ideal,
]
)
dx_fixed = np.array(
[
x_fixed[1],
x_fixed[2],
-PLANT_A @ x_fixed + B * u_fixed,
]
)
dx_model = np.array(
[
x_model[1],
x_model[2],
-MODEL_A @ x_model + BM * r,
]
)
return np.concatenate((dx_ideal, dx_fixed, dx_model))
def rk4_step(
fun,
t: float,
z: np.ndarray,
dt: float,
theta_ideal: np.ndarray,
theta_fixed: np.ndarray,
) -> np.ndarray:
"""One classical fourth-order Runge-Kutta step."""
k1 = fun(t, z, theta_ideal, theta_fixed)
k2 = fun(t + 0.5 * dt, z + 0.5 * dt * k1, theta_ideal, theta_fixed)
k3 = fun(t + 0.5 * dt, z + 0.5 * dt * k2, theta_ideal, theta_fixed)
k4 = fun(t + dt, z + dt * k3, theta_ideal, theta_fixed)
return z + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
def matching_residual(theta: np.ndarray) -> tuple[np.ndarray, float]:
"""Return coefficient and command-channel residuals."""
kx = theta[:3]
kr = theta[3]
coefficient_residual = -PLANT_A + B * kx + MODEL_A
input_residual = B * kr - BM
return coefficient_residual, float(input_residual)
def main() -> None:
theta_star = ideal_parameters(PLANT_A, MODEL_A, B, BM)
theta_fixed = theta_star + np.array([0.8, -0.5, 0.4, -0.8])
coeff_residual, input_residual = matching_residual(theta_star)
print("theta* =", theta_star)
print("||coefficient matching residual||_2 =", np.linalg.norm(coeff_residual))
print("|input matching residual| =", abs(input_residual))
dt = 0.002
final_time = 20.0
time = np.arange(0.0, final_time + dt, dt)
trajectory = np.zeros((time.size, 9), dtype=float)
for index in range(time.size - 1):
trajectory[index + 1] = rk4_step(
derivative,
time[index],
trajectory[index],
dt,
theta_star,
theta_fixed,
)
x_ideal = trajectory[:, 0]
x_fixed = trajectory[:, 3]
x_model = trajectory[:, 6]
ideal_rms = np.sqrt(np.mean((x_ideal - x_model) ** 2))
fixed_rms = np.sqrt(np.mean((x_fixed - x_model) ** 2))
print("RMS output error with theta*:", ideal_rms)
print("RMS output error with fixed mismatch:", fixed_rms)
np.savetxt(
"Chapter7_Lesson2_results.csv",
np.column_stack((time, x_ideal, x_fixed, x_model)),
delimiter=",",
header="time,x_ideal,x_fixed,x_model",
comments="",
)
plt.figure()
plt.plot(time, x_model, label="reference model")
plt.plot(time, x_ideal, "--", label="plant with theta*")
plt.plot(time, x_fixed, ":", label="plant with fixed mismatch")
plt.xlabel("Time (s)")
plt.ylabel("First canonical state")
plt.title("MRAC structural matching verification")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()
if __name__ == "__main__":
main()
12. C++ Implementation
The C++ version uses only the standard library, writes a CSV file, and implements the same Runge-Kutta integration from scratch.
Chapter7_Lesson2.cpp
/*
Chapter7_Lesson2.cpp
Structural verification of the MRAC control law parameterization for a
third-order SISO plant in controllable canonical form.
This lesson intentionally uses fixed parameter vectors. The adaptive update
law is developed in Chapter 7, Lesson 3.
*/
#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <stdexcept>
using Vec4 = std::array<double, 4>;
using Vec9 = std::array<double, 9>;
constexpr std::array<double, 3> PLANT_A{1.5, 2.0, 1.0};
constexpr std::array<double, 3> MODEL_A{6.0, 11.0, 6.0};
constexpr double B = 1.4;
constexpr double BM = 6.0;
double reference(double t) {
return 1.0 + 0.25 * std::sin(0.6 * t);
}
Vec4 idealParameters() {
if (std::abs(B) < 1.0e-12) {
throw std::runtime_error("The plant high-frequency gain b must be nonzero.");
}
return Vec4{
(PLANT_A[0] - MODEL_A[0]) / B,
(PLANT_A[1] - MODEL_A[1]) / B,
(PLANT_A[2] - MODEL_A[2]) / B,
BM / B
};
}
double control(const Vec4& theta, const Vec9& z, int offset, double r) {
return theta[0] * z[offset]
+ theta[1] * z[offset + 1]
+ theta[2] * z[offset + 2]
+ theta[3] * r;
}
Vec9 derivative(double t, const Vec9& z, const Vec4& thetaIdeal, const Vec4& thetaFixed) {
const double r = reference(t);
const double uIdeal = control(thetaIdeal, z, 0, r);
const double uFixed = control(thetaFixed, z, 3, r);
Vec9 dz{};
dz[0] = z[1];
dz[1] = z[2];
dz[2] = -PLANT_A[0] * z[0] - PLANT_A[1] * z[1] - PLANT_A[2] * z[2] + B * uIdeal;
dz[3] = z[4];
dz[4] = z[5];
dz[5] = -PLANT_A[0] * z[3] - PLANT_A[1] * z[4] - PLANT_A[2] * z[5] + B * uFixed;
dz[6] = z[7];
dz[7] = z[8];
dz[8] = -MODEL_A[0] * z[6] - MODEL_A[1] * z[7] - MODEL_A[2] * z[8] + BM * r;
return dz;
}
Vec9 addScaled(const Vec9& x, const Vec9& y, double scale) {
Vec9 result{};
for (std::size_t i = 0; i < result.size(); ++i) {
result[i] = x[i] + scale * y[i];
}
return result;
}
Vec9 rk4Step(double t, const Vec9& z, double dt, const Vec4& thetaIdeal, const Vec4& thetaFixed) {
const Vec9 k1 = derivative(t, z, thetaIdeal, thetaFixed);
const Vec9 k2 = derivative(t + 0.5 * dt, addScaled(z, k1, 0.5 * dt), thetaIdeal, thetaFixed);
const Vec9 k3 = derivative(t + 0.5 * dt, addScaled(z, k2, 0.5 * dt), thetaIdeal, thetaFixed);
const Vec9 k4 = derivative(t + dt, addScaled(z, k3, dt), thetaIdeal, thetaFixed);
Vec9 next{};
for (std::size_t i = 0; i < next.size(); ++i) {
next[i] = z[i] + (dt / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
}
return next;
}
int main() {
try {
const Vec4 thetaStar = idealParameters();
Vec4 thetaFixed = thetaStar;
thetaFixed[0] += 0.8;
thetaFixed[1] -= 0.5;
thetaFixed[2] += 0.4;
thetaFixed[3] -= 0.8;
double coefficientResidualSquared = 0.0;
for (int i = 0; i < 3; ++i) {
const double residual = -PLANT_A[i] + B * thetaStar[i] + MODEL_A[i];
coefficientResidualSquared += residual * residual;
}
const double inputResidual = B * thetaStar[3] - BM;
std::cout << std::setprecision(12);
std::cout << "theta* = [";
for (std::size_t i = 0; i < thetaStar.size(); ++i) {
std::cout << thetaStar[i] << (i + 1 == thetaStar.size() ? "]\n" : ", ");
}
std::cout << "||coefficient matching residual||_2 = "
<< std::sqrt(coefficientResidualSquared) << '\n';
std::cout << "|input matching residual| = " << std::abs(inputResidual) << '\n';
const double dt = 0.002;
const double finalTime = 20.0;
const int steps = static_cast<int>(std::round(finalTime / dt));
Vec9 z{};
std::ofstream csv("Chapter7_Lesson2_results.csv");
if (!csv) {
throw std::runtime_error("Could not create the CSV output file.");
}
csv << "time,x_ideal,x_fixed,x_model\n";
double idealErrorEnergy = 0.0;
double fixedErrorEnergy = 0.0;
for (int step = 0; step <= steps; ++step) {
const double t = step * dt;
csv << t << ',' << z[0] << ',' << z[3] << ',' << z[6] << '\n';
const double eIdeal = z[0] - z[6];
const double eFixed = z[3] - z[6];
idealErrorEnergy += eIdeal * eIdeal;
fixedErrorEnergy += eFixed * eFixed;
if (step < steps) {
z = rk4Step(t, z, dt, thetaStar, thetaFixed);
}
}
const double sampleCount = static_cast<double>(steps + 1);
std::cout << "RMS output error with theta* = "
<< std::sqrt(idealErrorEnergy / sampleCount) << '\n';
std::cout << "RMS output error with fixed mismatch = "
<< std::sqrt(fixedErrorEnergy / sampleCount) << '\n';
return 0;
} catch (const std::exception& error) {
std::cerr << "Error: " << error.what() << '\n';
return 1;
}
}
13. Java Implementation
The Java version mirrors the C++ structure with primitive arrays and produces the same CSV columns.
Chapter7_Lesson2.java
/*
Chapter7_Lesson2.java
Structural verification of the MRAC control law parameterization for a
third-order SISO plant in controllable canonical form.
This lesson intentionally uses fixed parameter vectors. The adaptive update
law is developed in Chapter 7, Lesson 3.
*/
import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.util.Arrays;
public final class Chapter7_Lesson2 {
private static final double[] PLANT_A = {1.5, 2.0, 1.0};
private static final double[] MODEL_A = {6.0, 11.0, 6.0};
private static final double B = 1.4;
private static final double BM = 6.0;
private Chapter7_Lesson2() {
}
private static double reference(double t) {
return 1.0 + 0.25 * Math.sin(0.6 * t);
}
private static double[] idealParameters() {
if (Math.abs(B) < 1.0e-12) {
throw new IllegalArgumentException("The plant high-frequency gain b must be nonzero.");
}
return new double[] {
(PLANT_A[0] - MODEL_A[0]) / B,
(PLANT_A[1] - MODEL_A[1]) / B,
(PLANT_A[2] - MODEL_A[2]) / B,
BM / B
};
}
private static double control(double[] theta, double[] z, int offset, double r) {
return theta[0] * z[offset]
+ theta[1] * z[offset + 1]
+ theta[2] * z[offset + 2]
+ theta[3] * r;
}
private static double[] derivative(
double t,
double[] z,
double[] thetaIdeal,
double[] thetaFixed
) {
double r = reference(t);
double uIdeal = control(thetaIdeal, z, 0, r);
double uFixed = control(thetaFixed, z, 3, r);
double[] dz = new double[9];
dz[0] = z[1];
dz[1] = z[2];
dz[2] = -PLANT_A[0] * z[0] - PLANT_A[1] * z[1] - PLANT_A[2] * z[2] + B * uIdeal;
dz[3] = z[4];
dz[4] = z[5];
dz[5] = -PLANT_A[0] * z[3] - PLANT_A[1] * z[4] - PLANT_A[2] * z[5] + B * uFixed;
dz[6] = z[7];
dz[7] = z[8];
dz[8] = -MODEL_A[0] * z[6] - MODEL_A[1] * z[7] - MODEL_A[2] * z[8] + BM * r;
return dz;
}
private static double[] addScaled(double[] x, double[] y, double scale) {
double[] result = new double[x.length];
for (int i = 0; i < x.length; i++) {
result[i] = x[i] + scale * y[i];
}
return result;
}
private static double[] rk4Step(
double t,
double[] z,
double dt,
double[] thetaIdeal,
double[] thetaFixed
) {
double[] k1 = derivative(t, z, thetaIdeal, thetaFixed);
double[] k2 = derivative(
t + 0.5 * dt,
addScaled(z, k1, 0.5 * dt),
thetaIdeal,
thetaFixed
);
double[] k3 = derivative(
t + 0.5 * dt,
addScaled(z, k2, 0.5 * dt),
thetaIdeal,
thetaFixed
);
double[] k4 = derivative(
t + dt,
addScaled(z, k3, dt),
thetaIdeal,
thetaFixed
);
double[] next = new double[z.length];
for (int i = 0; i < z.length; i++) {
next[i] = z[i] + (dt / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
}
return next;
}
public static void main(String[] args) {
double[] thetaStar = idealParameters();
double[] thetaFixed = thetaStar.clone();
thetaFixed[0] += 0.8;
thetaFixed[1] -= 0.5;
thetaFixed[2] += 0.4;
thetaFixed[3] -= 0.8;
double residualSquared = 0.0;
for (int i = 0; i < 3; i++) {
double residual = -PLANT_A[i] + B * thetaStar[i] + MODEL_A[i];
residualSquared += residual * residual;
}
double inputResidual = B * thetaStar[3] - BM;
System.out.println("theta* = " + Arrays.toString(thetaStar));
System.out.println("||coefficient matching residual||_2 = " + Math.sqrt(residualSquared));
System.out.println("|input matching residual| = " + Math.abs(inputResidual));
double dt = 0.002;
double finalTime = 20.0;
int steps = (int) Math.round(finalTime / dt);
double[] z = new double[9];
double idealErrorEnergy = 0.0;
double fixedErrorEnergy = 0.0;
try (BufferedWriter writer = new BufferedWriter(
new FileWriter("Chapter7_Lesson2_results.csv")
)) {
writer.write("time,x_ideal,x_fixed,x_model");
writer.newLine();
for (int step = 0; step <= steps; step++) {
double t = step * dt;
writer.write(t + "," + z[0] + "," + z[3] + "," + z[6]);
writer.newLine();
double eIdeal = z[0] - z[6];
double eFixed = z[3] - z[6];
idealErrorEnergy += eIdeal * eIdeal;
fixedErrorEnergy += eFixed * eFixed;
if (step < steps) {
z = rk4Step(t, z, dt, thetaStar, thetaFixed);
}
}
} catch (IOException error) {
System.err.println("Could not write the CSV output: " + error.getMessage());
System.exit(1);
}
double sampleCount = steps + 1.0;
System.out.println(
"RMS output error with theta* = "
+ Math.sqrt(idealErrorEnergy / sampleCount)
);
System.out.println(
"RMS output error with fixed mismatch = "
+ Math.sqrt(fixedErrorEnergy / sampleCount)
);
}
}
14. MATLAB Implementation
The MATLAB script uses ode45, verifies both matching
residuals, plots the trajectories, and exports the numerical results.
Chapter7_Lesson2.m
% Chapter7_Lesson2.m
% Structural verification of the MRAC control law parameterization for a
% third-order SISO plant in controllable canonical form.
%
% This lesson intentionally uses fixed parameter vectors. The adaptive update
% law is developed in Chapter 7, Lesson 3.
clear;
clc;
close all;
plantA = [1.5; 2.0; 1.0]; % [a0; a1; a2]
modelA = [6.0; 11.0; 6.0]; % [am0; am1; am2]
b = 1.4;
bm = 6.0;
if abs(b) < 1.0e-12
error('The plant high-frequency gain b must be nonzero.');
end
thetaStar = [(plantA - modelA) / b; bm / b];
thetaFixed = thetaStar + [0.8; -0.5; 0.4; -0.8];
coefficientResidual = -plantA + b * thetaStar(1:3) + modelA;
inputResidual = b * thetaStar(4) - bm;
fprintf('theta* = [%g, %g, %g, %g]^T\n', thetaStar);
fprintf('||coefficient matching residual||_2 = %.12g\n', norm(coefficientResidual));
fprintf('|input matching residual| = %.12g\n', abs(inputResidual));
reference = @(t) 1.0 + 0.25 * sin(0.6 * t);
dynamics = @(t, z) combinedDynamics( ...
t, z, thetaStar, thetaFixed, plantA, modelA, b, bm, reference);
timeGrid = linspace(0.0, 20.0, 10001);
initialState = zeros(9, 1);
options = odeset('RelTol', 1.0e-9, 'AbsTol', 1.0e-11);
[t, z] = ode45(dynamics, timeGrid, initialState, options);
xIdeal = z(:, 1);
xFixed = z(:, 4);
xModel = z(:, 7);
idealRms = sqrt(mean((xIdeal - xModel).^2));
fixedRms = sqrt(mean((xFixed - xModel).^2));
fprintf('RMS output error with theta* = %.12g\n', idealRms);
fprintf('RMS output error with fixed mismatch = %.12g\n', fixedRms);
results = table(t, xIdeal, xFixed, xModel);
writetable(results, 'Chapter7_Lesson2_results.csv');
figure;
plot(t, xModel, 'LineWidth', 1.6);
hold on;
plot(t, xIdeal, '--', 'LineWidth', 1.4);
plot(t, xFixed, ':', 'LineWidth', 1.4);
grid on;
xlabel('Time (s)');
ylabel('First canonical state');
title('MRAC structural matching verification');
legend('Reference model', 'Plant with theta*', ...
'Plant with fixed mismatch', 'Location', 'best');
function dz = combinedDynamics( ...
t, z, thetaIdeal, thetaFixed, plantA, modelA, b, bm, reference)
xIdeal = z(1:3);
xFixed = z(4:6);
xModel = z(7:9);
r = reference(t);
omegaIdeal = [xIdeal; r];
omegaFixed = [xFixed; r];
uIdeal = thetaIdeal.' * omegaIdeal;
uFixed = thetaFixed.' * omegaFixed;
dxIdeal = [ ...
xIdeal(2);
xIdeal(3);
-plantA.' * xIdeal + b * uIdeal
];
dxFixed = [ ...
xFixed(2);
xFixed(3);
-plantA.' * xFixed + b * uFixed
];
dxModel = [ ...
xModel(2);
xModel(3);
-modelA.' * xModel + bm * r
];
dz = [dxIdeal; dxFixed; dxModel];
end
15. Simulink Model Builder
The following MATLAB script programmatically constructs a Simulink model from standard blocks. It exposes the regressor construction and the matrix gain implementing \(u=\boldsymbol{\theta}^{*T}\boldsymbol{\omega}\).
Chapter7_Lesson2_Simulink.m
% Chapter7_Lesson2_Simulink.m
% Programmatically builds a Simulink model that exposes the parameterized
% control law u = theta_hat' * [x; r].
%
% Requirements: MATLAB, Simulink. The script uses standard Simulink blocks
% and does not require the Model Reference Adaptive Control block.
clear;
clc;
modelName = 'Chapter7_Lesson2_Simulink_Model';
if bdIsLoaded(modelName)
close_system(modelName, 0);
end
new_system(modelName);
open_system(modelName);
plantA = [1.5; 2.0; 1.0];
modelA = [6.0; 11.0; 6.0];
b = 1.4;
bm = 6.0;
thetaStar = [(plantA - modelA) / b; bm / b];
Ap = [0 1 0; 0 0 1; -plantA.'];
Bp = [0; 0; b];
Cp = [1 0 0];
Dp = 0;
Am = [0 1 0; 0 0 1; -modelA.'];
Bm = [0; 0; bm];
Cm = eye(3);
Dm = zeros(3, 1);
assignin('base', 'Ap', Ap);
assignin('base', 'Bp', Bp);
assignin('base', 'Cp', Cp);
assignin('base', 'Dp', Dp);
assignin('base', 'Am', Am);
assignin('base', 'Bm', Bm);
assignin('base', 'Cm', Cm);
assignin('base', 'Dm', Dm);
assignin('base', 'thetaStar', thetaStar);
add_block('simulink/Sources/Sine Wave', [modelName '/Reference'], ...
'Amplitude', '0.25', 'Bias', '1', 'Frequency', '0.6', ...
'Position', [40 120 90 150]);
add_block('simulink/Continuous/State-Space', [modelName '/Plant'], ...
'A', 'Ap', 'B', 'Bp', 'C', 'eye(3)', 'D', 'zeros(3,1)', ...
'Position', [420 60 540 130]);
add_block('simulink/Continuous/State-Space', [modelName '/Reference Model'], ...
'A', 'Am', 'B', 'Bm', 'C', 'Cm', 'D', 'Dm', ...
'Position', [420 210 540 280]);
add_block('simulink/Signal Routing/Mux', [modelName '/Regressor Mux'], ...
'Inputs', '2', 'Position', [190 50 195 150]);
add_block('simulink/Math Operations/Gain', [modelName '/thetaStar Transpose'], ...
'Gain', 'transpose(thetaStar)', 'Multiplication', 'Matrix(K*u)', ...
'Position', [260 75 360 125]);
add_block('simulink/Signal Routing/Selector', [modelName '/Plant Output'], ...
'Indices', '1', 'InputPortWidth', '3', ...
'Position', [590 70 635 115]);
add_block('simulink/Signal Routing/Selector', [modelName '/Model Output'], ...
'Indices', '1', 'InputPortWidth', '3', ...
'Position', [590 220 635 265]);
add_block('simulink/Math Operations/Sum', [modelName '/Tracking Error'], ...
'Inputs', '+-', 'Position', [690 130 720 180]);
add_block('simulink/Sinks/Scope', [modelName '/Scope'], ...
'NumInputPorts', '3', 'Position', [790 90 840 190]);
add_line(modelName, 'Reference/1', 'Regressor Mux/2', 'autorouting', 'on');
add_line(modelName, 'Reference/1', 'Reference Model/1', 'autorouting', 'on');
add_line(modelName, 'Plant/1', 'Regressor Mux/1', 'autorouting', 'on');
add_line(modelName, 'Regressor Mux/1', 'thetaStar Transpose/1', 'autorouting', 'on');
add_line(modelName, 'thetaStar Transpose/1', 'Plant/1', 'autorouting', 'on');
add_line(modelName, 'Plant/1', 'Plant Output/1', 'autorouting', 'on');
add_line(modelName, 'Reference Model/1', 'Model Output/1', 'autorouting', 'on');
add_line(modelName, 'Plant Output/1', 'Tracking Error/1', 'autorouting', 'on');
add_line(modelName, 'Model Output/1', 'Tracking Error/2', 'autorouting', 'on');
add_line(modelName, 'Plant Output/1', 'Scope/1', 'autorouting', 'on');
add_line(modelName, 'Model Output/1', 'Scope/2', 'autorouting', 'on');
add_line(modelName, 'Tracking Error/1', 'Scope/3', 'autorouting', 'on');
set_param(modelName, 'StopTime', '20', 'Solver', 'ode45');
save_system(modelName, [modelName '.slx']);
open_system(modelName);
16. Wolfram Mathematica Implementation
The notebook expression below can be saved directly as an
.nb file and opened in Wolfram Mathematica.
Chapter7_Lesson2.nb
Notebook[{
Cell["Chapter 7, Lesson 2: Control Law Structure and Parameterization", "Title"],
Cell[
"Structural verification of ideal matching for a third-order SISO plant in controllable canonical form. The parameter vector is fixed; the adaptive update law belongs to Lesson 3.",
"Text"
],
Cell[
BoxData @ ToBoxes @ Unevaluated[
ClearAll["Global`*"];
plantA = {1.5, 2.0, 1.0};
modelA = {6.0, 11.0, 6.0};
b = 1.4;
bm = 6.0;
If[Abs[b] < 10^-12,
Print["The plant high-frequency gain b must be nonzero."];
Abort[]
];
thetaStar = Join[(plantA - modelA)/b, {bm/b}];
thetaFixed = thetaStar + {0.8, -0.5, 0.4, -0.8};
coefficientResidual =
-plantA + b thetaStar[[1 ;; 3]] + modelA;
inputResidual = b thetaStar[[4]] - bm;
Print["theta* = ", thetaStar];
Print[
"||coefficient matching residual||_2 = ",
Norm[coefficientResidual]
];
Print["|input matching residual| = ", Abs[inputResidual]];
reference[t_] := 1.0 + 0.25 Sin[0.6 t];
solution = NDSolveValue[
{
xi1'[t] == xi2[t],
xi2'[t] == xi3[t],
xi3'[t] ==
-plantA.{xi1[t], xi2[t], xi3[t]} +
b (thetaStar.{xi1[t], xi2[t], xi3[t], reference[t]}),
xf1'[t] == xf2[t],
xf2'[t] == xf3[t],
xf3'[t] ==
-plantA.{xf1[t], xf2[t], xf3[t]} +
b (thetaFixed.{xf1[t], xf2[t], xf3[t], reference[t]}),
xm1'[t] == xm2[t],
xm2'[t] == xm3[t],
xm3'[t] ==
-modelA.{xm1[t], xm2[t], xm3[t]} + bm reference[t],
xi1[0] == 0, xi2[0] == 0, xi3[0] == 0,
xf1[0] == 0, xf2[0] == 0, xf3[0] == 0,
xm1[0] == 0, xm2[0] == 0, xm3[0] == 0
},
{xi1, xi2, xi3, xf1, xf2, xf3, xm1, xm2, xm3},
{t, 0, 20},
Method -> {"TimeIntegration" -> {"ExplicitRungeKutta"}}
];
idealErrorRMS = Sqrt[
NIntegrate[(solution[[1]][t] - solution[[7]][t])^2, {t, 0, 20}]/20
];
fixedErrorRMS = Sqrt[
NIntegrate[(solution[[4]][t] - solution[[7]][t])^2, {t, 0, 20}]/20
];
Print["RMS output error with theta* = ", idealErrorRMS];
Print["RMS output error with fixed mismatch = ", fixedErrorRMS];
Plot[
Evaluate[{solution[[7]][t], solution[[1]][t], solution[[4]][t]}],
{t, 0, 20},
PlotLegends -> {
"Reference model",
"Plant with theta*",
"Plant with fixed mismatch"
},
AxesLabel -> {"Time (s)", "First canonical state"},
PlotLabel -> "MRAC structural matching verification",
GridLines -> Automatic,
ImageSize -> Large
]
],
"Input"
]
}]
17. Problems and Solutions
Problem 1: Second-Order Matching
A plant has \(\ddot{x}_1=-2x_1-3\dot{x}_1+2u\). The reference model is \(\ddot{x}_{m1}=-6x_{m1}-5\dot{x}_{m1}+6r\). For \(u=k_1x_1+k_2\dot{x}_1+k_rr\), calculate the ideal controller parameters.
Solution.
Here \(\mathbf{a}=[2,3]^T\), \(\mathbf{a}_m=[6,5]^T\), \(b=2\), and \(b_m=6\). Therefore
\[ \mathbf{k}_x^* = \frac{\mathbf{a}-\mathbf{a}_m}{b} = \frac{1}{2} \begin{bmatrix}-4\\-2\end{bmatrix} = \begin{bmatrix}-2\\-1\end{bmatrix}, \qquad k_r^*=\frac{6}{2}=3. \]
Substitution yields \(\ddot{x}_1=-6x_1-5\dot{x}_1+6r\), exactly the reference-model dynamics.
Problem 2: Regressor and Parameter Dimensions
For a fifth-order canonical plant, write \(\hat{\boldsymbol{\theta}}\) and \(\boldsymbol{\omega}\), and state their dimensions.
Solution.
\[ \hat{\boldsymbol{\theta}} = \begin{bmatrix} \hat{k}_1&\hat{k}_2&\hat{k}_3&\hat{k}_4& \hat{k}_5&\hat{k}_r \end{bmatrix}^T, \]
\[ \boldsymbol{\omega} = \begin{bmatrix} x_1&x_2&x_3&x_4&x_5&r \end{bmatrix}^T. \]
Both vectors belong to \(\mathbb{R}^6\), so their inner product is a scalar control input.
Problem 3: Error-Model Derivation
Starting from \(\dot{\mathbf{x}}=\mathbf{A}_p\mathbf{x}+ \mathbf{b}_p\hat{\boldsymbol{\theta}}^T\boldsymbol{\omega}\), prove that the matching conditions imply \(\dot{\mathbf{e}}=\mathbf{A}_m\mathbf{e}+ \mathbf{b}_p\tilde{\boldsymbol{\theta}}^T\boldsymbol{\omega}\).
Solution.
Write
\[ \hat{\boldsymbol{\theta}} = \boldsymbol{\theta}^*+\tilde{\boldsymbol{\theta}}. \]
Then
\[ \dot{\mathbf{x}} = \mathbf{A}_p\mathbf{x}+ \mathbf{b}_p\boldsymbol{\theta}^{*T}\boldsymbol{\omega}+ \mathbf{b}_p\tilde{\boldsymbol{\theta}}^T\boldsymbol{\omega}. \]
The matching conditions replace the first two terms by \(\mathbf{A}_m\mathbf{x}+\mathbf{b}_m r\). Subtracting \(\dot{\mathbf{x}}_m=\mathbf{A}_m\mathbf{x}_m+ \mathbf{b}_m r\) gives
\[ \dot{\mathbf{e}} = \mathbf{A}_m(\mathbf{x}-\mathbf{x}_m)+ \mathbf{b}_p\tilde{\boldsymbol{\theta}}^T\boldsymbol{\omega} = \mathbf{A}_m\mathbf{e}+ \mathbf{b}_p\tilde{\boldsymbol{\theta}}^T\boldsymbol{\omega}. \]
Problem 4: Third-Order Numerical Residual
For the worked example, use the perturbed fixed vector \(\hat{\boldsymbol{\theta}}= \boldsymbol{\theta}^*+[0.8,-0.5,0.4,-0.8]^T\). Determine the final-row coefficient residual and input-gain residual.
Solution.
Since the ideal residual is zero, only the perturbation contributes:
\[ \mathbf{r}_a = b \begin{bmatrix}0.8\\-0.5\\0.4\end{bmatrix} = \begin{bmatrix}1.12\\-0.70\\0.56\end{bmatrix}, \qquad r_b=b(-0.8)=-1.12. \]
Thus the fixed perturbed controller no longer reproduces either the model polynomial or its command channel.
Problem 5: Negative-Feedback Convention
Suppose a text uses \(u=-\hat{\mathbf{q}}^T\mathbf{x}+\hat{q}_rr\). Derive the ideal \(\mathbf{q}^*\) for the canonical plant used in this lesson.
Solution.
The state-feedback coefficient in the present convention is \(\mathbf{k}_x^*=(\mathbf{a}-\mathbf{a}_m)/b\). Since \(\mathbf{k}_x^*=-\mathbf{q}^*\),
\[ \boxed{ \mathbf{q}^*=\frac{\mathbf{a}_m-\mathbf{a}}{b} }, \qquad \boxed{ q_r^*=\frac{b_m}{b} }. \]
This confirms why sign conventions must be handled algebraically rather than by copying a formula without checking the controller definition.
18. Summary
For a full-state SISO plant in controllable canonical form, the direct MRAC controller naturally separates into adjustable state-feedback and command-feedforward gains. Canonical structure reduces exact model-matching to scalar coefficient equations, yielding the conceptual ideal parameters \(\mathbf{k}_x^*=(\mathbf{a}-\mathbf{a}_m)/b\) and \(k_r^*=b_m/b\). Stacking the controller coefficients and measurable signals gives \(u=\hat{\boldsymbol{\theta}}^T\boldsymbol{\omega}\), a linear parameterization. Under the matching conditions, the tracking error satisfies \(\dot{\mathbf{e}}=\mathbf{A}_m\mathbf{e}+ \mathbf{b}_p\tilde{\boldsymbol{\theta}}^T\boldsymbol{\omega}\). Lesson 3 will use this equation to construct and prove a Lyapunov-based adaptive update law.
19. References
- 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. (1976). Stable adaptive observers and controllers. Proceedings of the IEEE, 64(8), 1198–1208.
- Feuer, A., & Morse, A. S. (1978). Adaptive control of single-input, single-output linear systems. IEEE Transactions on Automatic Control, 23(4), 557–569.
- Narendra, K. S., & Valavani, L. S. (1978). Stable adaptive controller design—Direct control. IEEE Transactions on Automatic Control, 23(4), 570–582.
- Morse, A. S. (1980). Global stability of parameter-adaptive control systems. IEEE Transactions on Automatic Control, 25(3), 433–439.
- 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.
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.
