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

• construct a stable unity-DC-gain second-order reference model;
• map damping ratio and natural frequency to poles and transient specifications;
• derive overshoot, peak-time, rise-time, settling-time, and bandwidth relations;
• translate performance specifications into admissible regions for model parameters;
• explain why excessive reference-model speed increases control demand and robustness risk;
• screen candidate models before deriving the adaptive law in the next lessons.

The plant tracking problem was established in Lesson 1. Here, no parameter-update law is designed. Instead, we define the target dynamics that the later MRAC law will attempt to reproduce.

For a scalar command \( r(t) \) and desired output \( y_m(t) \), use the normalized model

\[ G_m(s)=\frac{Y_m(s)}{R(s)}= \frac{\omega_n^2}{s^2+2\zeta\omega_n s+\omega_n^2}, \qquad \omega_n > 0,\; \zeta > 0. \]

The numerator makes the DC gain equal to one, so a constant command is tracked without steady-state scaling by the reference model:

\[ G_m(0)=1, \qquad \lim_{t→\infty}y_m(t)=r_0 \quad \text{for} \quad r(t)=r_0\mathbf{1}(t). \]

The equivalent differential equation is

\[ \ddot y_m+2\zeta\omega_n\dot y_m+\omega_n^2 y_m =\omega_n^2 r. \]

With \( x_m=[y_m\;\;\dot y_m]^T \), a convenient realization is

\[ \dot x_m=A_mx_m+B_mr,\qquad y_m=C_mx_m, \] \[ 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}. \]

The characteristic roots are

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

Because \( \omega_n > 0 \) and \( \zeta > 0 \), both roots lie in the open left half-plane. Thus \( A_m \) is Hurwitz. This stability property is mandatory: an adaptive controller must not be asked to reproduce an unstable desired model.

flowchart TD
  R["Command r(t)"] --> M["Second-order reference model"]
  M --> YM["Desired output y_m(t)"]
  YM --> E["Tracking comparison with plant output y(t)"]
  P["Design variables: zeta and omega_n"] --> M
  S["Specifications and feasibility limits"] --> P
        

For \( 0 < \zeta < 1 \), define

\[ \omega_d=\omega_n\sqrt{1-\zeta^2}, \qquad \sigma=\zeta\omega_n. \]

The poles become

\[ s_{1,2}=-\sigma\pm j\omega_d. \]

The damping ratio controls the pole angle and therefore oscillation and overshoot. The product \( \zeta\omega_n \) controls the exponential envelope. Increasing \( \omega_n \) at fixed \( \zeta \) scales the entire pole pair away from the origin, speeding up all model dynamics without changing the normalized response shape.

\[ y_m(t;\zeta,\alpha\omega_n)=y_m(\alpha t;\zeta,\omega_n) \quad \text{for a unit step and zero initial state}. \]

The scaling identity follows by introducing \( \tau=\omega_n t \) in the normalized equation:

\[ \frac{d^2y_m}{d\tau^2}+2\zeta\frac{dy_m}{d\tau}+y_m=r. \]

Therefore \( \zeta \) selects shape, while \( \omega_n \) selects time scale.

For a unit step, zero initial conditions, and \( 0 < \zeta < 1 \), the response is

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

Peak time. Differentiating and setting the velocity to zero gives

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

Maximum fractional overshoot. Evaluating the response at the first peak gives

\[ M_p=e^{-\pi\zeta/\sqrt{1-\zeta^2}}, \qquad \%OS=100M_p. \]

Solving this equation for the minimum damping ratio associated with an allowed overshoot \( 0 < M_{p,\max} < 1 \) yields

\[ \zeta_{\min}= \frac{-\ln M_{p,\max}}{\sqrt{\pi^2+(\ln M_{p,\max})^2}}. \]

Zero-to-one-hundred-percent rise time. The first crossing of unity is

\[ t_r=\frac{\pi-\cos^{-1}(\zeta)} {\omega_n\sqrt{1-\zeta^2}}. \]

Software tools often report a 10–90% rise time instead. That quantity is usually computed numerically because a compact elementary expression is not available for arbitrary damping ratio.

Two-percent settling time. The standard engineering approximation is

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

It follows from requiring the exponential envelope to become small. A more explicit envelope condition is

\[ \frac{e^{-\zeta\omega_n t_s}}{\sqrt{1-\zeta^2}} \leq 0.02, \]

which gives

\[ t_s\geq \frac{-\ln\left(0.02\sqrt{1-\zeta^2}\right)}{\zeta\omega_n}. \]

The familiar coefficient 4 is a useful approximation over common design values, especially near \( \zeta\approx0.7 \); it is not an exact identity.

