After completing this lesson, students should be able to:

• express a first-order scalar plant in coefficient form and in gain/time-constant form;
• derive the exact state response for arbitrary and constant inputs;
• distinguish pole uncertainty, steady-state-gain uncertainty, and unknown control direction;
• prove exponential stability of the unforced plant when \( a > 0 \);
• derive a bounded-input bound and interpret its dependence on \( a \) and \( b \);
• calculate the fixed controller gains that would be used if the plant parameters were known; and
• reproduce the numerical experiment in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

The lesson assumes familiarity with first-order differential equations, transfer functions, pole locations, state feedback, and Lyapunov stability. The concepts of parametric uncertainty, linearly parameterized models, matching conditions, ideal parameters, and the general MRAC architecture were introduced in earlier chapters.

Consider the measured scalar state/output \( x(t)\in\mathbb{R} \) and scalar control input \( u(t)\in\mathbb{R} \). The uncertain plant is

\[ \dot{x}(t)=-a\,x(t)+b\,u(t),\qquad x(0)=x_0, \]

where \( a \) and \( b \) are unknown constants. The coefficient \( a \) determines the open-loop pole, while \( b \) determines the magnitude and direction through which the input acts on the state.

flowchart LR
  U["Input u(t)"] --> P["Unknown first-order plant"]
  TH["Unknown constants: a and b"] --> P
  P --> X["Measured state/output x(t)"]
  P --> D["Dynamics: x_dot = -a*x + b*u"]
        

A common introductory MRAC assumption is that the parameters belong to a fixed uncertainty set:

\[ a\in\mathcal{A},\qquad b\in\mathcal{B},\qquad 0<b_{\min}\leq |b|\leq b_{\max}. \]

The last condition excludes \( b=0 \), because a plant with \( b=0 \) cannot be influenced by the input. In the basic MRAC development, the sign of \( b \) is normally assumed known even though its magnitude is unknown. The more difficult unknown-sign case is postponed to Chapter 16.

For a stable first-order plant, assume \( a > 0 \). Dividing the differential equation by \( a \) gives

\[ \frac{1}{a}\dot{x}(t)+x(t)=\frac{b}{a}u(t). \]

Define the physical parameters

\[ \tau=\frac{1}{a},\qquad K=\frac{b}{a}. \]

The plant can then be written as

\[ \tau\dot{x}(t)+x(t)=K\,u(t), \qquad \tau > 0. \]

Conversely, the coefficient parameters are recovered from \( K \) and \( \tau \) as

\[ a=\frac{1}{\tau},\qquad b=\frac{K}{\tau}. \]

The two parameterizations contain the same information, but they emphasize different control properties:

• \( (a,b) \) makes the differential equation linear in the unknown coefficients and is convenient for adaptive-control derivations.
• \( (K,\tau) \) directly exposes steady-state gain and transient speed and is convenient for physical interpretation.

With zero initial condition, the unilateral Laplace transform yields

\[ sX(s)=-aX(s)+bU(s), \]

\[ G(s)=\frac{X(s)}{U(s)} =\frac{b}{s+a} =\frac{K}{\tau s+1}. \]

Therefore, the plant has one pole at \( s=-a \), no finite zero, and relative degree one.

Rewrite the plant equation in standard linear differential-equation form:

\[ \dot{x}(t)+a\,x(t)=b\,u(t). \]

Multiplying both sides by the integrating factor \( e^{at} \) gives

\[ e^{at}\dot{x}(t)+a e^{at}x(t) =b e^{at}u(t). \]

The left-hand side is an exact derivative:

\[ \frac{d}{dt}\left(e^{at}x(t)\right)=b e^{at}u(t). \]

Integrating from \( 0 \) to \( t \),

\[ e^{at}x(t)-x_0 =b\int_0^t e^{a\sigma}u(\sigma)\,d\sigma. \]

Hence the exact solution is

\[ x(t)=e^{-at}x_0+ b\int_0^t e^{-a(t-\sigma)}u(\sigma)\,d\sigma. \]

This decomposition separates the zero-input response \( e^{-at}x_0 \) from the zero-state convolution response. It also shows exactly where the unknown parameters enter: \( a \) changes the memory kernel, while \( b \) scales and reverses the forced response.

For a constant input \( u(t)=u_0 \) and \( a\neq 0 \), the exact response becomes

\[ x(t)=e^{-at}x_0+ b u_0\int_0^t e^{-a(t-\sigma)}\,d\sigma. \]

Evaluating the integral gives

\[ x(t)=e^{-at}x_0+\frac{b}{a}u_0\left(1-e^{-at}\right). \]

Equivalently,

\[ x(t)=x_{\mathrm{ss}}+ \left(x_0-x_{\mathrm{ss}}\right)e^{-at}, \qquad x_{\mathrm{ss}}=\frac{b}{a}u_0. \]

If \( a > 0 \), then \( e^{-at} → 0 \) as \( t → \infty \), and therefore

\[ \lim_{t→\infty}x(t)=\frac{b}{a}u_0=K u_0. \]

At one time constant \( t=\tau=1/a \),

\[ \frac{x(\tau)-x_0}{x_{\mathrm{ss}}-x_0} =1-e^{-1}\approx 0.6321. \]

