In regulation, the control objective is usually to make the state approach the origin or to keep the output near a fixed operating point. In tracking, the desired signal is explicitly present. Therefore, the plant model must specify not only the state and control input, but also where the reference enters, where the disturbance enters, and which output must follow the reference. A convenient continuous-time LTI tracking model is

\[ \begin{aligned} \dot{x}(t) &= A x(t) + B u(t) + E d(t), \\ y(t) &= C x(t) + D u(t) + F d(t), \\ e(t) &= y(t) - r(t). \end{aligned} \]

Here \( x(t)\in\mathbb{R}^n \) is the state, \( u(t)\in\mathbb{R}^m \) is the manipulated input, \( d(t)\in\mathbb{R}^q \) is an external disturbance, \( y(t)\in\mathbb{R}^p \) is the controlled output, \( r(t)\in\mathbb{R}^p \) is the reference, and \( e(t)\in\mathbb{R}^p \) is the tracking error. The matrix \( E \) identifies disturbance directions in the state equation; \( F \) identifies direct disturbance feedthrough to the output.

A tracking problem is therefore not fully defined by \( (A,B,C,D) \) alone. It also requires the regulated error definition. For example, if a system has several measured outputs but only one output must track the command, then \( C \) in the tracking model should select the regulated output, not necessarily the entire sensor vector.

flowchart TD
  R["reference r(t)"] --> ERR["tracking error e(t) = y(t) - r(t)"]
  D["disturbance d(t)"] --> PLANT["state-space plant: xdot = A x + B u + E d"]
  U["control input u(t)"] --> PLANT
  PLANT --> X["state x(t)"]
  PLANT --> Y["regulated output y(t)"]
  Y --> ERR
  ERR --> CTRL["controller or design objective"]
  X --> CTRL
  CTRL --> U
        

In state-space form, the tracking goal is written directly in terms of the regulated error:

\[ \lim_{t\to\infty} e(t) = 0. \]

This statement is stronger than saying the closed-loop poles are stable. Stability makes transients decay; tracking also requires the steady-state component of the output to match the reference. For a constant reference and constant disturbance, this requires an equilibrium satisfying both the plant equilibrium condition and the output matching condition.

If the closed-loop input approaches a constant value \( \bar{u} \) and the state approaches \( \bar{x} \), then exact constant tracking requires

\[ \begin{aligned} 0 &= A\bar{x} + B\bar{u} + E\bar{d}, \\ \bar{r} &= C\bar{x} + D\bar{u} + F\bar{d}. \end{aligned} \]

Equivalently, the unknown steady-state pair \( (\bar{x},\bar{u}) \) must solve the stacked linear equation

\[ \begin{bmatrix} A & B \\ C & D \end{bmatrix} \begin{bmatrix} \bar{x} \\ \bar{u} \end{bmatrix} = \begin{bmatrix} -E\bar{d} \\ \bar{r}-F\bar{d} \end{bmatrix}. \]

This equation is a feasibility test for the desired output. If it has no solution, then no constant input can make the specified output exactly equal to the specified constant reference under the specified constant disturbance. Feedback may reduce the error, but it cannot remove an algebraic infeasibility in the steady-state equations.

Define the steady-state tracking matrix \( \mathcal{M} \) by

\[ \mathcal{M} = \begin{bmatrix} A & B \\ C & D \end{bmatrix}. \]

For given \( \bar{r} \) and \( \bar{d} \), exact constant tracking is feasible if and only if the right-hand side belongs to the column space of \( \mathcal{M} \):

\[ \begin{bmatrix} -E\bar{d} \\ \bar{r}-F\bar{d} \end{bmatrix} \in \operatorname{im}(\mathcal{M}). \]

In numerical computations, this condition can be tested by the residual of a least-squares solution:

\[ \rho = \left\| \mathcal{M}\mathcal{M}^\dagger \begin{bmatrix} -E\bar{d} \\ \bar{r}-F\bar{d} \end{bmatrix} - \begin{bmatrix} -E\bar{d} \\ \bar{r}-F\bar{d} \end{bmatrix} \right\|. \]

The reference is feasible precisely when \( \rho=0 \) in exact arithmetic. In finite precision arithmetic, feasibility is judged by whether \( \rho \) is small relative to the norms of \( \mathcal{M} \) and the right-hand side.

The number of controlled outputs also matters. If the plant has fewer independent input directions than regulated output constraints, then arbitrary constant tracking may be impossible. This is not a failure of pole placement; it is an output-reachability limitation at steady state.