The squared magnitude of the reference model is

\[ |G_m(j\omega)|^2= \frac{\omega_n^4}{(\omega_n^2-\omega^2)^2+ 4\zeta^2\omega_n^2\omega^2}. \]

Defining the minus-three-decibel bandwidth by \( |G_m(j\omega_{bw})|^2=1/2 \), and setting \( q=(\omega_{bw}/\omega_n)^2 \), gives

\[ q^2+(4\zeta^2-2)q-1=0. \]

The positive root produces

\[ \omega_{bw}=\omega_n \sqrt{1-2\zeta^2+\sqrt{2-4\zeta^2+4\zeta^4}}. \]

Hence model bandwidth is proportional to \( \omega_n \), with a damping-dependent proportionality factor. A reference model that has bandwidth close to neglected plant modes, actuator dynamics, sampling limitations, or noisy sensor content asks the adaptive loop to operate where the nominal model is least trustworthy.

Consider the second-order plant form introduced in Lesson 1:

\[ \ddot y+a_1\dot y+a_0y=bu, \qquad b\neq0. \]

If the plant were exactly on the desired trajectory, the control required to maintain \( y=y_m \) would satisfy

\[ u^*=\frac{1}{b}\left(\ddot y_m+a_1\dot y_m+a_0y_m\right). \]

Substitution of the reference-model equation gives

\[ u^*=\frac{1}{b}\left[ (a_1-2\zeta\omega_n)\dot y_m+(a_0-\omega_n^2)y_m+ \omega_n^2r\right]. \]

This expression is not yet an adaptive control law; it is a feasibility diagnostic. It shows that the command-related term grows with \( \omega_n^2 \). For a unit step and zero model state,

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

Doubling \( \omega_n \) therefore quadruples the initial desired acceleration. In a real loop this tends to increase peak control, excite neglected dynamics, amplify measurement noise through high-bandwidth feedback, and encounter actuator rate or magnitude limits. The best reference model is consequently constrained by physical feasibility, not only by desired settling time.

1. Speed versus control effort: increasing \( \omega_n \) reduces rise and settling times, but raises desired acceleration and normally raises peak control effort.
2. Speed versus modelling validity: a high-bandwidth model may lie close to actuator poles, flexible modes, transport delays, or unmodelled zeros.
3. Damping versus rise time: increasing \( \zeta \) suppresses overshoot, but an overdamped model can be slower in its early response for fixed \( \omega_n \).
4. Tracking sharp commands versus noise sensitivity: a fast model preserves more high-frequency command content; the resulting loop must react to higher-frequency measurement components as well.
5. Ideal specifications versus actuator constraints: model output, velocity, acceleration, and optionally jerk should be checked against physical limits before adaptation is designed.

A damping ratio around 0.7 is a common starting point because it gives modest overshoot and useful speed, but it is not a universal optimum. Monotonic response may require \( \zeta\geq1 \); stringent speed may justify a lower value if overshoot is acceptable and the plant can support it.

Reference-model selection can be posed as a feasibility problem. Let the specifications provide a lower bound on damping and a lower bound on natural frequency, while implementation constraints provide an upper bound on natural frequency:

\[ \zeta\geq\zeta_{\min}, \qquad \omega_n\geq\frac{4}{\zeta t_{s,\max}}, \qquad \omega_n\leq\omega_{n,\max}^{\text{physical}}. \]

The upper bound should be obtained from actuator testing, known plant bandwidth, the first neglected mode, sampling rate, sensor filtering, and permitted command acceleration. A conservative design leaves clear separation from these limits. There is no universal numerical separation factor valid for every plant.

flowchart TD
  A["Specify overshoot, settling time, command class"] --> B["Compute minimum zeta from overshoot"]
  B --> C["Compute minimum omega_n from settling time"]
  C --> D["Check plant and actuator bandwidth limits"]
  D --> E["Check peak model velocity and acceleration"]
  E --> F{"Candidate feasible?"}
  F -->|"No"| G["Relax speed, reshape command, \nor improve hardware"]
  G --> C
  F -->|"Yes"| H["Simulate candidate reference models"]
  H --> I["Select model and freeze it for adaptive-law design"]
        

Suppose the design requirements are:

• maximum step overshoot: \( 5\% \);
• two-percent settling time: \( t_s\leq1.5\;\text{s} \);
• second-order unity-gain reference model.

The overshoot condition gives

\[ \zeta_{\min}= \frac{-\ln(0.05)}{\sqrt{\pi^2+[\ln(0.05)]^2}} \approx0.6901. \]

Using the standard two-percent settling approximation,