Thus a stable first-order step response has reached approximately 63.2% of its total change after one time constant. This provides a direct physical interpretation of uncertainty in \( a \).

The special case \( a=0 \) is not a first-order lag but an integrator:

\[ \dot{x}=b u,\qquad x(t)=x_0+b\int_0^t u(\sigma)\,d\sigma. \]

Under a nonzero constant input, its state grows linearly and no finite steady state exists.

6.1 Unforced exponential stability

Set \( u(t)=0 \). The plant becomes

\[ \dot{x}=-a x. \]

Choose the quadratic Lyapunov function

\[ V(x)=\frac{1}{2}x^2. \]

Its derivative along plant trajectories is

\[ \dot{V}=x\dot{x}=-a x^2=-2aV. \]

For \( a > 0 \), \( \dot{V}=-2aV < 0 \) whenever \( x\neq 0 \). Solving the scalar differential inequality exactly gives

\[ V(t)=e^{-2at}V(0), \qquad |x(t)|=e^{-at}|x_0|. \]

Therefore, the origin is globally exponentially stable with decay rate \( a \). If \( a<0 \), the zero-input mode grows exponentially, and the open-loop plant is unstable.

6.2 Bounded-input bound for the stable plant

Assume \( a > 0 \) and \( |u(t)|\leq \bar{u} \). Starting from the exact solution,

\[ |x(t)|\leq e^{-at}|x_0| +|b|\int_0^t e^{-a(t-\sigma)}|u(\sigma)|\,d\sigma. \]

Using \( |u(\sigma)|\leq\bar{u} \),

\[ |x(t)|\leq e^{-at}|x_0| +|b|\bar{u}\int_0^t e^{-a(t-\sigma)}\,d\sigma, \]

\[ |x(t)|\leq e^{-at}|x_0| +\frac{|b|}{a}\bar{u}\left(1-e^{-at}\right). \]

Consequently,

\[ |x(t)|\leq |x_0|+\frac{|b|}{a}\bar{u}, \qquad \limsup_{t→\infty}|x(t)| \leq \frac{|b|}{a}\bar{u}. \]

The bound reveals two effects of uncertainty. Small \( a \) produces slow decay and a large asymptotic input-to-state gain, while large \( |b| \) increases the influence of the input.

For \( a > 0 \), the time constant and DC gain are

\[ \tau(a)=\frac{1}{a}, \qquad K(a,b)=\frac{b}{a}. \]

Their local sensitivities are

\[ \frac{\partial \tau}{\partial a} =-\frac{1}{a^2}, \qquad \frac{\partial K}{\partial a} =-\frac{b}{a^2}, \qquad \frac{\partial K}{\partial b} =\frac{1}{a}. \]

The normalized sensitivities, defined by \( S_q^p=(p/q)(\partial q/\partial p) \), are

\[ S_{\tau}^{a}=-1,\qquad S_K^a=-1,\qquad S_K^b=1. \]

Hence, to first order, a 1% increase in \( a \) produces a 1% decrease in both the time constant and the DC gain, while a 1% increase in \( b \) produces a 1% increase in the DC gain. These simple relationships explain why a single fixed controller may perform very differently across plants in the uncertainty set.

For a stable plant under a constant input \( u_0 \),

\[ x_{\mathrm{ss}}(a,b)=\frac{b}{a}u_0. \]

Therefore,

\[ \frac{\partial x_{\mathrm{ss}}}{\partial a} =-\frac{b u_0}{a^2}, \qquad \frac{\partial x_{\mathrm{ss}}}{\partial b} =\frac{u_0}{a}. \]

Before constructing an adaptive controller, it is useful to determine what controller would achieve the desired dynamics if \( a \) and \( b \) were known. Let the desired first-order command dynamics be

\[ \dot{x}=-a_d x+b_d r, \qquad a_d > 0, \]

where \( r(t) \) is a bounded command. Consider the static state-feedback/feedforward structure

\[ u=k_x x+k_r r. \]

Substitution into the actual plant gives

\[ \dot{x} =-a x+b(k_x x+k_r r) =(-a+b k_x)x+b k_r r. \]

Exact coefficient matching with the desired dynamics requires

\[ -a+b k_x=-a_d, \qquad b k_r=b_d. \]

If \( b\neq 0 \), the unique matching gains are

\[ k_x^\star=\frac{a-a_d}{b}, \qquad k_r^\star=\frac{b_d}{b}. \]

These are the ideal controller parameters for this structure. They exist because the scalar plant satisfies the matching conditions. However, they cannot be calculated directly when \( a \) and \( b \) are unknown. The remainder of Chapter 5 will develop mechanisms that adjust controller parameters without first knowing the true plant constants.

The sign of \( b \) matters immediately. Since \( k_r^\star=b_d/b \), reversing the sign of \( b \) reverses the required feedforward control direction. This is why basic MRAC designs assume the input-gain sign is known.

The implementations below conduct the same constant-input experiment:

\[ x_0=0.30,\qquad u_0=2.00,\qquad 0\leq t\leq 8\ \mathrm{s},\qquad \Delta t=0.002\ \mathrm{s}. \]

Five parameter cases are simulated:

flowchart TD
  A["Choose a, b, x0, and constant input u0"] --> B["Integrate x_dot = -a*x + b*u with RK4"]
  B --> C["Compute exact response for the nominal case"]
  C --> D["Check maximum numerical error"]
  D --> E["Extract tau, DC gain, and step equilibrium"]
  E --> F["Compare stable, reversed-direction, and unstable cases"]
  F --> G["Export CSV and response plot"]
        

Each implementation uses fourth-order Runge–Kutta integration and checks the nominal numerical trajectory against the exact solution. The expected nominal equilibrium is

\[ x_{\mathrm{ss}} =\frac{2.00}{1.25}(2.00) =3.20. \]

The cross-language files deliberately use the same names, parameters, time step, and output columns so their results can be compared directly.

Libraries: Python standard library, NumPy for arrays and exact vectorized evaluation, and Matplotlib for plotting. The differential equation is integrated from scratch with RK4.

Chapter5_Lesson1.py

"""
Chapter5_Lesson1.py
First-order scalar plant with unknown gain/time constant.

Model:
    x_dot = -a*x + b*u

The script compares several admissible and non-admissible parameter cases,
checks the RK4 solution against the exact constant-input response, exports a
CSV file, and saves a plot.
"""

from __future__ import annotations

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

import matplotlib.pyplot as plt
import numpy as np


@dataclass(frozen=True)
class PlantCase:
    name: str
    a: float
    b: float


def plant_rhs(x: float, u: float, a: float, b: float) -> float:
    """Continuous-time plant right-hand side."""
    return -a * x + b * u


def rk4_step(x: float, u: float, dt: float, a: float, b: float) -> float:
    """One Runge-Kutta fourth-order step for constant u over the step."""
    k1 = plant_rhs(x, u, a, b)
    k2 = plant_rhs(x + 0.5 * dt * k1, u, a, b)
    k3 = plant_rhs(x + 0.5 * dt * k2, u, a, b)
    k4 = plant_rhs(x + dt * k3, u, a, b)
    return x + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def simulate(case: PlantCase, time: np.ndarray, x0: float, u0: float) -> np.ndarray:
    """Simulate a constant-input experiment."""
    x = np.empty_like(time, dtype=float)
    x[0] = x0
    for k in range(time.size - 1):
        dt = float(time[k + 1] - time[k])
        x[k + 1] = rk4_step(x[k], u0, dt, case.a, case.b)
    return x


def exact_constant_input(
    time: np.ndarray, x0: float, u0: float, a: float, b: float
) -> np.ndarray:
    """Exact solution for a nonzero a and constant input u0."""
    if abs(a) < 1.0e-14:
        return x0 + b * u0 * time
    x_ss = (b / a) * u0
    return x_ss + (x0 - x_ss) * np.exp(-a * time)


def describe_case(case: PlantCase, u0: float) -> str:
    """Return the main open-loop characteristics."""
    if case.a > 0.0:
        tau = 1.0 / case.a
        dc_gain = case.b / case.a
        x_ss = dc_gain * u0
        return (
            f"stable; tau={tau:.6g}; K={dc_gain:.6g}; "
            f"x_ss(step)={x_ss:.6g}"
        )
    if case.a == 0.0:
        return "integrator; no finite time constant or constant-input equilibrium"
    return "open-loop unstable; the exponential mode grows"


def main() -> None:
    x0 = 0.30
    u0 = 2.00
    dt = 0.002
    t_final = 8.0
    time = np.arange(0.0, t_final + 0.5 * dt, dt)

    cases = [
        PlantCase("nominal_stable", a=1.25, b=2.00),
        PlantCase("slow_stable", a=0.60, b=1.20),
        PlantCase("fast_stable", a=2.00, b=1.20),
        PlantCase("reversed_input_direction", a=1.25, b=-2.00),
        PlantCase("unstable", a=-0.30, b=1.20),
    ]

    trajectories: dict[str, np.ndarray] = {}
    for case in cases:
        trajectories[case.name] = simulate(case, time, x0, u0)
        print(f"{case.name}: {describe_case(case, u0)}")

    nominal = cases[0]
    exact = exact_constant_input(time, x0, u0, nominal.a, nominal.b)
    max_error = float(np.max(np.abs(trajectories[nominal.name] - exact)))
    print(f"Maximum RK4-versus-exact error for nominal case: {max_error:.3e}")
    if max_error > 1.0e-9:
        raise RuntimeError("Numerical verification failed; reduce dt or inspect the code.")

    csv_path = Path("Chapter5_Lesson1_results.csv")
    with csv_path.open("w", newline="", encoding="utf-8") as stream:
        writer = csv.writer(stream)
        writer.writerow(["time", *[case.name for case in cases], "nominal_exact"])
        for k, t_value in enumerate(time):
            writer.writerow(
                [
                    f"{t_value:.8f}",
                    *[f"{trajectories[case.name][k]:.12g}" for case in cases],
                    f"{exact[k]:.12g}",
                ]
            )
    print(f"Wrote {csv_path.resolve()}")

    for case in cases:
        plt.plot(time, trajectories[case.name], label=case.name)
    plt.xlabel("Time (s)")
    plt.ylabel("State/output x(t)")
    plt.title("First-order plant responses under parameter uncertainty")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    figure_path = Path("Chapter5_Lesson1_responses.png")
    plt.savefig(figure_path, dpi=180)
    print(f"Wrote {figure_path.resolve()}")
    plt.show()