Suppose the constant tracking equilibrium \( (\bar{x},\bar{u}) \) exists. Define deviation variables

\[ \tilde{x} = x - \bar{x}, \qquad \tilde{u} = u - \bar{u}, \qquad \tilde{y} = y - \bar{r}. \]

Substituting \( x=\tilde{x}+\bar{x} \) and \( u=\tilde{u}+\bar{u} \) into the plant gives

\[ \begin{aligned} \dot{\tilde{x} } &= A(\tilde{x}+\bar{x})+B(\tilde{u}+\bar{u})+E\bar{d} \\ &= A\tilde{x}+B\tilde{u} +\left(A\bar{x}+B\bar{u}+E\bar{d}\right) \\ &= A\tilde{x}+B\tilde{u}. \end{aligned} \]

Similarly, the error output becomes

\[ \begin{aligned} e &= y-\bar{r} \\ &= C(\tilde{x}+\bar{x})+D(\tilde{u}+\bar{u}) +F\bar{d}-\bar{r} \\ &= C\tilde{x}+D\tilde{u}. \end{aligned} \]

Thus, once the correct equilibrium is identified, local tracking of a constant reference reduces to regulation of the deviation system:

\[ \dot{\tilde{x} } = A\tilde{x} + B\tilde{u}, \qquad e = C\tilde{x} + D\tilde{u}. \]

This reduction is central: the regulator problem is not independent of the tracking target, because the equilibrium \( (\bar{x},\bar{u}) \) changes with the reference and disturbance.

The reference signal does not usually enter the physical plant. It enters the control objective through the error \( e=y-r \). The disturbance, in contrast, may enter either the state equation, the output equation, or both:

\[ d(t) \mapsto \begin{cases} E d(t) & \text{state disturbance}, \\ F d(t) & \text{output disturbance}, \\ E d(t) \text{ and } F d(t) & \text{combined disturbance}. \end{cases} \]

A disturbance is often called matched when it enters through the same subspace as the control input, meaning \( \operatorname{im}(E)\subseteq\operatorname{im}(B) \). In that case a suitable constant control offset can directly counteract a constant disturbance in the state equation. If the disturbance is mismatched, the state can still be influenced by feedback, but the disturbance cannot be algebraically canceled at the input channel.

flowchart TD
  A["Choose regulated output y"] --> B["Define tracking error e = y - r"]
  B --> C["Specify disturbance path: state, output, or both"]
  C --> D["Write stacked steady-state equations"]
  D --> E["Check feasibility by rank or residual"]
  E --> F["Define deviation variables around equilibrium"]
  F --> G["Reduce constant tracking to deviation regulation"]
        

Constant references are only the simplest tracking signals. Step, ramp, sinusoidal, and disturbance signals can be represented by an autonomous signal generator called an exosystem:

\[ \dot{w} = S w, \qquad r = Q w, \qquad d = P w. \]

The plant and exosystem together form the augmented open-loop model

\[ \begin{bmatrix} \dot{x} \\ \dot{w} \end{bmatrix} = \begin{bmatrix} A & E P \\ 0 & S \end{bmatrix} \begin{bmatrix} x \\ w \end{bmatrix} + \begin{bmatrix} B \\ 0 \end{bmatrix} u, \qquad e = \begin{bmatrix} C & F P-Q \end{bmatrix} \begin{bmatrix} x \\ w \end{bmatrix} +D u. \]

The formulation asks whether there exist matrices \( \Pi \) and \( \Gamma \) such that the state and input can remain on the reference-dependent manifold

\[ x=\Pi w, \qquad u=\Gamma w. \]

Substitution gives the regulator equations

\[ \begin{aligned} A\Pi + B\Gamma + E P &= \Pi S, \\ C\Pi + D\Gamma + F P &= Q. \end{aligned} \]

In this lesson these equations are used only as a formulation device. Later lessons will discuss how these equations lead to feedforward tracking gains and how the internal model principle explains robust asymptotic tracking.

Chapters 22 through 26 introduced state feedback and integral action. Those ideas become tracking tools only after the error channel is precisely defined. A common constant-reference tracking structure is

\[ u = \bar{u} - K(x-\bar{x}), \]

where \( (\bar{x},\bar{u}) \) comes from the steady-state tracking equations. The resulting deviation dynamics are

\[ \dot{\tilde{x} } = (A-BK)\tilde{x}, \qquad e=(C-DK)\tilde{x}. \]

