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

• define a continuous-time reference model in state-space and transfer-function form;
• translate time-domain specifications into reference-model poles;
• distinguish a command signal from the desired model output;
• verify stability, properness, relative degree, DC gain, and bandwidth consistency;
• derive bounded-input bounded-state and input-to-state bounds for the model;
• explain why an arbitrarily fast reference model may be mathematically stable but physically infeasible;
• simulate and compare candidate reference models before introducing an adaptive law.

Students are assumed to know pole placement, controllability, transfer functions, state-space realization, second-order response metrics, and Lyapunov stability. The earlier chapters have also introduced parametric uncertainty, linearly parameterized models, ideal parameters, and matching conditions. No parameter-update law is designed in this lesson.

In MRAC, the command \( r(t) \) is not normally compared directly with the plant output. Instead, it drives a deliberately chosen stable dynamical system, called the reference model. Its state and output are denoted by \( x_m(t) \) and \( y_m(t) \). The plant output \( y(t) \) will later be required to approach \( y_m(t) \).

\[ \dot{x}_m(t)=A_m x_m(t)+B_m r(t), \qquad y_m(t)=C_m x_m(t)+D_m r(t). \]

The desired tracking error is therefore defined conceptually as \( e(t)=y(t)-y_m(t) \). The reference model is a performance contract: it encodes desired speed, damping, steady-state gain, relative degree, and command filtering. Adaptation is the later mechanism used to make the uncertain plant satisfy that contract.

flowchart TD
  R["Command r(t)"] --> RM["Stable reference model"]
  RM --> YM["Desired output ym(t)"]
  R --> AC["Adaptive controller and uncertain plant"]
  AC --> Y["Actual output y(t)"]
  Y --> E["Tracking comparison: e(t) = y(t) - ym(t)"]
  YM --> E
  E --> U["Adaptation mechanism introduced in later lessons"]
        

A reference model is not an estimate of the plant. It is also not an observer. It is a user-selected dynamical specification. The plant may have unknown parameters, while the reference model parameters must be known and fixed at this stage.

3.1 Internal stability

The matrix \( A_m \) must be Hurwitz in continuous time:

\[ \operatorname{Re}\{\lambda_i(A_m)\} < 0, \qquad i=1,\ldots,n_m. \]

Consequently, there exist constants \( k_m\ge 1 \) and \( \alpha_m>0 \) such that

\[ \left\|e^{A_m t}\right\| \le k_m e^{-\alpha_m t}, \qquad t\ge 0. \]

3.2 Properness and causality

For a SISO transfer model \( M(s)=Y_m(s)/R(s) \), causality requires \( M(s) \) to be proper. A strictly proper model has \( D_m=0 \) and filters command discontinuities rather than passing them instantaneously to the desired output.

\[ M(s)=C_m(sI-A_m)^{-1}B_m+D_m. \]

3.3 Desired steady-state gain

For constant commands and nonsingular \( A_m \), the equilibrium is

\[ x_{m,\mathrm{ss}}=-A_m^{-1}B_m r_0, \qquad y_{m,\mathrm{ss}}=\left(-C_mA_m^{-1}B_m+D_m\right)r_0. \]

Unit command tracking is usually represented by \( M(0)=1 \). A different DC gain is legitimate when the command and output have different scaling or when the desired map is not unity.

3.4 Minimality and interpretable order

A controllable and observable realization avoids hidden unstable or redundant modes. The model order should be no larger than needed to express the desired behavior. Additional poles used only for command filtering must be placed intentionally and documented.

4.1 First-order model

A common scalar model is

\[ \dot{y}_m=-a_m y_m+a_m r, \qquad a_m>0, \]

\[ M(s)=\frac{a_m}{s+a_m}, \qquad \tau_m=\frac{1}{a_m}. \]

For a unit step and zero initial condition,

\[ y_m(t)=1-e^{-a_m t}. \]

The two-percent settling-time approximation follows by requiring \( e^{-a_m t_s}=0.02 \):

\[ t_s=\frac{-\ln(0.02)}{a_m}\approx\frac{3.912}{a_m}\approx 4\tau_m. \]

4.2 Unit-gain second-order model

The standard model used throughout this lesson is

\[ M(s)=\frac{\omega_n^2}{s^2+2\zeta\omega_n s+\omega_n^2}, \qquad \omega_n>0,\quad \zeta>0. \]

One controllable realization is

\[ A_m=\begin{bmatrix}0&1\\-\omega_n^2&-2\zeta\omega_n\end{bmatrix}, \quad B_m=\begin{bmatrix}0\\\omega_n^2\end{bmatrix}, \quad C_m=\begin{bmatrix}1&0\end{bmatrix},\quad D_m=0. \]

The characteristic roots are

\[ \lambda_{1,2}=-\zeta\omega_n\pm \omega_n\sqrt{\zeta^2-1}. \]

For \( 0<\zeta<1 \), define the damped natural frequency \( \omega_d=\omega_n\sqrt{1-\zeta^2} \). The poles become

\[ \lambda_{1,2}=-\zeta\omega_n\pm j\omega_d. \]

5.1 Unit-step response

For zero initial conditions and \( 0<\zeta<1 \), the inverse Laplace transform gives