\[ \omega_n\geq\frac{4}{(0.6901)(1.5)} \approx3.864\;\text{rad/s}. \]

A convenient candidate is

\[ \zeta=0.70, \qquad \omega_n=4\;\text{rad/s}, \] \[ G_m(s)=\frac{16}{s^2+5.6s+16}. \]

Its approximate and exact-model indicators are

\[ t_s\approx\frac{4}{(0.70)(4)}=1.429\;\text{s}, \] \[ M_p=e^{-\pi(0.70)/\sqrt{1-0.70^2}}\approx0.0460, \] \[ \%OS\approx4.60\%, \qquad \ddot y_m(0^+)=16\;\text{units/s}^2. \]

The transient requirements are met. The candidate is accepted only if the plant and actuator can safely produce the associated velocity, acceleration, bandwidth, and control demand. If not, the settling-time requirement must be relaxed or the command must be shaped.

The implementations compare four candidate models. Each reports overshoot, 10–90% rise time, two-percent settling time, and peak reference-model acceleration. The acceleration metric is a useful proxy for aggressiveness before a complete adaptive controller is available.

• Python: NumPy and Matplotlib are used directly; the control package is a standard option for transfer-function and state-space studies.
• C++: the example uses a from-scratch RK4 solver; Eigen and Boost.Odeint are useful for larger adaptive-control simulations.
• Java: the example uses standard Java; Apache Commons Math provides matrix and ODE facilities.
• MATLAB/Simulink: Control System Toolbox provides tf, step, and stepinfo; the script also creates a Simulink realization programmatically.
• Wolfram Mathematica: TransferFunctionModel, OutputResponse, symbolic algebra, and numerical sampling are built in.

Chapter6_Lesson2.py

# Chapter6_Lesson2.py
# Reference-model selection for a second-order MRAC design.
# Dependencies: numpy, matplotlib. Optional library for further work: python-control.

from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path
from typing import Callable

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class ModelResult:
    zeta: float
    omega_n: float
    time: np.ndarray
    output: np.ndarray
    velocity: np.ndarray
    acceleration: np.ndarray
    overshoot_percent: float
    rise_time_10_90: float
    settling_time_2_percent: float
    peak_acceleration: float


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


def first_crossing(time: np.ndarray, signal: np.ndarray, level: float) -> float:
    indices = np.flatnonzero(signal >= level)
    return float(time[indices[0]]) if indices.size else float("nan")


def settling_time(time: np.ndarray, signal: np.ndarray, tolerance: float = 0.02) -> float:
    violations = np.flatnonzero(np.abs(signal - 1.0) > tolerance)
    if not violations.size:
        return 0.0
    last = int(violations[-1])
    return float(time[last + 1]) if last + 1 < time.size else float("nan")


def simulate_reference_model(
    zeta: float,
    omega_n: float,
    duration: float = 6.0,
    step: float = 0.001,
) -> ModelResult:
    if zeta <= 0.0 or omega_n <= 0.0:
        raise ValueError("zeta and omega_n must be positive")

    time = np.arange(0.0, duration + step, step)
    state = np.zeros(2, dtype=float)  # state = [y_m, y_m_dot]
    output = np.zeros_like(time)
    velocity = np.zeros_like(time)
    acceleration = np.zeros_like(time)

    def dynamics(_: float, x: np.ndarray) -> np.ndarray:
        reference = 1.0
        return np.array(
            [x[1], -2.0 * zeta * omega_n * x[1] - omega_n**2 * x[0] + omega_n**2 * reference],
            dtype=float,
        )

    for index, current_time in enumerate(time):
        output[index] = state[0]
        velocity[index] = state[1]
        acceleration[index] = dynamics(current_time, state)[1]
        if index + 1 < time.size:
            state = rk4_step(state, current_time, step, dynamics)

    peak = float(np.max(output))
    overshoot = max(0.0, (peak - 1.0) * 100.0)
    t10 = first_crossing(time, output, 0.10)
    t90 = first_crossing(time, output, 0.90)
    rise = t90 - t10 if np.isfinite(t10) and np.isfinite(t90) else float("nan")

    return ModelResult(
        zeta=zeta,
        omega_n=omega_n,
        time=time,
        output=output,
        velocity=velocity,
        acceleration=acceleration,
        overshoot_percent=overshoot,
        rise_time_10_90=rise,
        settling_time_2_percent=settling_time(time, output),
        peak_acceleration=float(np.max(np.abs(acceleration))),
    )