If \( A-BK \) is asymptotically stable, then \( \tilde{x}(t) \) converges to zero and therefore \( e(t) \) converges to zero for the modeled constant reference and disturbance. If the disturbance or plant parameters are uncertain, the equilibrium offset may be wrong; then integral action or an internal model is required to recover robust zero steady-state error.

Thus, the formulation separates three questions:

1. Is the reference/disturbance pair feasible at steady state?
2. Can the deviation dynamics be stabilized by feedback?
3. Is the controller robust to unknown constant or dynamic exogenous signals?

Consider the second-order plant

\[ \dot{x} = \begin{bmatrix} 0 & 1 \\ -2 & -0.8 \end{bmatrix}x + \begin{bmatrix} 0 \\ 1 \end{bmatrix}u + \begin{bmatrix} 0 \\ 1 \end{bmatrix}d, \qquad y = \begin{bmatrix} 1 & 0 \end{bmatrix}x. \]

Let the desired constant reference be \( \bar{r}=1 \) and the constant disturbance be \( \bar{d}=0.2 \). The steady-state equations are

\[ \begin{bmatrix} 0 & 1 & 0 \\ -2 & -0.8 & 1 \\ 1 & 0 & 0 \end{bmatrix} \begin{bmatrix} \bar{x}_1 \\ \bar{x}_2 \\ \bar{u} \end{bmatrix} = \begin{bmatrix} 0 \\ -0.2 \\ 1 \end{bmatrix}. \]

Solving gives \( \bar{x}_1=1 \), \( \bar{x}_2=0 \), and \( \bar{u}=1.8 \). With a stabilizing deviation gain \( K=\begin{bmatrix}4&2.2\end{bmatrix} \), the implemented input is

\[ u = 1.8 - \begin{bmatrix}4&2.2\end{bmatrix} \left(x- \begin{bmatrix}1\\0\end{bmatrix} \right). \]

This example is intentionally simple: the important lesson is not how the gain was chosen, but how the reference, disturbance, equilibrium, and error are written in state-space form.

Chapter27_Lesson1.py computes the stacked equilibrium equation, checks the residual, and simulates the deviation-feedback tracking formulation.

# Chapter27_Lesson1.py
# Formulating Tracking Problems in State-Space Form
#
# This script formulates a continuous-time tracking problem
#
#   x_dot = A x + B u + E d
#   y     = C x + D u + F d
#   e     = y - r
#
# It computes the steady-state equations for a constant reference r and
# constant disturbance d, checks reference feasibility, and simulates a
# stabilizing deviation-feedback law around the computed equilibrium.
#
# Required libraries:
#   pip install numpy scipy matplotlib

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt


def equilibrium_for_constant_reference(A, B, C, D, E, F, r, d):
    """Solve the constant tracking equilibrium equations.

    At steady state, x_dot = 0 and e = y - r = 0, hence

        A x_bar + B u_bar + E d = 0
        C x_bar + D u_bar + F d = r

    The stacked equation is solved in the least-squares sense so that
    underactuated or incompatible references can be detected by the residual.
    """
    A = np.asarray(A, dtype=float)
    B = np.asarray(B, dtype=float)
    C = np.asarray(C, dtype=float)
    D = np.asarray(D, dtype=float)
    E = np.asarray(E, dtype=float)
    F = np.asarray(F, dtype=float)
    r = np.asarray(r, dtype=float).reshape((-1, 1))
    d = np.asarray(d, dtype=float).reshape((-1, 1))

    M = np.block([[A, B],
                  [C, D]])
    rhs = np.vstack((-E @ d, r - F @ d))

    theta, *_ = np.linalg.lstsq(M, rhs, rcond=None)
    residual = M @ theta - rhs

    n = A.shape[0]
    x_bar = theta[:n]
    u_bar = theta[n:]
    return x_bar, u_bar, residual, M, rhs