\[ y_m(t)=1-e^{-\zeta\omega_n t} \left[\cos(\omega_d t)+\frac{\zeta}{\sqrt{1-\zeta^2}} \sin(\omega_d t)\right]. \]

5.2 Peak time

Differentiating the step response and setting its first positive extremum to zero gives \( \sin(\omega_d t_p)=0 \). The first peak therefore occurs at

\[ t_p=\frac{\pi}{\omega_d} =\frac{\pi}{\omega_n\sqrt{1-\zeta^2}}. \]

5.3 Percentage overshoot

Substituting \( t_p \) into the response yields the peak excess above one:

\[ M_p=\exp\left(-\frac{\pi\zeta}{\sqrt{1-\zeta^2}}\right), \qquad \%OS=100M_p. \]

Solving the overshoot equation for damping ratio gives a direct specification-to-parameter relation. With \( m=\%OS/100 \),

\[ \zeta=\frac{-\ln m}{\sqrt{\pi^2+(\ln m)^2}}. \]

5.4 Settling time and speed

The oscillatory envelope decays approximately as \( e^{-\zeta\omega_n t} \). For a tolerance \( \delta \),

\[ e^{-\zeta\omega_n t_s}\le\delta \quad\Longrightarrow\quad t_s\ge\frac{-\ln\delta}{\zeta\omega_n}. \]

Hence the familiar two-percent approximation is

\[ t_s\approx\frac{4}{\zeta\omega_n}. \]

Increasing \( \omega_n \) accelerates the model, whereas increasing \( \zeta \) suppresses oscillation. These parameters do not represent free performance: the plant, actuator, sampling, delay, and neglected dynamics constrain how large \( \omega_n \) can be selected.

6.1 Zero-input exponential stability

Let \( A_m \) be Hurwitz. For any symmetric \( Q=Q^T>0 \), the Lyapunov equation

\[ A_m^T P+PA_m=-Q \]

has a unique symmetric solution \( P=P^T>0 \). Choose \( V_m=x_m^T P x_m \). With \( r=0 \),

\[ \dot{V}_m=x_m^T(A_m^TP+PA_m)x_m=-x_m^TQx_m. \]

Using eigenvalue bounds,

\[ \lambda_{\min}(P)\|x_m\|^2\le V_m\le \lambda_{\max}(P)\|x_m\|^2, \]

\[ \dot V_m\le-\lambda_{\min}(Q)\|x_m\|^2 \le-\frac{\lambda_{\min}(Q)}{\lambda_{\max}(P)}V_m. \]

Therefore,

\[ V_m(t)\le V_m(0)e^{-ct}, \qquad c=\frac{\lambda_{\min}(Q)}{\lambda_{\max}(P)}, \]

\[ \|x_m(t)\|\le \sqrt{\frac{\lambda_{\max}(P)}{\lambda_{\min}(P)}} e^{-ct/2}\|x_m(0)\|. \]

6.2 Bounded-command result

For nonzero command,

\[ \dot V_m=-x_m^TQx_m+2x_m^TPB_m r. \]

Applying Cauchy-Schwarz and Young's inequality with any \( \varepsilon>0 \),

\[ 2\|x_m\|\|PB_m\||r| \le\varepsilon\|x_m\|^2+ \frac{\|PB_m\|^2}{\varepsilon}|r|^2. \]

Choose \( \varepsilon=\lambda_{\min}(Q)/2 \). Then

\[ \dot V_m\le- rac{\lambda_{\min}(Q)}{2}\|x_m\|^2 +\frac{2\|PB_m\|^2}{\lambda_{\min}(Q)}|r|^2. \]

Thus a bounded command produces a bounded model state, and the model has an input-to-state stability estimate. This property is indispensable: an MRAC objective cannot reasonably require the plant to follow an internally unstable desired trajectory.

A stable reference model can still be a poor MRAC specification. The model must be compatible with the known plant structure and with the matching conditions introduced previously.

7.1 Relative-degree compatibility

Let the plant and model transfer functions have relative degrees \( n_p-m_p \) and \( n_m-m_m \). A controller cannot generally make a causal plant exhibit an arbitrarily smaller relative degree without differentiating commands or cancelling uncertain dynamics. Therefore, classical MRAC formulations impose a known and compatible relative degree.

\[ n_p-m_p=n_m-m_m \quad\text{in the standard exact-matching formulation}. \]

7.2 Bandwidth separation

The reference bandwidth should lie below unmodeled flexible modes, neglected actuator poles, communication delays, and severe sensor-noise regions. A qualitative engineering rule is

\[ \omega_{\mathrm{model}}\ll \omega_{\mathrm{unmodeled\ dynamics}}, \]

where the separation factor must be chosen from available model validation, frequency-response data, and implementation constraints rather than from a universal numerical constant.

7.3 Control authority

A faster model demands larger acceleration, rate, and input. For the second-order model and a unit step, the initial desired acceleration is

\[ \ddot y_m(0^+)=\omega_n^2. \]

Doubling \( \omega_n \) quadruples this initial acceleration demand. If actuator saturation is ignored, the desired model may command behavior the physical system cannot reproduce. Saturation and adaptive modifications are treated later in the course.

7.4 Nonminimum-phase restrictions

Exact stable inversion of an uncertain right-half-plane zero is not admissible. Consequently, a reference model that implicitly demands cancellation of a nonminimum-phase plant zero is structurally infeasible. The desired response must retain the limitations imposed by known zero locations and delays.