if __name__ == "__main__":
    main()

Libraries: only the C++17 standard library. The program performs the same RK4 integration, exact-solution verification, and CSV export as the Python implementation.

Chapter5_Lesson1.cpp

/*
Chapter5_Lesson1.cpp
First-order scalar plant with unknown gain/time constant.

Build:
    g++ -std=c++17 -O2 Chapter5_Lesson1.cpp -o Chapter5_Lesson1
*/

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

struct PlantCase {
    std::string name;
    double a;
    double b;
};

double plantRhs(double x, double u, double a, double b) {
    return -a * x + b * u;
}

double rk4Step(double x, double u, double dt, double a, double b) {
    const double k1 = plantRhs(x, u, a, b);
    const double k2 = plantRhs(x + 0.5 * dt * k1, u, a, b);
    const double k3 = plantRhs(x + 0.5 * dt * k2, u, a, b);
    const double k4 = plantRhs(x + dt * k3, u, a, b);
    return x + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
}

std::vector<double> simulate(
    const PlantCase& plant,
    const std::vector<double>& time,
    double x0,
    double u0
) {
    std::vector<double> x(time.size(), 0.0);
    x.front() = x0;
    for (std::size_t k = 0; k + 1 < time.size(); ++k) {
        const double dt = time[k + 1] - time[k];
        x[k + 1] = rk4Step(x[k], u0, dt, plant.a, plant.b);
    }
    return x;
}

std::vector<double> exactConstantInput(
    const std::vector<double>& time,
    double x0,
    double u0,
    double a,
    double b
) {
    std::vector<double> exact(time.size(), 0.0);
    if (std::abs(a) < 1.0e-14) {
        for (std::size_t k = 0; k < time.size(); ++k) {
            exact[k] = x0 + b * u0 * time[k];
        }
        return exact;
    }

    const double xSs = (b / a) * u0;
    for (std::size_t k = 0; k < time.size(); ++k) {
        exact[k] = xSs + (x0 - xSs) * std::exp(-a * time[k]);
    }
    return exact;
}

void printDescription(const PlantCase& plant, double u0) {
    std::cout << plant.name << ": ";
    if (plant.a > 0.0) {
        const double tau = 1.0 / plant.a;
        const double gain = plant.b / plant.a;
        std::cout << "stable; tau=" << tau
                  << "; K=" << gain
                  << "; x_ss(step)=" << gain * u0 << '\n';
    } else if (plant.a == 0.0) {
        std::cout << "integrator; no finite constant-input equilibrium\n";
    } else {
        std::cout << "open-loop unstable; the exponential mode grows\n";
    }
}

int main() {
    try {
        const double x0 = 0.30;
        const double u0 = 2.00;
        const double dt = 0.002;
        const double tFinal = 8.0;
        const std::size_t sampleCount =
            static_cast<std::size_t>(std::llround(tFinal / dt)) + 1U;

        std::vector<double> time(sampleCount, 0.0);
        for (std::size_t k = 0; k < sampleCount; ++k) {
            time[k] = static_cast<double>(k) * dt;
        }

        const std::vector<PlantCase> cases = {
            {"nominal_stable", 1.25, 2.00},
            {"slow_stable", 0.60, 1.20},
            {"fast_stable", 2.00, 1.20},
            {"reversed_input_direction", 1.25, -2.00},
            {"unstable", -0.30, 1.20}
        };

        std::vector<std::vector<double>> trajectories;
        trajectories.reserve(cases.size());
        for (const PlantCase& plant : cases) {
            trajectories.push_back(simulate(plant, time, x0, u0));
            printDescription(plant, u0);
        }

        const std::vector<double> exact =
            exactConstantInput(time, x0, u0, cases.front().a, cases.front().b);
        double maxError = 0.0;
        for (std::size_t k = 0; k < time.size(); ++k) {
            maxError = std::max(
                maxError,
                std::abs(trajectories.front()[k] - exact[k])
            );
        }
        std::cout << std::scientific
                  << "Maximum RK4-versus-exact error: " << maxError << '\n';
        if (maxError > 1.0e-9) {
            throw std::runtime_error("Numerical verification failed.");
        }

        std::ofstream csv("Chapter5_Lesson1_results.csv");
        if (!csv) {
            throw std::runtime_error("Cannot open output CSV file.");
        }
        csv << "time";
        for (const PlantCase& plant : cases) {
            csv << ',' << plant.name;
        }
        csv << ",nominal_exact\n";
        csv << std::setprecision(12);

        for (std::size_t k = 0; k < time.size(); ++k) {
            csv << time[k];
            for (const auto& trajectory : trajectories) {
                csv << ',' << trajectory[k];
            }
            csv << ',' << exact[k] << '\n';
        }

        std::cout << "Wrote Chapter5_Lesson1_results.csv\n";
        return 0;
    } catch (const std::exception& error) {
        std::cerr << "Error: " << error.what() << '\n';
        return 1;
    }
}

Libraries: Java standard library only. A Java record stores each plant case, and the results are exported with UTF-8 file I/O.