def main() -> None:
    candidates = [(0.50, 4.0), (0.70, 4.0), (1.00, 4.0), (0.70, 7.0)]
    results = [simulate_reference_model(zeta, omega_n) for zeta, omega_n in candidates]

    print(" zeta   omega_n   overshoot(%)   rise_10_90(s)   settle_2%(s)   peak_accel")
    for result in results:
        print(
            f" {result.zeta:4.2f}    {result.omega_n:6.2f}       "
            f"{result.overshoot_percent:8.3f}        {result.rise_time_10_90:8.3f}       "
            f"{result.settling_time_2_percent:8.3f}      {result.peak_acceleration:8.3f}"
        )

    for result in results:
        label = f"zeta={result.zeta:.2f}, omega_n={result.omega_n:.1f} rad/s"
        plt.plot(result.time, result.output, label=label)

    plt.axhline(1.0, linewidth=1.0)
    plt.xlabel("Time (s)")
    plt.ylabel("Reference-model output y_m")
    plt.title("Second-order reference-model trade-offs")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    figure_path = Path(__file__).with_name("Chapter6_Lesson2.png")
    plt.savefig(figure_path, dpi=180)
    print(f"Saved figure: {figure_path}")
    plt.show()


if __name__ == "__main__":
    main()

Chapter6_Lesson2.cpp

// Chapter6_Lesson2.cpp
// Standard-library RK4 study of second-order MRAC reference models.
// Optional libraries for larger projects: Eigen and Boost.Odeint.

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

struct Result {
    double zeta{};
    double omega_n{};
    double overshoot_percent{};
    double rise_time_10_90{};
    double settling_time_2_percent{};
    double peak_acceleration{};
    std::vector<double> time;
    std::vector<double> output;
};

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

State derivative(const State& x, double zeta, double omega_n) {
    const double reference = 1.0;
    return {x[1], -2.0 * zeta * omega_n * x[1]
                      - omega_n * omega_n * x[0]
                      + omega_n * omega_n * reference};
}

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

State rk4_step(const State& x, double step, double zeta, double omega_n) {
    const State k1 = derivative(x, zeta, omega_n);
    const State k2 = derivative(add(x, k1, 0.5 * step), zeta, omega_n);
    const State k3 = derivative(add(x, k2, 0.5 * step), zeta, omega_n);
    const State k4 = derivative(add(x, k3, step), zeta, omega_n);
    return {
        x[0] + step * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
        x[1] + step * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0
    };
}

double first_crossing(const std::vector<double>& time,
                      const std::vector<double>& signal,
                      double level) {
    for (std::size_t i = 0; i < signal.size(); ++i) {
        if (signal[i] >= level) {
            return time[i];
        }
    }
    return std::numeric_limits<double>::quiet_NaN();
}

double settling_time(const std::vector<double>& time,
                     const std::vector<double>& signal,
                     double tolerance = 0.02) {
    std::size_t last_violation = 0;
    bool found = false;
    for (std::size_t i = 0; i < signal.size(); ++i) {
        if (std::abs(signal[i] - 1.0) > tolerance) {
            last_violation = i;
            found = true;
        }
    }
    if (!found) {
        return 0.0;
    }
    return (last_violation + 1 < time.size())
        ? time[last_violation + 1]
        : std::numeric_limits<double>::quiet_NaN();
}

Result simulate(double zeta, double omega_n,
                double duration = 6.0, double step = 0.001) {
    if (zeta <= 0.0 || omega_n <= 0.0) {
        throw std::invalid_argument("zeta and omega_n must be positive");
    }

    const std::size_t count = static_cast<std::size_t>(std::floor(duration / step)) + 1;
    Result result;
    result.zeta = zeta;
    result.omega_n = omega_n;
    result.time.resize(count);
    result.output.resize(count);

    State state{0.0, 0.0};
    double peak = 0.0;
    double peak_acceleration = 0.0;

    for (std::size_t i = 0; i < count; ++i) {
        const double time = static_cast<double>(i) * step;
        const State dx = derivative(state, zeta, omega_n);
        result.time[i] = time;
        result.output[i] = state[0];
        peak = std::max(peak, state[0]);
        peak_acceleration = std::max(peak_acceleration, std::abs(dx[1]));
        if (i + 1 < count) {
            state = rk4_step(state, step, zeta, omega_n);
        }
    }

    result.overshoot_percent = std::max(0.0, (peak - 1.0) * 100.0);
    const double t10 = first_crossing(result.time, result.output, 0.10);
    const double t90 = first_crossing(result.time, result.output, 0.90);
    result.rise_time_10_90 = t90 - t10;
    result.settling_time_2_percent = settling_time(result.time, result.output);
    result.peak_acceleration = peak_acceleration;
    return result;
}