8.1 Initial-condition consistency

Even with an ideal controller, a mismatch between plant and model initial states can create a transient tracking error. When measured plant states and model states have the same physical interpretation, choosing \( x_m(0)=x(0) \) may isolate parameter-induced tracking effects. When such initialization is impossible, the resulting homogeneous model transient should be recognized explicitly.

\[ x_m(t)=e^{A_m t}x_m(0)+ \int_0^t e^{A_m(t-\tau)}B_m r(\tau)\,d\tau. \]

8.2 Command smoothness

Step commands are useful for analysis, but they may be unsuitable for hardware with rate and acceleration limits. A command prefilter can be included as part of the desired model. For example,

\[ F_c(s)=\frac{\omega_f}{s+\omega_f}, \qquad M_{\mathrm{total}}(s)=M(s)F_c(s). \]

This increases the model order and changes its relative degree, so the controller structure and matching assumptions must remain consistent.

8.3 Open-loop and closed-loop reference models

The classical reference model evolves independently of tracking error. A later extension is a closed-loop reference model, in which the model receives an error-dependent feedback term to reshape adaptive transients. At this point, retain the classical form \( \dot x_m=A_mx_m+B_mr \); closed-loop reference models are mentioned only as an advanced research direction.

flowchart TD
  A["State the desired command-to-output behavior"] --> B["Choose model order and relative degree"]
  B --> C["Select DC gain and output scaling"]
  C --> D["Translate overshoot and settling specifications into poles"]
  D --> E["Verify Am is Hurwitz and realization is minimal"]
  E --> F["Check matching structure and known plant limitations"]
  F --> G["Check actuator, sampling, delay, and unmodeled-mode bandwidth"]
  G --> H["Simulate command classes and initial conditions"]
  H --> I["Freeze the reference model before designing adaptation"]
        

The final step is important. If reference-model parameters are repeatedly changed while adaptation gains are tuned, it becomes difficult to determine whether an observed transient is caused by the desired dynamics, by the controller structure, or by the parameter-update mechanism.

Compare two models with identical natural frequency \( \omega_n=4\ \mathrm{rad/s} \) but different damping ratios:

\[ M_1(s)=\frac{16}{s^2+5.6s+16},\qquad \zeta_1=0.7, \]

\[ M_2(s)=\frac{16}{s^2+2.4s+16},\qquad \zeta_2=0.3. \]

Theoretical percentage overshoots are

\[ \%OS_1=100\exp\left(-\frac{0.7\pi}{\sqrt{1-0.7^2}}\right) \approx 4.60\%, \]

\[ \%OS_2=100\exp\left(-\frac{0.3\pi}{\sqrt{1-0.3^2}}\right) \approx 37.23\%. \]

The approximate two-percent settling times are

\[ t_{s,1}\approx\frac{4}{0.7(4)}=1.429\ \mathrm{s},\qquad t_{s,2}\approx\frac{4}{0.3(4)}=3.333\ \mathrm{s}. \]

The comparison demonstrates that equal natural frequency does not imply equal transient quality. The lightly damped model asks the future adaptive loop to reproduce a large overshoot and a longer oscillatory transient. A reference model should encode behavior that is genuinely desired, not merely a convenient polynomial.

The Python implementation uses NumPy for vector operations and Matplotlib for visualization. The fourth-order Runge-Kutta integrator is implemented from scratch so that the numerical state update is visible. For larger control studies, scipy.signal, scipy.integrate.solve_ivp, and the python-control package are useful alternatives.

Chapter4_Lesson1.py

"""Chapter 4, Lesson 1: reference-model design and simulation.

Dependencies:
    numpy
    matplotlib

The program simulates a unit-DC-gain second-order reference model using an
explicit fourth-order Runge-Kutta integrator, evaluates standard step-response
metrics, and compares two damping-ratio choices.
"""

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 ReferenceModel:
    natural_frequency: float
    damping_ratio: float

    def matrices(self) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
        wn = self.natural_frequency
        zeta = self.damping_ratio
        a_matrix = np.array([[0.0, 1.0], [-wn * wn, -2.0 * zeta * wn]])
        b_matrix = np.array([0.0, wn * wn])
        c_matrix = np.array([1.0, 0.0])
        return a_matrix, b_matrix, c_matrix