def simulate_tracking_example():
    # Second-order plant: position output, force input, matched constant disturbance.
    A = np.array([[0.0, 1.0],
                  [-2.0, -0.8]])
    B = np.array([[0.0],
                  [1.0]])
    C = np.array([[1.0, 0.0]])
    D = np.array([[0.0]])
    E = np.array([[0.0],
                  [1.0]])
    F = np.array([[0.0]])

    r = np.array([1.0])      # constant reference position
    d = np.array([0.2])      # constant disturbance acceleration

    x_bar, u_bar, residual, M, rhs = equilibrium_for_constant_reference(
        A, B, C, D, E, F, r, d
    )

    print("Stacked equilibrium matrix [A B; C D]:")
    print(M)
    print("\nRight-hand side [-E d; r - F d]:")
    print(rhs)
    print("\nComputed x_bar:")
    print(x_bar)
    print("Computed u_bar:")
    print(u_bar)
    print("Equilibrium residual norm:", np.linalg.norm(residual))

    # The gain K is not designed in this lesson. It is used only to close
    # the loop around the formulated tracking equilibrium.
    K = np.array([[4.0, 2.2]])

    def closed_loop(t, x_flat):
        x = x_flat.reshape((-1, 1))
        u = u_bar - K @ (x - x_bar)
        x_dot = A @ x + B @ u + E @ d.reshape((-1, 1))
        return x_dot.ravel()

    t_span = (0.0, 8.0)
    t_eval = np.linspace(t_span[0], t_span[1], 801)
    x0 = np.array([0.0, 0.0])

    sol = solve_ivp(closed_loop, t_span, x0, t_eval=t_eval, rtol=1e-9, atol=1e-11)
    y = (C @ sol.y).ravel()
    e = y - r[0]

    print("\nFinal output y(T):", y[-1])
    print("Final tracking error e(T):", e[-1])

    plt.figure()
    plt.plot(sol.t, y, label="y(t)")
    plt.plot(sol.t, r[0] * np.ones_like(sol.t), "--", label="r")
    plt.plot(sol.t, e, label="e(t) = y(t) - r")
    plt.xlabel("time [s]")
    plt.ylabel("signals")
    plt.grid(True)
    plt.legend()
    plt.title("Chapter 27 Lesson 1: Tracking formulation simulation")
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    simulate_tracking_example()

Chapter27_Lesson1.cpp implements the same formulation from scratch, including a small Gaussian-elimination solver and RK4 integration.

// Chapter27_Lesson1.cpp
// Formulating Tracking Problems in State-Space Form
//
// Compile:
//   g++ -std=c++17 Chapter27_Lesson1.cpp -o Chapter27_Lesson1
//
// This from-scratch implementation solves the constant tracking equilibrium
// equations for a second-order plant and simulates a deviation-feedback
// tracking formulation using RK4 integration.

#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>

using Vec2 = std::array<double, 2>;
using Vec3 = std::array<double, 3>;

struct Equilibrium {
    double x1_bar;
    double x2_bar;
    double u_bar;
    double residual_norm;
};

Vec3 solve3x3(double M[3][3], Vec3 b) {
    // Gaussian elimination with partial pivoting.
    for (int k = 0; k < 3; ++k) {
        int pivot = k;
        double best = std::abs(M[k][k]);
        for (int i = k + 1; i < 3; ++i) {
            if (std::abs(M[i][k]) > best) {
                best = std::abs(M[i][k]);
                pivot = i;
            }
        }
        if (pivot != k) {
            for (int j = 0; j < 3; ++j) std::swap(M[k][j], M[pivot][j]);
            std::swap(b[k], b[pivot]);
        }

        const double diag = M[k][k];
        if (std::abs(diag) < 1e-12) {
            throw std::runtime_error("Singular equilibrium matrix.");
        }

        for (int j = k; j < 3; ++j) M[k][j] /= diag;
        b[k] /= diag;

        for (int i = 0; i < 3; ++i) {
            if (i == k) continue;
            const double factor = M[i][k];
            for (int j = k; j < 3; ++j) M[i][j] -= factor * M[k][j];
            b[i] -= factor * b[k];
        }
    }
    return b;
}

Equilibrium compute_equilibrium(double reference, double disturbance) {
    // Plant:
    //   x1_dot = x2
    //   x2_dot = -2 x1 - 0.8 x2 + u + d
    //   y      = x1
    //
    // Constant tracking equilibrium:
    //   0 = x2
    //   0 = -2 x1 - 0.8 x2 + u + d
    //   r = x1

    double M[3][3] = {
        {0.0, 1.0, 0.0},
        {-2.0, -0.8, 1.0},
        {1.0, 0.0, 0.0}
    };
    Vec3 rhs = {0.0, -disturbance, reference};

    Vec3 theta = solve3x3(M, rhs);
    const double x1 = theta[0];
    const double x2 = theta[1];
    const double u = theta[2];

    const double res1 = x2;
    const double res2 = -2.0 * x1 - 0.8 * x2 + u + disturbance;
    const double res3 = x1 - reference;
    const double norm = std::sqrt(res1 * res1 + res2 * res2 + res3 * res3);

    return {x1, x2, u, norm};
}