int main() {
    const std::vector<std::pair<double, double>> candidates{
        {0.50, 4.0}, {0.70, 4.0}, {1.00, 4.0}, {0.70, 7.0}
    };

    std::vector<Result> results;
    for (const auto& candidate : candidates) {
        results.push_back(simulate(candidate.first, candidate.second));
    }

    std::cout << " zeta   omega_n   overshoot(%)   rise_10_90(s)   settle_2%(s)   peak_accel\n";
    std::cout << std::fixed << std::setprecision(3);
    for (const Result& result : results) {
        std::cout << std::setw(5) << result.zeta << "   "
                  << std::setw(7) << result.omega_n << "       "
                  << std::setw(9) << result.overshoot_percent << "        "
                  << std::setw(9) << result.rise_time_10_90 << "       "
                  << std::setw(9) << result.settling_time_2_percent << "      "
                  << std::setw(9) << result.peak_acceleration << '\n';
    }

    std::ofstream csv("Chapter6_Lesson2.csv");
    if (!csv) {
        std::cerr << "Could not create Chapter6_Lesson2.csv\n";
        return 1;
    }
    csv << "time";
    for (const Result& result : results) {
        csv << ",zeta_" << result.zeta << "_wn_" << result.omega_n;
    }
    csv << '\n';
    for (std::size_t i = 0; i < results.front().time.size(); ++i) {
        csv << results.front().time[i];
        for (const Result& result : results) {
            csv << ',' << result.output[i];
        }
        csv << '\n';
    }
    std::cout << "Saved data: Chapter6_Lesson2.csv\n";
    return 0;
}

Chapter6_Lesson2.java

// Chapter6_Lesson2.java
// Pure-Java RK4 study of second-order MRAC reference models.
// Optional library for larger projects: Apache Commons Math.

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

public final class Chapter6_Lesson2 {
    private Chapter6_Lesson2() { }

    private record Result(
            double zeta,
            double omegaN,
            double[] time,
            double[] output,
            double overshootPercent,
            double riseTime1090,
            double settlingTime2Percent,
            double peakAcceleration) { }

    private static double[] derivative(double[] state, double zeta, double omegaN) {
        double reference = 1.0;
        return new double[] {
            state[1],
            -2.0 * zeta * omegaN * state[1]
                - omegaN * omegaN * state[0]
                + omegaN * omegaN * reference
        };
    }

    private static double[] shifted(double[] state, double[] slope, double scale) {
        return new double[] {
            state[0] + scale * slope[0],
            state[1] + scale * slope[1]
        };
    }

    private static double[] rk4Step(double[] state, double step, double zeta, double omegaN) {
        double[] k1 = derivative(state, zeta, omegaN);
        double[] k2 = derivative(shifted(state, k1, 0.5 * step), zeta, omegaN);
        double[] k3 = derivative(shifted(state, k2, 0.5 * step), zeta, omegaN);
        double[] k4 = derivative(shifted(state, k3, step), zeta, omegaN);
        return new double[] {
            state[0] + step * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
            state[1] + step * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0
        };
    }

    private static double firstCrossing(double[] time, double[] signal, double level) {
        for (int i = 0; i < signal.length; i++) {
            if (signal[i] >= level) {
                return time[i];
            }
        }
        return Double.NaN;
    }

    private static double settlingTime(double[] time, double[] signal, double tolerance) {
        int lastViolation = -1;
        for (int i = 0; i < signal.length; i++) {
            if (Math.abs(signal[i] - 1.0) > tolerance) {
                lastViolation = i;
            }
        }
        if (lastViolation < 0) {
            return 0.0;
        }
        return lastViolation + 1 < time.length ? time[lastViolation + 1] : Double.NaN;
    }

    private static Result simulate(double zeta, double omegaN, double duration, double step) {
        if (zeta <= 0.0 || omegaN <= 0.0) {
            throw new IllegalArgumentException("zeta and omegaN must be positive");
        }

        int count = (int) Math.floor(duration / step) + 1;
        double[] time = new double[count];
        double[] output = new double[count];
        double[] state = {0.0, 0.0};
        double peak = 0.0;
        double peakAcceleration = 0.0;

        for (int i = 0; i < count; i++) {
            time[i] = i * step;
            output[i] = state[0];
            double acceleration = derivative(state, zeta, omegaN)[1];
            peak = Math.max(peak, state[0]);
            peakAcceleration = Math.max(peakAcceleration, Math.abs(acceleration));
            if (i + 1 < count) {
                state = rk4Step(state, step, zeta, omegaN);
            }
        }

        double t10 = firstCrossing(time, output, 0.10);
        double t90 = firstCrossing(time, output, 0.90);
        return new Result(
            zeta,
            omegaN,
            time,
            output,
            Math.max(0.0, (peak - 1.0) * 100.0),
            t90 - t10,
            settlingTime(time, output, 0.02),
            peakAcceleration
        );
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);
        double[][] candidates = {
            {0.50, 4.0}, {0.70, 4.0}, {1.00, 4.0}, {0.70, 7.0}
        };