def rk4_step(
    dynamics: Callable[[float, np.ndarray], np.ndarray],
    time: float,
    state: np.ndarray,
    step_size: float,
) -> np.ndarray:
    k1 = dynamics(time, state)
    k2 = dynamics(time + 0.5 * step_size, state + 0.5 * step_size * k1)
    k3 = dynamics(time + 0.5 * step_size, state + 0.5 * step_size * k2)
    k4 = dynamics(time + step_size, state + step_size * k3)
    return state + step_size * (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0


def simulate(
    model: ReferenceModel,
    command: Callable[[float], float],
    final_time: float = 8.0,
    step_size: float = 0.001,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    a_matrix, b_matrix, c_matrix = model.matrices()
    times = np.arange(0.0, final_time + step_size, step_size)
    states = np.zeros((times.size, 2), dtype=float)
    outputs = np.zeros(times.size, dtype=float)

    def dynamics(time: float, state: np.ndarray) -> np.ndarray:
        return a_matrix @ state + b_matrix * command(time)

    for index, time in enumerate(times[:-1]):
        outputs[index] = float(c_matrix @ states[index])
        states[index + 1] = rk4_step(dynamics, time, states[index], step_size)
    outputs[-1] = float(c_matrix @ states[-1])
    commands = np.array([command(time) for time in times])
    return times, outputs, commands


def step_metrics(
    times: np.ndarray,
    outputs: np.ndarray,
    target: float = 1.0,
    settling_band: float = 0.02,
) -> dict[str, float]:
    if target == 0.0:
        raise ValueError("The target must be nonzero for normalized metrics.")

    normalized = outputs / target
    peak_index = int(np.argmax(normalized))
    overshoot = max(0.0, (normalized[peak_index] - 1.0) * 100.0)

    above_10 = np.flatnonzero(normalized >= 0.10)
    above_90 = np.flatnonzero(normalized >= 0.90)
    rise_time = float("nan")
    if above_10.size and above_90.size:
        rise_time = float(times[above_90[0]] - times[above_10[0]])

    outside = np.flatnonzero(np.abs(normalized - 1.0) > settling_band)
    settling_time = 0.0 if outside.size == 0 else float(times[min(outside[-1] + 1, times.size - 1)])

    return {
        "rise_time_10_90_s": rise_time,
        "settling_time_2_percent_s": settling_time,
        "overshoot_percent": overshoot,
        "peak_time_s": float(times[peak_index]),
        "steady_state_value": float(outputs[-1]),
    }


def main() -> None:
    command = lambda _time: 1.0
    models = {
        "well_damped": ReferenceModel(natural_frequency=4.0, damping_ratio=0.70),
        "lightly_damped": ReferenceModel(natural_frequency=4.0, damping_ratio=0.30),
    }

    plt.figure(figsize=(9, 5))
    for label, model in models.items():
        times, outputs, commands = simulate(model, command)
        metrics = step_metrics(times, outputs)
        print(f"\n{label}: wn={model.natural_frequency:.2f}, zeta={model.damping_ratio:.2f}")
        for name, value in metrics.items():
            print(f"  {name}: {value:.6f}")
        plt.plot(times, outputs, label=label.replace("_", " "))

    plt.plot(times, commands, "--", label="command")
    plt.xlabel("Time (s)")
    plt.ylabel("Output")
    plt.title("Chapter 4 Lesson 1: Reference-model step responses")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig("Chapter4_Lesson1_reference_models.png", dpi=180)
    print("\nSaved Chapter4_Lesson1_reference_models.png")


if __name__ == "__main__":
    main()

This version uses only the C++17 standard library. It implements RK4, computes response metrics, and writes CSV data. In larger projects, Eigen is appropriate for matrix algebra and libraries such as Boost.Odeint can provide reusable numerical integrators.

Chapter4_Lesson1.cpp

// Chapter 4, Lesson 1: second-order reference-model simulation.
// Build: g++ -std=c++17 -O2 Chapter4_Lesson1.cpp -o Chapter4_Lesson1

#include <algorithm>
#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <vector>

struct ReferenceModel {
    double naturalFrequency;
    double dampingRatio;
};

using State = std::array<double, 2>;

State derivative(const ReferenceModel& model, const State& state, double command) {
    const double wn2 = model.naturalFrequency * model.naturalFrequency;
    return {state[1], -wn2 * state[0] - 2.0 * model.dampingRatio * model.naturalFrequency * state[1] + wn2 * command};
}

State addScaled(const State& a, const State& b, double scale) {
    return {a[0] + scale * b[0], a[1] + scale * b[1]};
}

State rk4Step(const ReferenceModel& model, const State& state, double command, double stepSize) {
    const State k1 = derivative(model, state, command);
    const State k2 = derivative(model, addScaled(state, k1, 0.5 * stepSize), command);
    const State k3 = derivative(model, addScaled(state, k2, 0.5 * stepSize), command);
    const State k4 = derivative(model, addScaled(state, k3, stepSize), command);

    return {
        state[0] + stepSize * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
        state[1] + stepSize * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0
    };
}

struct Metrics {
    double riseTime;
    double settlingTime;
    double overshootPercent;
    double peakTime;
    double steadyStateValue;
};

Metrics calculateMetrics(const std::vector<double>& times, const std::vector<double>& outputs) {
    if (times.size() != outputs.size() || times.empty()) {
        throw std::invalid_argument("Time and output vectors must have equal nonzero length.");
    }

    const auto peakIterator = std::max_element(outputs.begin(), outputs.end());
    const std::size_t peakIndex = static_cast<std::size_t>(std::distance(outputs.begin(), peakIterator));
    const double overshoot = std::max(0.0, (*peakIterator - 1.0) * 100.0);

    std::size_t index10 = outputs.size();
    std::size_t index90 = outputs.size();
    for (std::size_t i = 0; i < outputs.size(); ++i) {
        if (index10 == outputs.size() && outputs[i] >= 0.10) index10 = i;
        if (index90 == outputs.size() && outputs[i] >= 0.90) index90 = i;
    }
    const double riseTime = (index10 < outputs.size() && index90 < outputs.size())
        ? times[index90] - times[index10]
        : std::numeric_limits<double>::quiet_NaN();

    std::size_t lastOutside = 0;
    bool foundOutside = false;
    for (std::size_t i = 0; i < outputs.size(); ++i) {
        if (std::abs(outputs[i] - 1.0) > 0.02) {
            lastOutside = i;
            foundOutside = true;
        }
    }
    const std::size_t settlingIndex = foundOutside ? std::min(lastOutside + 1, times.size() - 1) : 0;

    return {riseTime, times[settlingIndex], overshoot, times[peakIndex], outputs.back()};
}

void runCase(const ReferenceModel& model, const std::string& name) {
    constexpr double stepSize = 0.001;
    constexpr double finalTime = 8.0;
    constexpr double command = 1.0;
    const std::size_t sampleCount = static_cast<std::size_t>(finalTime / stepSize) + 1;

    std::vector<double> times(sampleCount);
    std::vector<double> outputs(sampleCount);
    State state{0.0, 0.0};

    std::ofstream csv("Chapter4_Lesson1_" + name + ".csv");
    if (!csv) throw std::runtime_error("Could not create CSV output.");
    csv << "time,command,output,output_rate\n";

    for (std::size_t i = 0; i < sampleCount; ++i) {
        const double time = static_cast<double>(i) * stepSize;
        times[i] = time;
        outputs[i] = state[0];
        csv << std::setprecision(10) << time << ',' << command << ',' << state[0] << ',' << state[1] << '\n';
        if (i + 1 < sampleCount) state = rk4Step(model, state, command, stepSize);
    }

    const Metrics metrics = calculateMetrics(times, outputs);
    std::cout << "\n" << name << ": wn=" << model.naturalFrequency << ", zeta=" << model.dampingRatio << '\n'
              << "  rise time (10-90%): " << metrics.riseTime << " s\n"
              << "  settling time (2%): " << metrics.settlingTime << " s\n"
              << "  overshoot: " << metrics.overshootPercent << " %\n"
              << "  peak time: " << metrics.peakTime << " s\n"
              << "  final value: " << metrics.steadyStateValue << '\n';
}

int main() {
    try {
        runCase({4.0, 0.70}, "well_damped");
        runCase({4.0, 0.30}, "lightly_damped");
        std::cout << "\nCSV files created for both reference models.\n";
        return 0;
    } catch (const std::exception& exception) {
        std::cerr << "Error: " << exception.what() << '\n';
        return 1;
    }
}

The Java implementation uses records for immutable model and state data and the standard I/O library for CSV output. Apache Commons Math can replace the hand-coded integrator when general ODE solvers and linear algebra are needed.

Chapter4_Lesson1.java

// Chapter 4, Lesson 1: second-order reference-model simulation.
// Build and run: javac Chapter4_Lesson1.java && java Chapter4_Lesson1

import java.io.BufferedWriter;
import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.Locale;

public final class Chapter4_Lesson1 {
    private Chapter4_Lesson1() {}

    private record ReferenceModel(double naturalFrequency, double dampingRatio) {}
    private record State(double position, double velocity) {}
    private record Metrics(double riseTime, double settlingTime, double overshootPercent,
                           double peakTime, double steadyStateValue) {}

    private static State derivative(ReferenceModel model, State state, double command) {
        double wn2 = model.naturalFrequency() * model.naturalFrequency();
        return new State(
            state.velocity(),
            -wn2 * state.position()
                - 2.0 * model.dampingRatio() * model.naturalFrequency() * state.velocity()
                + wn2 * command
        );
    }

    private static State addScaled(State a, State b, double scale) {
        return new State(a.position() + scale * b.position(), a.velocity() + scale * b.velocity());
    }

    private static State rk4Step(ReferenceModel model, State state, double command, double stepSize) {
        State k1 = derivative(model, state, command);
        State k2 = derivative(model, addScaled(state, k1, 0.5 * stepSize), command);
        State k3 = derivative(model, addScaled(state, k2, 0.5 * stepSize), command);
        State k4 = derivative(model, addScaled(state, k3, stepSize), command);

        return new State(
            state.position() + stepSize * (k1.position() + 2.0 * k2.position() + 2.0 * k3.position() + k4.position()) / 6.0,
            state.velocity() + stepSize * (k1.velocity() + 2.0 * k2.velocity() + 2.0 * k3.velocity() + k4.velocity()) / 6.0
        );
    }

    private static Metrics calculateMetrics(double[] times, double[] outputs) {
        int peakIndex = 0;
        int index10 = -1;
        int index90 = -1;
        int lastOutside = -1;

        for (int i = 0; i < outputs.length; i++) {
            if (outputs[i] > outputs[peakIndex]) peakIndex = i;
            if (index10 < 0 && outputs[i] >= 0.10) index10 = i;
            if (index90 < 0 && outputs[i] >= 0.90) index90 = i;
            if (Math.abs(outputs[i] - 1.0) > 0.02) lastOutside = i;
        }

        double riseTime = (index10 >= 0 && index90 >= 0) ? times[index90] - times[index10] : Double.NaN;
        int settlingIndex = Math.min(lastOutside + 1, times.length - 1);
        double overshoot = Math.max(0.0, (outputs[peakIndex] - 1.0) * 100.0);
        return new Metrics(riseTime, times[settlingIndex], overshoot, times[peakIndex], outputs[outputs.length - 1]);
    }

    private static void runCase(ReferenceModel model, String name) throws IOException {
        final double stepSize = 0.001;
        final double finalTime = 8.0;
        final double command = 1.0;
        final int sampleCount = (int) Math.round(finalTime / stepSize) + 1;

        double[] times = new double[sampleCount];
        double[] outputs = new double[sampleCount];
        State state = new State(0.0, 0.0);
        Path csvPath = Path.of("Chapter4_Lesson1_" + name + ".csv");

        try (BufferedWriter writer = Files.newBufferedWriter(csvPath)) {
            writer.write("time,command,output,output_rate\n");
            for (int i = 0; i < sampleCount; i++) {
                double time = i * stepSize;
                times[i] = time;
                outputs[i] = state.position();
                writer.write(String.format(Locale.ROOT, "%.9f,%.3f,%.9f,%.9f%n",
                    time, command, state.position(), state.velocity()));
                if (i + 1 < sampleCount) state = rk4Step(model, state, command, stepSize);
            }
        }

        Metrics metrics = calculateMetrics(times, outputs);
        System.out.printf(Locale.ROOT,
            "%n%s: wn=%.2f, zeta=%.2f%n" +
            "  rise time (10-90%%): %.6f s%n" +
            "  settling time (2%%): %.6f s%n" +
            "  overshoot: %.6f %%%n" +
            "  peak time: %.6f s%n" +
            "  final value: %.6f%n",
            name, model.naturalFrequency(), model.dampingRatio(),
            metrics.riseTime(), metrics.settlingTime(), metrics.overshootPercent(),
            metrics.peakTime(), metrics.steadyStateValue());
    }

    public static void main(String[] args) {
        try {
            runCase(new ReferenceModel(4.0, 0.70), "well_damped");
            runCase(new ReferenceModel(4.0, 0.30), "lightly_damped");
            System.out.println("\nCSV files created for both reference models.");
        } catch (IOException exception) {
            System.err.println("I/O error: " + exception.getMessage());
            System.exit(1);
        }
    }
}

MATLAB uses ode45 for state-space simulation. When Control System Toolbox is available, the script also constructs a transfer function. When Simulink is installed, it programmatically creates a step source, a second-order transfer-function block, and a scope, then saves Chapter4_Lesson1_Simulink.slx.

Chapter4_Lesson1.m

%% Chapter4_Lesson1.m
% Reference Models and Desired Closed-Loop Behavior
% This script:
%   1) builds two unit-DC-gain second-order reference models,
%   2) simulates their step responses with ode45,
%   3) computes response metrics,
%   4) optionally creates an equivalent Simulink model.