Vec2 dynamics(const Vec2& x, const Equilibrium& eq, double disturbance) {
    const double k1 = 4.0;
    const double k2 = 2.2;
    const double u = eq.u_bar - k1 * (x[0] - eq.x1_bar) - k2 * (x[1] - eq.x2_bar);

    Vec2 dx;
    dx[0] = x[1];
    dx[1] = -2.0 * x[0] - 0.8 * x[1] + u + disturbance;
    return dx;
}

Vec2 add_scaled(const Vec2& x, const Vec2& k, double scale) {
    return {x[0] + scale * k[0], x[1] + scale * k[1]};
}

Vec2 rk4_step(const Vec2& x, double h, const Equilibrium& eq, double disturbance) {
    Vec2 k1 = dynamics(x, eq, disturbance);
    Vec2 k2 = dynamics(add_scaled(x, k1, 0.5 * h), eq, disturbance);
    Vec2 k3 = dynamics(add_scaled(x, k2, 0.5 * h), eq, disturbance);
    Vec2 k4 = dynamics(add_scaled(x, k3, h), eq, disturbance);

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

int main() {
    const double reference = 1.0;
    const double disturbance = 0.2;
    Equilibrium eq = compute_equilibrium(reference, disturbance);

    std::cout << std::fixed << std::setprecision(8);
    std::cout << "x_bar = [" << eq.x1_bar << ", " << eq.x2_bar << "]\n";
    std::cout << "u_bar = " << eq.u_bar << "\n";
    std::cout << "equilibrium residual norm = " << eq.residual_norm << "\n";

    Vec2 x = {0.0, 0.0};
    const double h = 0.001;
    const double tf = 8.0;
    const int steps = static_cast<int>(tf / h);

    for (int i = 0; i < steps; ++i) {
        x = rk4_step(x, h, eq, disturbance);
    }

    const double y = x[0];
    const double e = y - reference;

    std::cout << "final y(T) = " << y << "\n";
    std::cout << "final e(T) = " << e << "\n";
    return 0;
}

Chapter27_Lesson1.java mirrors the C++ implementation using plain Java arrays.

// Chapter27_Lesson1.java
// Formulating Tracking Problems in State-Space Form
//
// Compile and run:
//   javac Chapter27_Lesson1.java
//   java Chapter27_Lesson1
//
// This from-scratch implementation solves the constant tracking equilibrium
// equations for a second-order plant and simulates a deviation-feedback
// tracking formulation using RK4 integration.

public class Chapter27_Lesson1 {
    static class Equilibrium {
        double x1Bar;
        double x2Bar;
        double uBar;
        double residualNorm;

        Equilibrium(double x1Bar, double x2Bar, double uBar, double residualNorm) {
            this.x1Bar = x1Bar;
            this.x2Bar = x2Bar;
            this.uBar = uBar;
            this.residualNorm = residualNorm;
        }
    }

    static double[] solve3x3(double[][] M, double[] b) {
        for (int k = 0; k < 3; k++) {
            int pivot = k;
            double best = Math.abs(M[k][k]);
            for (int i = k + 1; i < 3; i++) {
                if (Math.abs(M[i][k]) > best) {
                    best = Math.abs(M[i][k]);
                    pivot = i;
                }
            }

            if (pivot != k) {
                double[] tmpRow = M[k];
                M[k] = M[pivot];
                M[pivot] = tmpRow;

                double tmp = b[k];
                b[k] = b[pivot];
                b[pivot] = tmp;
            }

            double diag = M[k][k];
            if (Math.abs(diag) < 1e-12) {
                throw new RuntimeException("Singular equilibrium matrix.");
            }

            for (int j = k; j < 3; j++) {
                M[k][j] /= diag;
            }
            b[k] /= diag;

            for (int i = 0; i < 3; i++) {
                if (i == k) {
                    continue;
                }
                double factor = M[i][k];
                for (int j = k; j < 3; j++) {
                    M[i][j] -= factor * M[k][j];
                }
                b[i] -= factor * b[k];
            }
        }
        return b;
    }

    static Equilibrium computeEquilibrium(double reference, double disturbance) {
        // Plant:
        //   x1_dot = x2
        //   x2_dot = -2 x1 - 0.8 x2 + u + d
        //   y      = x1
        //
        // Constant tracking equilibrium:
        //   0 = x2
        //   0 = -2 x1 - 0.8 x2 + u + d
        //   r = x1

        double[][] M = {
            {0.0, 1.0, 0.0},
            {-2.0, -0.8, 1.0},
            {1.0, 0.0, 0.0}
        };
        double[] rhs = {0.0, -disturbance, reference};

        double[] theta = solve3x3(M, rhs);
        double x1 = theta[0];
        double x2 = theta[1];
        double u = theta[2];

        double res1 = x2;
        double res2 = -2.0 * x1 - 0.8 * x2 + u + disturbance;
        double res3 = x1 - reference;
        double norm = Math.sqrt(res1 * res1 + res2 * res2 + res3 * res3);

        return new Equilibrium(x1, x2, u, norm);
    }

    static double[] dynamics(double[] x, Equilibrium eq, double disturbance) {
        double k1 = 4.0;
        double k2 = 2.2;
        double u = eq.uBar - k1 * (x[0] - eq.x1Bar) - k2 * (x[1] - eq.x2Bar);

        return new double[] {
            x[1],
            -2.0 * x[0] - 0.8 * x[1] + u + disturbance
        };
    }

    static double[] addScaled(double[] x, double[] k, double scale) {
        return new double[] {x[0] + scale * k[0], x[1] + scale * k[1]};
    }

    static double[] rk4Step(double[] x, double h, Equilibrium eq, double disturbance) {
        double[] k1 = dynamics(x, eq, disturbance);
        double[] k2 = dynamics(addScaled(x, k1, 0.5 * h), eq, disturbance);
        double[] k3 = dynamics(addScaled(x, k2, 0.5 * h), eq, disturbance);
        double[] k4 = dynamics(addScaled(x, k3, h), eq, disturbance);

        return new double[] {
            x[0] + h * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0,
            x[1] + h * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0
        };
    }

    public static void main(String[] args) {
        double reference = 1.0;
        double disturbance = 0.2;
        Equilibrium eq = computeEquilibrium(reference, disturbance);

        System.out.printf("x_bar = [%.8f, %.8f]%n", eq.x1Bar, eq.x2Bar);
        System.out.printf("u_bar = %.8f%n", eq.uBar);
        System.out.printf("equilibrium residual norm = %.8e%n", eq.residualNorm);

        double[] x = {0.0, 0.0};
        double h = 0.001;
        double tf = 8.0;
        int steps = (int) (tf / h);

        for (int i = 0; i < steps; i++) {
            x = rk4Step(x, h, eq, disturbance);
        }

        double y = x[0];
        double e = y - reference;
        System.out.printf("final y(T) = %.8f%n", y);
        System.out.printf("final e(T) = %.8e%n", e);
    }
}

Chapter27_Lesson1.m runs the MATLAB simulation and, if Simulink is available, builds a small deviation-state-space model.

% Chapter27_Lesson1.m
% Formulating Tracking Problems in State-Space Form
%
% This script formulates a continuous-time tracking problem
%
%   x_dot = A*x + B*u + E*d
%   y     = C*x + D*u + F*d
%   e     = y - r
%
% It computes constant-reference equilibrium equations, checks residual
% feasibility, and simulates a stabilizing deviation-feedback law.
%
% Toolboxes:
%   - The numerical simulation works in base MATLAB.
%   - The optional Simulink section requires Simulink.

clear; clc; close all;

A = [0 1; -2 -0.8];
B = [0; 1];
C = [1 0];
D = 0;
E = [0; 1];
F = 0;

r = 1.0;       % constant reference
d = 0.2;       % constant disturbance

% Steady-state tracking equations:
%   A*x_bar + B*u_bar + E*d = 0
%   C*x_bar + D*u_bar + F*d = r
M = [A B; C D];
rhs = [-E*d; r - F*d];

theta = pinv(M) * rhs;
x_bar = theta(1:2);
u_bar = theta(3);
residual = M * theta - rhs;

disp('Stacked equilibrium matrix [A B; C D]:');
disp(M);
disp('Computed x_bar:');
disp(x_bar);
disp('Computed u_bar:');
disp(u_bar);
fprintf('Equilibrium residual norm = %.8e\n', norm(residual));

% The gain K is not the topic of this lesson; it only stabilizes the
% formulated tracking equilibrium for demonstration.
K = [4.0 2.2];

closed_loop = @(t, x) A*x + B*(u_bar - K*(x - x_bar)) + E*d;

tspan = [0 8];
x0 = [0; 0];
[t, x] = ode45(closed_loop, tspan, x0);

y = (C*x')';
e = y - r;

fprintf('Final y(T) = %.8f\n', y(end));
fprintf('Final e(T) = %.8e\n', e(end));

figure;
plot(t, y, 'LineWidth', 1.5); hold on;
plot(t, r*ones(size(t)), '--', 'LineWidth', 1.2);
plot(t, e, 'LineWidth', 1.2);
grid on;
xlabel('time [s]');
ylabel('signals');
legend('y(t)', 'r', 'e(t)=y(t)-r', 'Location', 'best');
title('Chapter 27 Lesson 1: Tracking formulation simulation');

% Optional Simulink construction:
% The Simulink interpretation is the deviation model
%
%   x_tilde_dot = (A - B*K)*x_tilde
%   y_tilde     = C*x_tilde
%
% where x_tilde = x - x_bar and y_tilde = y - r.
if exist('simulink', 'file') == 4
    modelName = 'Chapter27_Lesson1_Simulink_Formulation';
    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    end

    new_system(modelName);
    open_system(modelName);

    Acl = A - B*K;
    Bcl = zeros(2, 1);
    Ccl = C;
    Dcl = 0;

    add_block('simulink/Sources/Constant', [modelName '/zero input'], ...
        'Value', '0', 'Position', [80 80 130 110]);
    add_block('simulink/Continuous/State-Space', [modelName '/deviation state-space'], ...
        'A', mat2str(Acl), 'B', mat2str(Bcl), ...
        'C', mat2str(Ccl), 'D', mat2str(Dcl), ...
        'X0', mat2str(x0 - x_bar), 'Position', [190 65 340 125]);
    add_block('simulink/Sinks/Scope', [modelName '/tracking error scope'], ...
        'Position', [410 75 460 115]);

    add_line(modelName, 'zero input/1', 'deviation state-space/1');
    add_line(modelName, 'deviation state-space/1', 'tracking error scope/1');

    set_param(modelName, 'StopTime', '8');
    save_system(modelName);
    disp(['Optional Simulink model created: ' modelName '.slx']);
else
    disp('Simulink is not available; skipped optional Simulink model creation.');
end

Chapter27_Lesson1.nb provides the same computation in a Mathematica notebook expression.

Notebook[{
Cell["Chapter27_Lesson1.nb", "Title"],
Cell["Formulating Tracking Problems in State-Space Form", "Subtitle"],
Cell["This notebook formulates xdot = A x + B u + E d, y = C x + D u + F d, e = y - r; computes the constant tracking equilibrium; and simulates a deviation-feedback tracking formulation.", "Text"],
Cell[BoxData["ClearAll[\"Global`*\"];"], "Input"],
Cell[BoxData["A = { {0, 1}, {-2, -0.8} };\nB = { {0}, {1} };\nCmat = { {1, 0} };\nDmat = { {0} };\nEcon = { {0}, {1} };\nFmat = { {0} };\nr = { {1.0} };\nd = { {0.2} };"], "Input"],
Cell[BoxData["M = ArrayFlatten[{ {A, B}, {Cmat, Dmat} }];\nrhs = Join[-Econ.d, r - Fmat.d];\ntheta = LeastSquares[M, rhs];\nxbar = theta[[1 ;; 2]];\nubar = theta[[3 ;; 3]];\nresidual = M.theta - rhs;\n{MatrixForm[M], MatrixForm[xbar], MatrixForm[ubar], Norm[residual]}"], "Input"],
Cell[BoxData["K = { {4.0, 2.2} };\nx1bar = xbar[[1, 1]];\nx2bar = xbar[[2, 1]];\nubarScalar = ubar[[1, 1]];\ndScalar = d[[1, 1]];"], "Input"],
Cell[BoxData["sol = NDSolve[\n  {\n    x1'[t] == x2[t],\n    x2'[t] == -2*x1[t] - 0.8*x2[t] + ubarScalar - K[[1, 1]]*(x1[t] - x1bar) - K[[1, 2]]*(x2[t] - x2bar) + dScalar,\n    x1[0] == 0,\n    x2[0] == 0\n  },\n  {x1, x2}, {t, 0, 8}\n];"], "Input"],
Cell[BoxData["Plot[\n  Evaluate[{x1[t] /. sol[[1]], 1, (x1[t] /. sol[[1]]) - 1}],\n  {t, 0, 8},\n  PlotLegends -> {\"y(t)\", \"r\", \"e(t)=y(t)-r\"},\n  AxesLabel -> {\"time [s]\", \"signals\"},\n  PlotLabel -> \"Chapter 27 Lesson 1: Tracking formulation simulation\"\n]"], "Input"],
Cell[BoxData["finalY = x1[8] /. sol[[1]];\nfinalError = finalY - 1;\n{finalY, finalError}"], "Input"]
}]

Problem 1 (Constant tracking feasibility): Consider \( \dot{x}=Ax+Bu \) and \( y=Cx \). Derive the algebraic condition for exact tracking of a constant reference \( \bar{r} \).

Solution: Exact constant tracking requires a constant pair \( (\bar{x},\bar{u}) \) such that \( \dot{x}=0 \) and \( y=\bar{r} \). Therefore,

\[ A\bar{x}+B\bar{u}=0, \qquad C\bar{x}=\bar{r}. \]

In stacked form,

\[ \begin{bmatrix} A&B\\C&0 \end{bmatrix} \begin{bmatrix} \bar{x}\\\bar{u} \end{bmatrix} = \begin{bmatrix}0\\\bar{r}\end{bmatrix}. \]

Hence the reference is feasible if and only if \( \begin{bmatrix}0\\\bar{r}\end{bmatrix} \) belongs to the image of \( \begin{bmatrix} A&B\\C&0 \end{bmatrix} \).

Problem 2 (Deviation model proof): Suppose \( (\bar{x},\bar{u}) \) satisfies \( A\bar{x}+B\bar{u}+E\bar{d}=0 \) and \( C\bar{x}+D\bar{u}+F\bar{d}=\bar{r} \). Prove that \( \tilde{x}=x-\bar{x} \) and \( \tilde{u}=u-\bar{u} \) satisfy \( \dot{\tilde{x} }=A\tilde{x}+B\tilde{u} \) for constant \( \bar{d} \).

Solution: Since \( \bar{x} \) is constant, \( \dot{\tilde{x} }=\dot{x} \). Substituting \( x=\tilde{x}+\bar{x} \) and \( u=\tilde{u}+\bar{u} \) gives

\[ \dot{\tilde{x} } =A\tilde{x}+B\tilde{u} +A\bar{x}+B\bar{u}+E\bar{d} =A\tilde{x}+B\tilde{u}. \]

The equilibrium term vanishes by assumption.

Problem 3 (Incompatible output constraint): Let \( A=0 \), \( B=1 \), \( C=0 \), and \( D=0 \). Can the output track \( \bar{r}=1 \)?

Solution: The output is \( y=Cx+Du=0 \) for all states and inputs. Therefore \( y=1 \) is impossible. The stacked equation is

\[ \begin{bmatrix}0&1\\0&0\end{bmatrix} \begin{bmatrix}\bar{x}\\\bar{u}\end{bmatrix} = \begin{bmatrix}0\\1\end{bmatrix}, \]

whose second row reads \( 0=1 \). Thus the reference is algebraically infeasible.

Problem 4 (Exosystem dimensions): Suppose \( x\in\mathbb{R}^n \), \( u\in\mathbb{R}^m \), \( e\in\mathbb{R}^p \), and \( w\in\mathbb{R}^s \). What are the dimensions of \( \Pi \) and \( \Gamma \) in the regulator equations?

Solution: The manifold equations are \( x=\Pi w \) and \( u=\Gamma w \). Therefore \( \Pi\in\mathbb{R}^{n\times s} \) and \( \Gamma\in\mathbb{R}^{m\times s} \). The first regulator equation \( A\Pi+B\Gamma+EP=\Pi S \) has dimension \( n\times s \), while the second equation \( C\Pi+D\Gamma+FP=Q \) has dimension \( p\times s \).

Problem 5 (Matched constant disturbance): Assume \( E=B L \) for some matrix \( L \). Show how a known constant disturbance \( \bar{d} \) can be absorbed into a shifted input.

Solution: Since \( E=B L \),

\[ Ax+Bu+E\bar{d} = Ax+B\left(u+L\bar{d}\right). \]

Defining \( u_a=u+L\bar{d} \) gives \( \dot{x}=Ax+Bu_a \). Thus a known matched constant disturbance can be represented as an input offset. If \( \bar{d} \) is unknown, feedback with integral action or an internal model is needed for robust rejection.

This lesson established the state-space formulation of reference tracking and disturbance rejection. The key objects are the regulated output, tracking error, disturbance channels, steady-state feasibility equation, deviation variables, and exosystem representation. These objects define the problem before any controller is designed. In the next lesson, the steady-state equations will be used to construct feedforward gains for reference tracking.