Chapter5_Lesson1.java

/*
Chapter5_Lesson1.java
First-order scalar plant with unknown gain/time constant.

Build and run:
    javac Chapter5_Lesson1.java
    java Chapter5_Lesson1
*/

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

public final class Chapter5_Lesson1 {
    private record PlantCase(String name, double a, double b) {}

    private static double plantRhs(double x, double u, double a, double b) {
        return -a * x + b * u;
    }

    private static double rk4Step(double x, double u, double dt, double a, double b) {
        double k1 = plantRhs(x, u, a, b);
        double k2 = plantRhs(x + 0.5 * dt * k1, u, a, b);
        double k3 = plantRhs(x + 0.5 * dt * k2, u, a, b);
        double k4 = plantRhs(x + dt * k3, u, a, b);
        return x + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
    }

    private static double[] simulate(
        PlantCase plant,
        double[] time,
        double x0,
        double u0
    ) {
        double[] x = new double[time.length];
        x[0] = x0;
        for (int k = 0; k < time.length - 1; k++) {
            double dt = time[k + 1] - time[k];
            x[k + 1] = rk4Step(x[k], u0, dt, plant.a(), plant.b());
        }
        return x;
    }

    private static double[] exactConstantInput(
        double[] time,
        double x0,
        double u0,
        double a,
        double b
    ) {
        double[] exact = new double[time.length];
        if (Math.abs(a) < 1.0e-14) {
            for (int k = 0; k < time.length; k++) {
                exact[k] = x0 + b * u0 * time[k];
            }
            return exact;
        }

        double xSs = (b / a) * u0;
        for (int k = 0; k < time.length; k++) {
            exact[k] = xSs + (x0 - xSs) * Math.exp(-a * time[k]);
        }
        return exact;
    }

    private static void printDescription(PlantCase plant, double u0) {
        System.out.print(plant.name() + ": ");
        if (plant.a() > 0.0) {
            double tau = 1.0 / plant.a();
            double gain = plant.b() / plant.a();
            System.out.printf(
                Locale.US,
                "stable; tau=%.6g; K=%.6g; x_ss(step)=%.6g%n",
                tau,
                gain,
                gain * u0
            );
        } else if (plant.a() == 0.0) {
            System.out.println("integrator; no finite constant-input equilibrium");
        } else {
            System.out.println("open-loop unstable; the exponential mode grows");
        }
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);

        final double x0 = 0.30;
        final double u0 = 2.00;
        final double dt = 0.002;
        final double tFinal = 8.0;
        final int sampleCount = (int) Math.round(tFinal / dt) + 1;

        double[] time = new double[sampleCount];
        for (int k = 0; k < sampleCount; k++) {
            time[k] = k * dt;
        }

        List<PlantCase> cases = List.of(
            new PlantCase("nominal_stable", 1.25, 2.00),
            new PlantCase("slow_stable", 0.60, 1.20),
            new PlantCase("fast_stable", 2.00, 1.20),
            new PlantCase("reversed_input_direction", 1.25, -2.00),
            new PlantCase("unstable", -0.30, 1.20)
        );

        List<double[]> trajectories = new ArrayList<>();
        for (PlantCase plant : cases) {
            trajectories.add(simulate(plant, time, x0, u0));
            printDescription(plant, u0);
        }

        PlantCase nominal = cases.get(0);
        double[] exact = exactConstantInput(
            time,
            x0,
            u0,
            nominal.a(),
            nominal.b()
        );

        double maxError = 0.0;
        for (int k = 0; k < sampleCount; k++) {
            maxError = Math.max(maxError, Math.abs(trajectories.get(0)[k] - exact[k]));
        }
        System.out.printf(
            Locale.US,
            "Maximum RK4-versus-exact error: %.3e%n",
            maxError
        );
        if (maxError > 1.0e-9) {
            throw new IllegalStateException("Numerical verification failed.");
        }

        Path output = Path.of("Chapter5_Lesson1_results.csv");
        try (BufferedWriter writer = Files.newBufferedWriter(
            output,
            StandardCharsets.UTF_8
        )) {
            writer.write("time");
            for (PlantCase plant : cases) {
                writer.write("," + plant.name());
            }
            writer.write(",nominal_exact\n");

            for (int k = 0; k < sampleCount; k++) {
                writer.write(String.format(Locale.US, "%.8f", time[k]));
                for (double[] trajectory : trajectories) {
                    writer.write(String.format(Locale.US, ",%.12g", trajectory[k]));
                }
                writer.write(String.format(Locale.US, ",%.12g%n", exact[k]));
            }
        } catch (IOException error) {
            throw new RuntimeException("Failed to write CSV output.", error);
        }

        System.out.println("Wrote " + output.toAbsolutePath());
    }
}

Libraries: base MATLAB is used for numerical integration, tables, plotting, and export. If Simulink is installed, the same script programmatically creates Chapter5_Lesson1_Simulink.slx with the block equation x_dot = b*u - a*x.

Chapter5_Lesson1.m

% Chapter5_Lesson1.m
% First-order scalar plant with unknown gain/time constant.
%
% Model:
%   x_dot = -a*x + b*u
%
% The script verifies RK4 against the exact constant-input response, exports
% results, plots uncertainty cases, and optionally creates a Simulink model.