clear; close all; clc;

models = [struct('wn',4.0,'zeta',0.70,'name','well damped'), ...
          struct('wn',4.0,'zeta',0.30,'name','lightly damped')];

tspan = [0 8];
x0 = [0; 0];
r = @(t) 1.0;

figure('Name','Chapter 4 Lesson 1 Reference Models');
hold on; grid on;

for k = 1:numel(models)
    wn = models(k).wn;
    zeta = models(k).zeta;
    A = [0 1; -wn^2 -2*zeta*wn];
    B = [0; wn^2];
    C = [1 0];

    ode = @(t,x) A*x + B*r(t);
    [t,x] = ode45(ode,tspan,x0);
    y = x*C.';

    metrics = localStepMetrics(t,y,1.0,0.02);
    fprintf('\n%s: wn=%.2f, zeta=%.2f\n',models(k).name,wn,zeta);
    disp(metrics);

    plot(t,y,'LineWidth',1.5,'DisplayName',models(k).name);
end

yline(1.0,'--','DisplayName','command');
xlabel('Time (s)'); ylabel('Output');
title('Reference-model step responses');
legend('Location','best');

% Control System Toolbox equivalent, when available:
if license('test','Control_Toolbox')
    wn = 4.0; zeta = 0.70;
    referenceTF = tf(wn^2,[1 2*zeta*wn wn^2]);
    fprintf('\nControl System Toolbox model:\n');
    disp(referenceTF);
    fprintf('DC gain = %.6f\n',dcgain(referenceTF));