        List<Result> results = new ArrayList<>();
        for (double[] candidate : candidates) {
            results.add(simulate(candidate[0], candidate[1], 6.0, 0.001));
        }

        System.out.println(" zeta   omega_n   overshoot(%)   rise_10_90(s)   settle_2%(s)   peak_accel");
        for (Result result : results) {
            System.out.printf(
                " %4.2f    %6.2f       %8.3f        %8.3f       %8.3f      %8.3f%n",
                result.zeta(), result.omegaN(), result.overshootPercent(),
                result.riseTime1090(), result.settlingTime2Percent(),
                result.peakAcceleration()
            );
        }

        Path csvPath = Path.of("Chapter6_Lesson2.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(csvPath)) {
            writer.write("time");
            for (Result result : results) {
                writer.write(",zeta_" + result.zeta() + "_wn_" + result.omegaN());
            }
            writer.newLine();
            for (int i = 0; i < results.get(0).time().length; i++) {
                writer.write(Double.toString(results.get(0).time()[i]));
                for (Result result : results) {
                    writer.write("," + result.output()[i]);
                }
                writer.newLine();
            }
        }
        System.out.println("Saved data: " + csvPath.toAbsolutePath());
    }
}

Chapter6_Lesson2.m

% Chapter6_Lesson2.m
% Reference-model selection for a second-order MRAC design.
% Uses Control System Toolbox. The final section optionally creates a
% programmatic Simulink model when Simulink is installed.

clear; clc; close all;

candidates = [0.50 4.0; 0.70 4.0; 1.00 4.0; 0.70 7.0];
t = (0:0.001:6).';
results = table('Size', [size(candidates, 1), 6], ...
    'VariableTypes', repmat("double", 1, 6), ...
    'VariableNames', {'Zeta','OmegaN','OvershootPercent', ...
                      'RiseTime1090','SettlingTime2Percent','PeakAcceleration'});

figure('Name', 'Chapter 6 Lesson 2: Reference-model trade-offs');
hold on; grid on;

for k = 1:size(candidates, 1)
    zeta = candidates(k, 1);
    omegaN = candidates(k, 2);

    Gm = tf(omegaN^2, [1, 2*zeta*omegaN, omegaN^2]);
    [y, tout] = step(Gm, t);
    info = stepinfo(y, tout, 1, 'SettlingTimeThreshold', 0.02, ...
                    'RiseTimeLimits', [0.10 0.90]);

    velocity = gradient(y, tout);
    acceleration = gradient(velocity, tout);

    results.Zeta(k) = zeta;
    results.OmegaN(k) = omegaN;
    results.OvershootPercent(k) = info.Overshoot;
    results.RiseTime1090(k) = info.RiseTime;
    results.SettlingTime2Percent(k) = info.SettlingTime;
    results.PeakAcceleration(k) = max(abs(acceleration));

    plot(tout, y, 'DisplayName', ...
        sprintf('zeta=%.2f, omega_n=%.1f rad/s', zeta, omegaN));
end

yline(1.0, '--', 'Unit command');
xlabel('Time (s)');
ylabel('Reference-model output y_m');
title('Second-order reference-model trade-offs');
legend('Location', 'best');
exportgraphics(gcf, 'Chapter6_Lesson2.png', 'Resolution', 180);
disp(results);
writetable(results, 'Chapter6_Lesson2_metrics.csv');

%% Specification-driven design example
maximumOvershoot = 0.05;
maximumSettlingTime = 1.5;
zetaMinimum = -log(maximumOvershoot) / ...
    sqrt(pi^2 + log(maximumOvershoot)^2);
omegaNMinimum = 4 / (zetaMinimum * maximumSettlingTime);
fprintf('Minimum zeta from 5%% overshoot: %.4f\n', zetaMinimum);
fprintf('Minimum omega_n from 1.5 s settling: %.4f rad/s\n', omegaNMinimum);