clear; clc; close all;

x0 = 0.30;
u0 = 2.00;
dt = 0.002;
tFinal = 8.0;
t = (0:dt:tFinal)';

cases = struct( ...
    'name', {'nominal_stable', 'slow_stable', 'fast_stable', ...
             'reversed_input_direction', 'unstable'}, ...
    'a', {1.25, 0.60, 2.00, 1.25, -0.30}, ...
    'b', {2.00, 1.20, 1.20, -2.00, 1.20});

x = zeros(numel(t), numel(cases));
x(1, :) = x0;

for j = 1:numel(cases)
    for k = 1:numel(t)-1
        x(k+1, j) = rk4Step( ...
            x(k, j), u0, dt, cases(j).a, cases(j).b);
    end

    if cases(j).a > 0
        tau = 1 / cases(j).a;
        K = cases(j).b / cases(j).a;
        fprintf('%s: stable; tau=%.6g; K=%.6g; x_ss(step)=%.6g\n', ...
            cases(j).name, tau, K, K*u0);
    elseif cases(j).a == 0
        fprintf('%s: integrator; no finite constant-input equilibrium\n', ...
            cases(j).name);
    else
        fprintf('%s: open-loop unstable; the exponential mode grows\n', ...
            cases(j).name);
    end
end

aNom = cases(1).a;
bNom = cases(1).b;
xSsNom = (bNom/aNom)*u0;
xExact = xSsNom + (x0 - xSsNom).*exp(-aNom*t);
maxError = max(abs(x(:, 1) - xExact));
fprintf('Maximum RK4-versus-exact error: %.3e\n', maxError);
assert(maxError <= 1.0e-9, ...
    'Numerical verification failed; reduce dt or inspect the code.');

resultTable = table(t, 'VariableNames', {'time'});
for j = 1:numel(cases)
    resultTable.(cases(j).name) = x(:, j);
end
resultTable.nominal_exact = xExact;
writetable(resultTable, 'Chapter5_Lesson1_results.csv');

figure;
plot(t, x, 'LineWidth', 1.2);
grid on;
xlabel('Time (s)');
ylabel('State/output x(t)');
title('First-order plant responses under parameter uncertainty');
legend({cases.name}, 'Interpreter', 'none', 'Location', 'best');
exportgraphics(gcf, 'Chapter5_Lesson1_responses.png', 'Resolution', 180);

% Optional Simulink construction. The generated model implements
% x_dot = b*u - a*x for the nominal parameter values.
if license('test', 'Simulink') && exist('new_system', 'file') == 2
    model = 'Chapter5_Lesson1_Simulink';

    if bdIsLoaded(model)
        close_system(model, 0);
    end
    if isfile([model '.slx'])
        delete([model '.slx']);
    end

    new_system(model);
    open_system(model);

    add_block('simulink/Sources/Step', [model '/Step']);
    set_param([model '/Step'], ...
        'Time', '0', ...
        'Before', '0', ...
        'After', num2str(u0));

    add_block('simulink/Math Operations/Gain', [model '/Gain_b']);
    set_param([model '/Gain_b'], 'Gain', num2str(bNom));

    add_block('simulink/Math Operations/Gain', [model '/Gain_a']);
    set_param([model '/Gain_a'], 'Gain', num2str(aNom));

    add_block('simulink/Math Operations/Sum', [model '/Sum']);
    set_param([model '/Sum'], 'Inputs', '+-');

    add_block('simulink/Continuous/Integrator', [model '/Plant_state']);
    set_param([model '/Plant_state'], 'InitialCondition', num2str(x0));

    add_block('simulink/Sinks/Scope', [model '/Scope']);

    add_line(model, 'Step/1', 'Gain_b/1');
    add_line(model, 'Gain_b/1', 'Sum/1');
    add_line(model, 'Sum/1', 'Plant_state/1');
    add_line(model, 'Plant_state/1', 'Gain_a/1');
    add_line(model, 'Gain_a/1', 'Sum/2');
    add_line(model, 'Plant_state/1', 'Scope/1');

    set_param([model '/Step'], 'Position', [35 70 65 100]);
    set_param([model '/Gain_b'], 'Position', [105 65 165 105]);
    set_param([model '/Sum'], 'Position', [225 70 250 100]);
    set_param([model '/Plant_state'], 'Position', [305 65 345 105]);
    set_param([model '/Gain_a'], 'Position', [305 150 365 190]);
    set_param([model '/Scope'], 'Position', [430 65 460 105]);

    set_param(model, 'StopTime', num2str(tFinal));
    save_system(model);
    fprintf('Created %s.slx\n', model);
end

function xNext = rk4Step(x, u, dt, a, b)
    f = @(state) -a*state + b*u;
    k1 = f(x);
    k2 = f(x + 0.5*dt*k1);
    k3 = f(x + 0.5*dt*k2);
    k4 = f(x + dt*k3);
    xNext = x + (dt/6)*(k1 + 2*k2 + 2*k3 + k4);
end

Libraries: Wolfram Language built-ins only. The notebook uses associations for parameter cases, a from-scratch RK4 function, symbolic-form exact evaluation, CSV export, and ListLinePlot.