end

% Create a small Simulink realization when Simulink is installed.
if license('test','Simulink')
    localBuildSimulinkModel(4.0,0.70);
end

function metrics = localStepMetrics(t,y,target,band)
    normalized = y/target;
    [peakValue,peakIndex] = max(normalized);
    index10 = find(normalized >= 0.10,1,'first');
    index90 = find(normalized >= 0.90,1,'first');
    if isempty(index10) || isempty(index90)
        riseTime = NaN;
    else
        riseTime = t(index90)-t(index10);
    end

    lastOutside = find(abs(normalized-1.0) > band,1,'last');
    if isempty(lastOutside)
        settlingTime = 0.0;
    else
        settlingTime = t(min(lastOutside+1,numel(t)));
    end

    metrics = table(riseTime,settlingTime,max(0,(peakValue-1)*100), ...
        t(peakIndex),y(end), ...
        'VariableNames',{'RiseTime10To90','SettlingTime2Percent', ...
        'OvershootPercent','PeakTime','SteadyStateValue'});
end

function localBuildSimulinkModel(wn,zeta)
    modelName = 'Chapter4_Lesson1_Simulink';
    if bdIsLoaded(modelName)
        close_system(modelName,0);
    end
    if exist([modelName '.slx'],'file')
        delete([modelName '.slx']);
    end

    new_system(modelName);
    open_system(modelName);

    add_block('simulink/Sources/Step',[modelName '/Command'], ...
        'Time','0','Before','0','After','1','Position',[40 80 70 110]);
    add_block('simulink/Continuous/Transfer Fcn',[modelName '/Reference Model'], ...
        'Numerator',mat2str(wn^2), ...
        'Denominator',mat2str([1 2*zeta*wn wn^2]), ...
        'Position',[135 70 275 120]);
    add_block('simulink/Sinks/Scope',[modelName '/Scope'], ...
        'Position',[340 75 370 105]);

    add_line(modelName,'Command/1','Reference Model/1','autorouting','on');
    add_line(modelName,'Reference Model/1','Scope/1','autorouting','on');
    set_param(modelName,'StopTime','8');
    save_system(modelName);
    fprintf('Created %s.slx\n',modelName);