%% Optional: create a Simulink realization of one selected reference model
if license('test', 'Simulink')
    modelName = 'Chapter6_Lesson2_Simulink';
    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    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', 'omegaN^2', ...
        'Denominator', '[1 2*zeta*omegaN omegaN^2]', ...
        'Position', [130 70 270 120]);
    add_block('simulink/Sinks/Scope', [modelName '/Scope'], ...
        'Position', [330 75 360 115]);
    add_line(modelName, 'Command/1', 'Reference Model/1');
    add_line(modelName, 'Reference Model/1', 'Scope/1');

    zeta = 0.70; %#ok<NASGU>
    omegaN = 4.0; %#ok<NASGU>
    set_param(modelName, 'StopTime', '6');
    save_system(modelName, [modelName '.slx']);
    fprintf('Saved Simulink model: %s.slx\n', modelName);
else
    fprintf('Simulink is unavailable; skipped programmatic model creation.\n');
end

Chapter6_Lesson2.nb

Notebook[{
  Cell["Chapter 6, Lesson 2: Reference Model Selection and Performance Trade-Offs", "Title"],
  Cell["Evaluate the following cell to compare second-order reference models.", "Text"],
  Cell[BoxData[
"ClearAll[\"Global`*\"];

candidates = {{0.50, 4.0}, {0.70, 4.0}, {1.00, 4.0}, {0.70, 7.0}};

referenceModel[zeta_, omegaN_] :=
  TransferFunctionModel[{{omegaN^2}}, {1, 2 zeta omegaN, omegaN^2}, s];

responses = Table[
  With[{zeta = candidate[[1]], omegaN = candidate[[2]]},
    {zeta, omegaN,
     OutputResponse[referenceModel[zeta, omegaN], UnitStep[t], {t, 0, 6}]}],
  {candidate, candidates}
];

plot = Plot[
  Evaluate[responses[[All, 3]]], {t, 0, 6},
  PlotRange -> All,
  GridLines -> Automatic,
  AxesLabel -> {\"Time (s)\", \"y_m\"},
  PlotLegends -> (\"zeta=\" <> ToString[NumberForm[#[[1]], {2, 2}]] <>
      \", omega_n=\" <> ToString[#[[2]]] <> \" rad/s\" & /@ candidates),
  PlotLabel -> \"Second-order reference-model trade-offs\"
];

Export[\"Chapter6_Lesson2.png\", plot, ImageResolution -> 180];

sampleTimes = Range[0, 6, 0.001];

sampled = Table[
  With[{zeta = response[[1]], omegaN = response[[2]], y = response[[3]]},
    Module[
      {values, peak, overshoot, t10, t90, settlingCandidates,
       settling, acceleration},

      values = y /. t -> # & /@ sampleTimes;
      peak = Max[values];
      overshoot = Max[0, 100 (peak - 1)];

      t10 = FirstCase[
        Transpose[{sampleTimes, values}],
        {time_, value_} /; value >= 0.1 :> time,
        Missing[\"NotFound\"]
      ];

      t90 = FirstCase[
        Transpose[{sampleTimes, values}],
        {time_, value_} /; value >= 0.9 :> time,
        Missing[\"NotFound\"]
      ];

      settlingCandidates =
        Select[Range[Length[values]], Abs[values[[#]] - 1] > 0.02 &];

      settling = If[
        settlingCandidates === {}, 0,
        With[{last = Last[settlingCandidates]},
          If[last < Length[sampleTimes],
            sampleTimes[[last + 1]], Indeterminate]]
      ];

      acceleration = Table[
        -2 zeta omegaN D[y, t] - omegaN^2 y + omegaN^2 /. t -> time,
        {time, sampleTimes}
      ];

      {zeta, omegaN, overshoot, t90 - t10, settling,
       Max[Abs[acceleration]]}
    ]
  ],
  {response, responses}
];

Grid[
  Prepend[sampled, {\"zeta\", \"omega_n\", \"Overshoot (%)\",
    \"Rise 10-90 (s)\", \"Settling 2% (s)\", \"Peak acceleration\"}],
  Frame -> All
]

Export[
  \"Chapter6_Lesson2_metrics.csv\",
  Prepend[sampled, {\"zeta\", \"omega_n\", \"overshoot_percent\",
    \"rise_time_10_90\", \"settling_time_2_percent\", \"peak_acceleration\"}]
];

maximumOvershoot = 0.05;
maximumSettlingTime = 1.5;

zetaMinimum =
  -Log[maximumOvershoot]/Sqrt[Pi^2 + Log[maximumOvershoot]^2];

omegaNMinimum = 4/(zetaMinimum maximumSettlingTime);

N[{zetaMinimum, omegaNMinimum}]
"
  ], "Input"]
},
WindowTitle -> "Chapter6_Lesson2"
]        

Comparing \( (\zeta,\omega_n)=(0.5,4) \) and \( (0.7,4) \) isolates the effect of damping: the lower-damped model rises rapidly but has substantially greater overshoot. Comparing \( (0.7,4) \) and \( (0.7,7) \) isolates time-scale compression: the faster model has essentially the same normalized overshoot, but markedly shorter rise and settling times and a much larger initial acceleration.

The critically damped candidate \( (1,4) \) is monotonic, but its early rise differs from the underdamped candidates. These comparisons show why a single scalar specification such as settling time is insufficient for reference model selection.

Problem 1 (Model Construction): Construct a unity-DC-gain second-order reference model with poles at \( -3\pm4j \). Determine \( \omega_n \), \( \zeta \), and \( G_m(s) \).

Solution: The pole magnitude and real part give

\[ \omega_n=\sqrt{3^2+4^2}=5, \qquad \zeta\omega_n=3, \qquad \zeta=\frac{3}{5}=0.6. \]

Therefore

\[ G_m(s)=\frac{25}{s^2+6s+25}. \]

Problem 2 (Overshoot Constraint): Determine the minimum damping ratio for a maximum step overshoot of \( 10\% \).

Solution:

\[ \zeta_{\min}= \frac{-\ln(0.10)}{\sqrt{\pi^2+[\ln(0.10)]^2}} \approx0.5912. \]

Any underdamped candidate with \( \zeta\geq0.5912 \) meets the overshoot limit.

Problem 3 (Feasible Natural-Frequency Interval): A design requires \( \%OS\leq5\% \) and \( t_s\leq2\;\text{s} \). Hardware analysis imposes \( \omega_n\leq5\;\text{rad/s} \). Find an admissible interval using the standard settling approximation.

Solution: The overshoot requirement gives \( \zeta\geq0.6901 \). Selecting the smallest admissible damping produces the most demanding lower bound:

\[ \omega_n\geq\frac{4}{(0.6901)(2)}\approx2.898. \]

Hence one feasible interval is

\[ 2.898\leq\omega_n\leq5\;\text{rad/s}, \qquad \zeta\geq0.6901. \]

A practical choice might be \( \zeta=0.7 \) and \( \omega_n=3.2\;\text{rad/s} \), followed by actuator and bandwidth simulation.

Problem 4 (Acceleration Scaling Proof): Show that increasing \( \omega_n \) from \( \omega_1 \) to \( \omega_2 \) changes the initial unit-step acceleration by the factor \( (\omega_2/\omega_1)^2 \).

Solution: At zero initial state, the model equation gives

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

Therefore

\[ \frac{\ddot y_{m,2}(0^+)}{\ddot y_{m,1}(0^+)} =\frac{\omega_2^2}{\omega_1^2}= \left(\frac{\omega_2}{\omega_1}\right)^2. \]

Thus a 50% increase in natural frequency raises initial desired acceleration by a factor of \( 1.5^2=2.25 \).

Problem 5 (Bandwidth Derivation): Starting from \( |G_m(j\omega_{bw})|^2=1/2 \), derive the bandwidth formula.

Solution: Let \( q=(\omega_{bw}/\omega_n)^2 \). Dividing numerator and denominator by \( \omega_n^4 \) gives

\[ \frac{1}{(1-q)^2+4\zeta^2q}=\frac{1}{2}. \]

Rearrangement yields

\[ q^2+(4\zeta^2-2)q-1=0. \]

The physically meaningful positive root is

\[ q=1-2\zeta^2+\sqrt{2-4\zeta^2+4\zeta^4}, \]

and therefore

\[ \omega_{bw}=\omega_n\sqrt{q}. \]

Problem 6 (Design Judgment): Two candidates meet the same overshoot requirement: Model A has \( \zeta=0.7,\omega_n=4 \), while Model B has \( \zeta=0.7,\omega_n=8 \). The plant has a lightly damped neglected mode near \( 10\;\text{rad/s} \). Which model is safer?

Solution: Model A is safer. Model B places desired closed-loop activity much closer to the neglected mode and asks for four times the initial acceleration of Model A. Meeting nominal step specifications does not compensate for inadequate separation from unmodelled dynamics. The exact safe model must be confirmed using the known actuator and uncertainty information.

A second-order MRAC reference model is normally parameterized by damping ratio and natural frequency. Damping primarily determines normalized response shape and overshoot; natural frequency determines time scale and bandwidth. Time-domain specifications generate lower bounds on these parameters, while actuator, bandwidth, noise, sampling, and unmodelled-dynamics considerations generate upper limits on aggressiveness. The selected model must be stable, physically feasible, and fixed before the error dynamics and Lyapunov-based update laws are derived in the following lessons.