Chapter5_Lesson1.nb

Notebook[{
  Cell["Chapter 5, Lesson 1: First-Order Scalar Plant with Unknown Gain/Time Constant", "Title"],
  Cell["This notebook studies x_dot = -a x + b u, compares uncertainty cases, verifies a fourth-order Runge-Kutta implementation against the exact constant-input solution, and exports the results.", "Text"],
  Cell[
"ClearAll[plantRhs, rk4Step, simulate, exactConstantInput, describeCase];
plantRhs[x_, u_, a_, b_] := -a*x + b*u;
rk4Step[x_, u_, dt_, a_, b_] := Module[{k1, k2, k3, k4},
  k1 = plantRhs[x, u, a, b]; k2 = plantRhs[x + (dt/2)*k1, u, a, b];
  k3 = plantRhs[x + (dt/2)*k2, u, a, b]; k4 = plantRhs[x + dt*k3, u, a, b];
  x + (dt/6)*(k1 + 2*k2 + 2*k3 + k4)
];
simulate[case_Association, time_List, x0_, u0_] := FoldList[
  rk4Step[#1, u0, #2[[2]] - #2[[1]], case[\"a\"], case[\"b\"]] &,
  x0, Partition[time, 2, 1]
];
exactConstantInput[time_List, x0_, u0_, a_, b_] := If[
  Abs[a] < 10^-14, x0 + b*u0*time,
  With[{xSs = (b/a)*u0}, xSs + (x0 - xSs)*Exp[-a*time]]
];
x0 = 0.30; u0 = 2.00; dt = 0.002; tFinal = 8.0;
time = N@Range[0, tFinal, dt];
cases = {
  <|\"name\" -> \"nominal_stable\", \"a\" -> 1.25, \"b\" -> 2.00|>,
  <|\"name\" -> \"slow_stable\", \"a\" -> 0.60, \"b\" -> 1.20|>,
  <|\"name\" -> \"fast_stable\", \"a\" -> 2.00, \"b\" -> 1.20|>,
  <|\"name\" -> \"reversed_input_direction\", \"a\" -> 1.25, \"b\" -> -2.00|>,
  <|\"name\" -> \"unstable\", \"a\" -> -0.30, \"b\" -> 1.20|>
};
trajectories = Association@Table[
  case[\"name\"] -> simulate[case, time, x0, u0], {case, cases}
];
Do[
  If[case[\"a\"] > 0,
    Print[case[\"name\"], \": stable; tau=\", 1/case[\"a\"], \"; K=\", case[\"b\"]/case[\"a\"],
      \"; x_ss(step)=\", (case[\"b\"]/case[\"a\"])*u0],
    If[case[\"a\"] == 0, Print[case[\"name\"], \": integrator; no finite equilibrium\"],
      Print[case[\"name\"], \": open-loop unstable\"]]
  ], {case, cases}
];
nominal = First[cases];
exact = exactConstantInput[time, x0, u0, nominal[\"a\"], nominal[\"b\"]];
maxError = Max@Abs[trajectories[\"nominal_stable\"] - exact];
Print[\"Maximum RK4-versus-exact error: \", ScientificForm[maxError, 4]];
If[maxError > 10^-9, Print[\"Warning: numerical verification tolerance exceeded.\"]];
headers = Join[{\"time\"}, cases[[All, \"name\"]], {\"nominal_exact\"}];
data = Transpose@Join[{time}, (trajectories[#] & /@ cases[[All, \"name\"]]), {exact}];
Export[\"Chapter5_Lesson1_results.csv\", Prepend[data, headers]];
responsePlot = ListLinePlot[
  Table[Transpose[{time, trajectories[case[\"name\"]]}], {case, cases}],
  PlotLegends -> cases[[All, \"name\"]], Frame -> True,
  FrameLabel -> {\"Time (s)\", \"State/output x(t)\"},
  PlotLabel -> \"First-order plant responses under parameter uncertainty\", ImageSize -> Large
];
Export[\"Chapter5_Lesson1_responses.png\", responsePlot, ImageResolution -> 180];
responsePlot", "Input"]
}, WindowTitle -> "Chapter5_Lesson1.nb", StyleDefinitions -> "Default.nb"]        

For the nominal case \( a=1.25 \) and \( b=2.00 \),

\[ \tau=\frac{1}{1.25}=0.8\ \mathrm{s}, \qquad K=\frac{2.00}{1.25}=1.6, \qquad x_{\mathrm{ss}}=1.6(2.00)=3.2. \]

At \( t=\tau=0.8\ \mathrm{s} \),

\[ x(0.8)=3.2+(0.3-3.2)e^{-1}\approx 2.1331. \]

The slow stable case has \( \tau=1/0.6\approx1.6667\ \mathrm{s} \) and \( K=1.2/0.6=2 \), so it converges more slowly but to a larger equilibrium \( x_{\mathrm{ss}}=4 \). The fast stable case has \( \tau=0.5\ \mathrm{s} \) and \( K=0.6 \), so it converges faster but only to \( x_{\mathrm{ss}}=1.2 \).

The reversed-input-direction case retains the same stable pole as the nominal plant but has \( K=-1.6 \). It therefore moves toward \( -3.2 \) under the same positive input. The unstable case has pole \( s=0.30 \), because \( -a=0.30 \), and its homogeneous response grows as \( e^{0.30t} \).

The verification tolerance is intentionally strict. With \( \Delta t=0.002 \), the tested Python, C++, and Java implementations produce a nominal RK4-versus-exact maximum error of approximately \( 3.45\times10^{-13} \).

Problem 1: Conversion Between Parameterizations

A stable first-order plant is described by \( \tau\dot{x}+x=K u \), where \( \tau=0.4 \) and \( K=-1.5 \). Determine \( a \) and \( b \) in \( \dot{x}=-a x+b u \). State the pole and the input direction.

Solution:

\[ a=\frac{1}{\tau}=\frac{1}{0.4}=2.5, \qquad b=\frac{K}{\tau}=\frac{-1.5}{0.4}=-3.75. \]

The pole is \( s=-a=-2.5 \), so the plant is exponentially stable. Since \( b<0 \), a positive input initially drives the state in the negative direction when the state contribution is absent.

Problem 2: Exact Response and Time Constant

Let \( a=2 \), \( b=3 \), \( x_0=-1 \), and \( u(t)=2 \). Derive \( x(t) \), calculate the equilibrium, and find the state at one time constant.

Solution:

\[ x_{\mathrm{ss}}=\frac{b}{a}u_0 =\frac{3}{2}(2)=3, \qquad \tau=\frac{1}{a}=0.5. \]

\[ x(t)=3+(-1-3)e^{-2t} =3-4e^{-2t}. \]

At one time constant,

\[ x(0.5)=3-4e^{-1}\approx1.5285. \]

The total state change is \( 3-(-1)=4 \), and the achieved change after one time constant is \( 1.5285-(-1)=2.5285\approx0.6321(4) \).

Problem 3: Lyapunov and Bounded-Input Proof

For \( \dot{x}=-a x+b u \) with \( a\geq a_{\min}>0 \), \( |b|\leq b_{\max} \), and \( |u(t)|\leq\bar{u} \), derive a uniform bound that does not depend on the exact values of \( a \) and \( b \).

Solution:

From the exact-response inequality,

\[ |x(t)|\leq e^{-at}|x_0| +\frac{|b|}{a}\bar{u}\left(1-e^{-at}\right). \]

Since \( a\geq a_{\min} \) and \( |b|\leq b_{\max} \),

\[ e^{-at}\leq e^{-a_{\min}t}, \qquad \frac{|b|}{a}\leq\frac{b_{\max}}{a_{\min}}. \]

Therefore,

\[ |x(t)|\leq e^{-a_{\min}t}|x_0| +\frac{b_{\max}}{a_{\min}}\bar{u} \left(1-e^{-a_{\min}t}\right). \]

In particular,

\[ \limsup_{t→\infty}|x(t)| \leq\frac{b_{\max}}{a_{\min}}\bar{u}. \]

Problem 4: Known-Parameter Matching Gains

The plant is \( \dot{x}=-0.8x+1.6u \). The desired command dynamics are \( \dot{x}=-2x+3r \). For \( u=k_x x+k_r r \), calculate the exact matching gains.

Solution:

\[ k_x^\star=\frac{a-a_d}{b} =\frac{0.8-2}{1.6} =-0.75, \qquad k_r^\star=\frac{b_d}{b} =\frac{3}{1.6} =1.875. \]

Verification by substitution gives

\[ \dot{x} =-0.8x+1.6(-0.75x+1.875r) =-2x+3r. \]

Problem 5: Consequences of Using Nominal Gains on the Wrong Plant

A controller is designed for nominal parameters \( a_n=1 \) and \( b_n=2 \) to match \( \dot{x}=-3x+2r \). The resulting nominal gains are then applied to the true plant with \( a=0.5 \) and \( b=1 \). Determine the actual closed-loop dynamics.

Solution:

The nominal matching gains are

\[ k_{x,n}=\frac{a_n-a_d}{b_n} =\frac{1-3}{2}=-1, \qquad k_{r,n}=\frac{b_d}{b_n} =\frac{2}{2}=1. \]

Applying them to the true plant gives

\[ \dot{x} =-0.5x+1(-x+r) =-1.5x+r. \]

The actual pole is \( -1.5 \), not \( -3 \), and the command coefficient is \( 1 \), not \( 2 \). This mismatch is the central motivation for adapting the controller parameters rather than fixing them at nominal values.

The uncertain scalar plant \( \dot{x}=-a x+b u \) is the simplest setting in which the essential MRAC issues already appear. For \( a>0 \), its time constant and DC gain are \( \tau=1/a \) and \( K=b/a \). Its exact response separates the decaying initial-condition mode from the input convolution. A quadratic Lyapunov function proves global exponential stability of the unforced stable plant, and the exact solution gives a transparent bounded-input bound.

If the plant parameters were known, the controller \( u=k_xx+k_rr \) could match desired first-order dynamics with \( k_x^\star=(a-a_d)/b \) and \( k_r^\star=b_d/b \). Because these ideal gains depend on unknown plant constants, subsequent lessons will construct the reference model, derive the tracking-error dynamics, and then introduce adaptive update laws.