end

The Mathematica notebook uses symbolic pole calculation and NDSolveValue for simulation. It also evaluates the analytical overshoot and settling-time formulas. The file is stored as a notebook expression and can be opened directly by Wolfram Mathematica.

Chapter4_Lesson1.nb

Notebook[{
 Cell["Chapter 4, Lesson 1: Reference Models and Desired Closed-Loop Behavior", "Title"],
 Cell["This notebook constructs a unit-DC-gain second-order reference model, proves its pole placement symbolically, simulates two damping choices, and evaluates standard step-response quantities.", "Text"],
 Cell[BoxData[
  RowBox[{"ClearAll", "[", "\"Global`*\"", "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"referenceTransfer", "[", RowBox[{"wn_", ",", "zeta_"}], "]"}], ":=",
   FractionBox[SuperscriptBox["wn", "2"],
    RowBox[{SuperscriptBox["s", "2"], "+", RowBox[{"2", " ", "zeta", " ", "wn", " ", "s"}], "+", SuperscriptBox["wn", "2"]}]]
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"poles", "=", RowBox[{"Solve", "[",
    RowBox[{RowBox[{SuperscriptBox["s", "2"], "+", RowBox[{"2", " ", "zeta", " ", "wn", " ", "s"}], "+", SuperscriptBox["wn", "2"]}], "==", "0"}], ",", "s", "]"}]}], ";",
   RowBox[{"FullSimplify", "[", "poles", ",", RowBox[{"0", "<", "zeta", "<", "1", "&&", "wn", ">", "0"}], "]"}]
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"simulate", "[", RowBox[{"wnValue_", ",", "zetaValue_"}], "]"}], ":=",
   RowBox[{"NDSolveValue", "[",
    RowBox[{
     RowBox[{"{",
      RowBox[{
       RowBox[{RowBox[{SuperscriptBox["y", "\[Prime]\[Prime]"], "[", "t", "]"}], "+",
        RowBox[{"2", " ", "zetaValue", " ", "wnValue", " ", RowBox[{SuperscriptBox["y", "\[Prime]"], "[", "t", "]"}]}], "+",
        RowBox[{SuperscriptBox["wnValue", "2"], " ", RowBox[{"y", "[", "t", "]"}]}], "==", SuperscriptBox["wnValue", "2"]}], ",",
       RowBox[{RowBox[{"y", "[", "0", "]"}], "==", "0"}], ",",
       RowBox[{RowBox[{SuperscriptBox["y", "\[Prime]"], "[", "0", "]"}], "==", "0"}]}
      ], "}"}], ",", "y", ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "8"}], "}"}]
    }], "]"}]
  }]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"wellDamped", "=", RowBox[{"simulate", "[", RowBox[{"4.0", ",", "0.70"}], "]"}]}], ";",
   RowBox[{"lightlyDamped", "=", RowBox[{"simulate", "[", RowBox[{"4.0", ",", "0.30"}], "]"}]}], ";"
  }]], "Input"],
 Cell[BoxData[
  RowBox[{"Plot", "[",
   RowBox[{
    RowBox[{"{", RowBox[{RowBox[{"wellDamped", "[", "t", "]"}], ",", RowBox[{"lightlyDamped", "[", "t", "]"}], ",", "1"}], "}"}], ",",
    RowBox[{"{", RowBox[{"t", ",", "0", ",", "8"}], "}"}], ",",
    RowBox[{"PlotLegends", "->", RowBox[{"{", RowBox[{"\"well damped\"", ",", "\"lightly damped\"", ",", "\"command\""}], "}"}]}], ",",
    RowBox[{"PlotRange", "->", "All"}], ",", RowBox[{"GridLines", "->", "Automatic"}], ",",
    RowBox[{"AxesLabel", "->", RowBox[{"{", RowBox[{"\"time (s)\"", ",", "\"output\""}], "}"}]}]
   }], "]"}]], "Input"],
 Cell[BoxData[
  RowBox[{
   RowBox[{"overshootPercent", "[", "zeta_", "]"}], ":=",
   RowBox[{"100", " ", RowBox[{"Exp", "[", RowBox[{"-", FractionBox[RowBox[{"Pi", " ", "zeta"}], SqrtBox[RowBox[{"1", "-", SuperscriptBox["zeta", "2"]}]]}], "]"}]}]
  }]], "Input"],
 Cell[BoxData[
  RowBox[{"Table", "[",
   RowBox[{
    RowBox[{"{", RowBox[{"z", ",", RowBox[{"overshootPercent", "[", "z", "]"}], ",", FractionBox["4", RowBox[{"z", " ", "4.0"}]]}], "}"}], ",",
    RowBox[{"{", RowBox[{"z", ",", RowBox[{"{", RowBox[{"0.30", ",", "0.70"}], "}"}]}], "}"}]
   }], "]"}]], "Input"],
 Cell["The last table reports damping ratio, theoretical percentage overshoot, and the common 2-percent settling-time approximation 4/(zeta wn).", "Text"]
}, WindowSize->{1100, 760}, FrontEndVersion->"14.0", StyleDefinitions->"Default.nb"]

Problem 1 (First-Order Specification): A desired model must reach and remain within two percent of a unit step in approximately \( 0.8\ \mathrm{s} \). Choose the pole of a unit-gain first-order reference model.

Solution: Use

\[ M(s)=\frac{a_m}{s+a_m}, \qquad t_s\approx\frac{4}{a_m}. \]

Therefore,

\[ a_m\approx\frac{4}{0.8}=5\ \mathrm{rad/s}. \]

The selected model is \( M(s)=5/(s+5) \). The exact two-percent time is \( -\ln(0.02)/5\approx0.782\ \mathrm{s} \).

Problem 2 (Second-Order Pole Selection): Select a unit-gain second-order reference model with approximately five-percent overshoot and a two-percent settling time of \( 2\ \mathrm{s} \).

Solution: Let \( m=0.05 \). Then

\[ \zeta=\frac{-\ln(0.05)}{\sqrt{\pi^2+[\ln(0.05)]^2}} \approx0.690. \]

From \( t_s\approx4/(\zeta\omega_n) \),

\[ \omega_n\approx\frac{4}{(0.690)(2)}\approx2.899\ \mathrm{rad/s}. \]

Thus,

\[ M(s)=\frac{8.404}{s^2+4.000s+8.404}. \]

The coefficients are rounded. A numerical step simulation should be used to verify the actual rise and settling times because textbook relations are approximations.

Problem 3 (Lyapunov Equation): Consider the scalar model \( \dot x_m=-3x_m+r \). For \( Q=2 \), solve the Lyapunov equation and derive a bound for zero input.

Solution: The scalar Lyapunov equation is

\[ (-3)P+P(-3)=-2, \]

so \( P=1/3 \). With \( V=Px_m^2 \),

\[ \dot V=-2x_m^2=-6V. \]

Hence

\[ V(t)=V(0)e^{-6t}, \qquad |x_m(t)|=|x_m(0)|e^{-3t}. \]

The Lyapunov result exactly recovers the scalar solution and confirms exponential stability.

Problem 4 (Feasibility and Control Demand): Two candidate second-order models use the same damping ratio \( \zeta=0.7 \), but natural frequencies \( 4 \) and \( 12\ \mathrm{rad/s} \). Compare their initial desired accelerations for a unit step.

Solution: Since

\[ \ddot y_m(0^+)=\omega_n^2, \]

the two demands are

\[ \ddot y_{m,1}(0^+)=16,\qquad \ddot y_{m,2}(0^+)=144. \]

The faster model demands nine times the initial acceleration. Even though both models are stable, the second may be incompatible with actuator force, current, torque, or rate limits.

Problem 5 (Steady-State Gain): Let

\[ A_m=\begin{bmatrix}0&1\\-9&-4.2\end{bmatrix},\quad B_m=\begin{bmatrix}0\\18\end{bmatrix},\quad C_m=\begin{bmatrix}1&0\end{bmatrix},\quad D_m=0. \]

Determine the DC gain from \( r \) to \( y_m \).

Solution: The equilibrium gain is

\[ M(0)=-C_mA_m^{-1}B_m+D_m. \]

Equivalently, the transfer function is

\[ M(s)=\frac{18}{s^2+4.2s+9}. \]

Therefore,

\[ M(0)=\frac{18}{9}=2. \]

A unit command produces a desired steady output of two. To obtain unit DC gain, replace the lower entry of \( B_m \) by nine.

Problem 6 (Relative-Degree Check): A known plant structure has relative degree two, while a proposed model is \( M(s)=5/(s+5) \). Explain the structural concern.

Solution:

The model has relative degree one. Exact matching would require the closed-loop plant to respond with one fewer net integration than its known input-output structure permits. Achieving this may require noncausal command differentiation or cancellation of uncertain dynamics. A conventional MRAC design should instead choose a compatible relative-degree-two model, unless a specifically derived augmented controller structure establishes otherwise.

The reference model defines the desired command-to-output dynamics before any adaptive law is introduced. A valid model is stable, causal, minimally realized, correctly scaled, and structurally compatible with the plant. Its poles translate performance specifications into mathematical dynamics, but its bandwidth must respect control authority, delays, sampling, noise, and unmodeled modes. Lyapunov analysis proves that a Hurwitz model is exponentially stable and bounded under bounded commands. In the next lesson, the reference-model state and plant state will be combined to derive tracking error and parameter-error dynamics.