In the previous two lessons, transmission zeros were introduced through the Rosenbrock system matrix and computed from a realization \( (A,B,C,D) \). In this lesson we ask a deeper internal question:

If the output is forced to remain zero for all time, what internal motion of the state is still possible?

For a continuous-time LTI system

\[ \dot{x}(t)=Ax(t)+Bu(t), \qquad y(t)=Cx(t)+Du(t), \]

the zero dynamics are the internal dynamics compatible with the constraint \( y(t)\equiv 0 \). They reveal behavior that may be invisible at the output. This distinction is essential because input-output stability does not automatically imply internal stability unless the realization is minimal and the zero dynamics are stable.

flowchart TD
  A["State-space system"] --> B["Impose output constraint: y(t) = 0"]
  B --> C["Find compatible state-input trajectories"]
  C --> D["Internal motion remains"]
  D --> E["Zero dynamics"]
  E --> F["Stable internal behavior: \nminimum phase"]
  E --> G["Unstable internal behavior: \nnonminimum phase"]
        

A trajectory pair \( (x(t),u(t)) \) is called output-nulling if it satisfies both the state equation and the zero-output condition:

\[ \dot{x}(t)=Ax(t)+Bu(t), \qquad Cx(t)+Du(t)=0 \quad \text{for all } t\ge 0. \]

The set of initial states from which such a trajectory exists is called an output-nulling controlled invariant subspace. In geometric control language, a subspace \( \mathcal{V} \subseteq \mathbb{R}^n \) is output-nulling controlled invariant if for every \( x\in\mathcal{V} \), one can choose an input \( u \) such that

\[ Ax+Bu \in \mathcal{V}, \qquad Cx+Du=0. \]

Therefore, zero dynamics are not simply the dynamics of \( A \). They are the dynamics of the state constrained to an output-nulling manifold and driven by the special input that keeps the output zero.

For a transmission zero \( z \in \mathbb{C} \), the Rosenbrock matrix loses rank:

\[ \operatorname{rank} \begin{bmatrix} zI-A & -B\\ C & D \end{bmatrix} < n+\min(m,p). \]

Equivalently, for a nonzero vector \( x_z \) and some vector \( u_z \), one has

\[ \begin{bmatrix} zI-A & -B\\ C & D \end{bmatrix} \begin{bmatrix} x_z\\ u_z \end{bmatrix} = \begin{bmatrix} 0\\ 0 \end{bmatrix}. \]

This means \( x(t)=e^{zt}x_z \) and \( u(t)=e^{zt}u_z \) generate \( y(t)=0 \). Hence, invariant zeros are exponential modes of the zero-output internal motion.

Consider a strictly proper SISO realization

\[ \dot{x}=Ax+bu,\qquad y=cx,\qquad D=0. \]

Suppose the system has relative degree \( r \). This means

\[ cb=0,\quad cAb=0,\quad \dots,\quad cA^{r-2}b=0, \qquad cA^{r-1}b\ne 0. \]

Differentiating the output gives

\[ \begin{aligned} y &= cx,\\ \dot{y} &= cAx+cbu,\\ \ddot{y} &= cA^2x+cAbu,\\ &\vdots\\ y^{(r)} &= cA^r x+cA^{r-1}b\,u. \end{aligned} \]

To keep the output zero, impose \( y=\dot{y}=\cdots=y^{(r-1)}=0 \) and choose \( u \) from \( y^{(r)}=0 \):

\[ u = -\frac{cA^r x}{cA^{r-1}b}. \]

Substituting this feedback into the state equation produces motion constrained to the zero-output manifold:

\[ \dot{x} = \left( A-b\frac{cA^r}{cA^{r-1}b} \right)x. \]

The state is not free to move in all directions; it must satisfy

\[ cx=0,\quad cAx=0,\quad \dots,\quad cA^{r-1}x=0. \]

The restriction of the above closed internal dynamics to this manifold is the SISO zero dynamics. If the system is minimal and has order \( n \), then the dimension of the zero dynamics is \( n-r \).

For a SISO system with relative degree \( r \), one can locally transform the state into external coordinates \( \xi \) and internal coordinates \( \eta \):

\[ \xi_1=y,\quad \xi_2=\dot{y},\quad \dots,\quad \xi_r=y^{(r-1)}. \]

In these coordinates, the system takes the form

\[ \begin{aligned} \dot{\xi}_1 &= \xi_2,\\ \dot{\xi}_2 &= \xi_3,\\ &\vdots\\ \dot{\xi}_r &= a(\xi,\eta)+b(\xi,\eta)u,\\ \dot{\eta} &= q(\xi,\eta). \end{aligned} \]

The zero-output condition is \( \xi=0 \). Choosing \( u \) so that \( \dot{\xi}_r=0 \) yields

\[ u_0(\eta)=-\frac{a(0,\eta)}{b(0,\eta)}. \]

Therefore the zero dynamics are

\[ \dot{\eta}=q(0,\eta). \]

For linear systems, this internal system is linear:

\[ \dot{\eta}=A_0\eta, \]

and the eigenvalues of \( A_0 \) are precisely the finite transmission zeros of a minimal strictly proper SISO system.

Consider the third-order SISO system

\[ A= \begin{bmatrix} 0&1&0\\ 0&0&1\\ -6&-11&-6 \end{bmatrix}, \quad b= \begin{bmatrix} 0\\0\\1 \end{bmatrix}, \quad c= \begin{bmatrix} -1&1&0 \end{bmatrix}. \]

The transfer function is

\[ G(s)=c(sI-A)^{-1}b = \frac{s-1}{(s+1)(s+2)(s+3)}. \]

Thus the system has a transmission zero at \( z=+1 \). The poles are \( -1,-2,-3 \), so the open-loop internal matrix \( A \) is asymptotically stable. However, the zero is in the right half-plane.

The output is

\[ y=-x_1+x_2. \]

To impose \( y=0 \), we require \( x_2=x_1 \). Also,

\[ \dot{y}=-\dot{x}_1+\dot{x}_2=-x_2+x_3. \]

Therefore \( \dot{y}=0 \) requires \( x_3=x_1 \). On the output-nulling manifold,

\[ x= \begin{bmatrix} \eta\\ \eta\\ \eta \end{bmatrix}. \]

The state equation gives

\[ \dot{x}_1=x_2=\eta, \qquad \dot{x}_3=-6x_1-11x_2-6x_3+u=-23\eta+u. \]

Since \( x_3=x_1 \) must remain true, we impose \( \dot{x}_3=\dot{x}_1 \). Hence

\[ -23\eta+u=\eta \quad \Longrightarrow \quad u=24\eta. \]

Therefore the zero dynamics are

\[ \dot{\eta}=\eta. \]

This example demonstrates the core phenomenon: the output can be held identically zero while the internal state grows exponentially. The growth rate is the right-half-plane zero \( z=+1 \).

flowchart TD
  A["Output constraint: y = -x1 + x2 = 0"] --> B["Therefore x2 = x1"]
  B --> C["Derivative constraint: ydot = -x2 + x3 = 0"]
  C --> D["Therefore x3 = x1"]
  D --> E["State lies on manifold: x = eta times [1,1,1]"]
  E --> F["Choose input to keep manifold invariant"]
  F --> G["Zero dynamics: eta dot = eta"]
  G --> H["Internal state grows while output stays zero"]
        

Important note: With the realization above, direct calculation of the invariant-manifold input gives \( u=24\eta \). If an alternative companion convention or numerator normalization is used, the required output-nulling input may have a different coefficient, but the internal zero-dynamics eigenvalue remains \( +1 \).

For a continuous-time minimal SISO system, the zero dynamics are asymptotically stable if and only if all finite transmission zeros lie in the open left half-plane:

\[ \operatorname{Re}(z_i)<0 \quad \text{for all finite transmission zeros } z_i. \]

Such a system is called minimum phase. If at least one zero satisfies

\[ \operatorname{Re}(z_i)>0, \]

then the zero dynamics are unstable and the system is non-minimum phase.

The phrase "minimum phase" originates from frequency-domain theory, but the state-space interpretation is internal: stable zero dynamics mean that the hidden internal motion associated with zero output decays rather than grows.

For MIMO systems, zero dynamics are generally more subtle because the number of inputs and outputs may differ, and there may be multiple input directions capable of maintaining zero output. The invariant-zero condition is still expressed by the Rosenbrock matrix:

\[ \mathcal{R}(s)= \begin{bmatrix} sI-A & -B\\ C & D \end{bmatrix}. \]

A finite complex number \( z \) is an invariant zero if

\[ \operatorname{rank}\mathcal{R}(z) < \operatorname{normalrank}\mathcal{R}(s). \]

If there exist nonzero \( x_z \) and some \( u_z \) satisfying

\[ (zI-A)x_z-Bu_z=0,\qquad Cx_z+Du_z=0, \]

then the exponential state-input pair

\[ x(t)=e^{zt}x_z,\qquad u(t)=e^{zt}u_z \]

produces \( y(t)=0 \). Hence invariant zeros are internal output-nulling modes even for MIMO systems.

This Python script uses scipy.signal.ss2tf to compute the transfer numerator and denominator, extracts the zero, and simulates the output-nulling internal motion.

# Chapter21_Lesson3.py
import numpy as np
from scipy import signal
from scipy.integrate import solve_ivp

A = np.array([[0.0, 1.0, 0.0],
              [0.0, 0.0, 1.0],
              [-6.0, -11.0, -6.0]])

B = np.array([[0.0],
              [0.0],
              [1.0]])

C = np.array([[-1.0, 1.0, 0.0]])
D = np.array([[0.0]])

def transfer_zeros_from_ss(A, B, C, D):
    num, den = signal.ss2tf(A, B, C, D)
    num = np.trim_zeros(num[0], trim="f")
    den = np.trim_zeros(den, trim="f")
    zeros = np.roots(num)
    poles = np.roots(den)
    return zeros, poles, num, den

def output_nulling_feedback(x):
    # For the manifold x = eta*[1,1,1], invariance gives u = 24*eta.
    return 24.0 * x[0]

def closed_loop_zero_dynamics(t, x):
    u = output_nulling_feedback(x)
    return A @ x + B.flatten() * u

def main():
    zeros, poles, num, den = transfer_zeros_from_ss(A, B, C, D)

    print("Numerator coefficients:", num)
    print("Denominator coefficients:", den)
    print("Transmission zero(s):", zeros)
    print("Poles:", poles)

    eta0 = 0.02
    x0 = np.array([eta0, eta0, eta0])

    sol = solve_ivp(
        closed_loop_zero_dynamics,
        t_span=(0.0, 5.0),
        y0=x0,
        t_eval=np.linspace(0.0, 5.0, 101),
        rtol=1e-9,
        atol=1e-11,
    )

    y = np.array([(C @ sol.y[:, k])[0] for k in range(sol.y.shape[1])])
    eta_exact = eta0 * np.exp(sol.t)

    print("\nZero-output simulation")
    print("max |y(t)|:", np.max(np.abs(y)))
    print("eta(5) numerical:", sol.y[0, -1])
    print("eta(5) exact:", eta_exact[-1])
    print("interpretation: output is zero, internal state grows because zero = +1")

if __name__ == "__main__":
    main()
      

The C++ implementation uses Eigen and a fourth-order Runge-Kutta integrator to simulate the zero-output trajectory.

// Chapter21_Lesson3.cpp
#include <Eigen/Dense>
#include <cmath>
#include <iomanip>
#include <iostream>

using Vec3 = Eigen::Vector3d;
using Mat3 = Eigen::Matrix3d;

const Mat3 A = (Mat3() << 0.0, 1.0, 0.0,
                          0.0, 0.0, 1.0,
                         -6.0, -11.0, -6.0).finished();

const Vec3 B(0.0, 0.0, 1.0);
const Eigen::RowVector3d C(-1.0, 1.0, 0.0);

double outputNullingInput(const Vec3& x) {
    return 24.0 * x(0);
}

Vec3 f(const Vec3& x) {
    double u = outputNullingInput(x);
    return A * x + B * u;
}

Vec3 rk4Step(const Vec3& x, double h) {
    Vec3 k1 = f(x);
    Vec3 k2 = f(x + 0.5 * h * k1);
    Vec3 k3 = f(x + 0.5 * h * k2);
    Vec3 k4 = f(x + h * k3);
    return x + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
}

int main() {
    std::cout << "Known transfer function zero: z = +1\n";
    std::cout << "Known poles: -1, -2, -3\n\n";

    double eta0 = 0.02;
    Vec3 x(eta0, eta0, eta0);

    double h = 0.01;
    int steps = static_cast<int>(5.0 / h);

    double maxAbsY = 0.0;

    for (int k = 0; k <= steps; ++k) {
        double t = k * h;
        double u = outputNullingInput(x);
        double y = (C * x)(0);
        maxAbsY = std::max(maxAbsY, std::abs(y));

        if (k % 100 == 0) {
            std::cout << std::fixed << std::setprecision(6)
                      << "t=" << t
                      << ", eta=x1=" << x(0)
                      << ", u=" << u
                      << ", y=" << y << "\n";
        }

        if (k < steps) {
            x = rk4Step(x, h);
        }
    }

    std::cout << "\nmax |y(t)| = " << maxAbsY << "\n";
    std::cout << "eta(5) numerical = " << x(0) << "\n";
    std::cout << "eta(5) exact     = " << eta0 * std::exp(5.0) << "\n";
    return 0;
}
      

The Java version is written from scratch using arrays and RK4 integration. For larger projects, Apache Commons Math can be used for matrix and eigenvalue operations.

// Chapter21_Lesson3.java
public class Chapter21_Lesson3 {
    static final double[][] A = {
            {0.0, 1.0, 0.0},
            {0.0, 0.0, 1.0},
            {-6.0, -11.0, -6.0}
    };

    static final double[] B = {0.0, 0.0, 1.0};
    static final double[] C = {-1.0, 1.0, 0.0};

    static double outputNullingInput(double[] x) {
        return 24.0 * x[0];
    }

    static double[] dynamics(double[] x) {
        double u = outputNullingInput(x);
        double[] dx = new double[3];

        for (int i = 0; i < 3; i++) {
            dx[i] = B[i] * u;
            for (int j = 0; j < 3; j++) {
                dx[i] += A[i][j] * x[j];
            }
        }
        return dx;
    }

    static double output(double[] x) {
        return C[0] * x[0] + C[1] * x[1] + C[2] * x[2];
    }

    static double[] add(double[] x, double scale, double[] v) {
        return new double[] {
                x[0] + scale * v[0],
                x[1] + scale * v[1],
                x[2] + scale * v[2]
        };
    }

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

        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,
                x[2] + h * (k1[2] + 2.0 * k2[2] + 2.0 * k3[2] + k4[2]) / 6.0
        };
    }

    public static void main(String[] args) {
        double eta0 = 0.02;
        double[] x = {eta0, eta0, eta0};

        double h = 0.01;
        int steps = (int)(5.0 / h);
        double maxAbsY = 0.0;

        for (int k = 0; k <= steps; k++) {
            double y = output(x);
            maxAbsY = Math.max(maxAbsY, Math.abs(y));

            if (k % 100 == 0) {
                System.out.printf("t=%.2f, eta=x1=%.6f, u=%.6f, y=%.6e%n",
                        k * h, x[0], outputNullingInput(x), y);
            }

            if (k < steps) {
                x = rk4Step(x, h);
            }
        }

        System.out.printf("%nmax |y(t)| = %.6e%n", maxAbsY);
        System.out.printf("eta(5) numerical = %.9f%n", x[0]);
        System.out.printf("eta(5) exact     = %.9f%n", eta0 * Math.exp(5.0));
    }
}
      

MATLAB's Control System Toolbox directly computes poles, zeros, and transfer functions. In Simulink, the same model can be built using a State-Space block and a feedback block implementing the output-nulling input.

% Chapter21_Lesson3.m
clear; clc; close all;

A = [0 1 0;
     0 0 1;
    -6 -11 -6];

B = [0; 0; 1];
C = [-1 1 0];
D = 0;

sys = ss(A,B,C,D);
disp('Transfer function:')
G = tf(sys)

disp('Transmission zero(s):')
try
    z = tzero(sys)
catch
    z = zero(sys)
end

disp('Poles:')
p = pole(sys)

eta0 = 0.02;
x0 = [eta0; eta0; eta0];

odefun = @(t,x) A*x + B*(24*x(1));
tspan = linspace(0,5,501);
[t,x] = ode45(odefun, tspan, x0);

u = 24*x(:,1);
y = (C*x.').';

fprintf('max |y(t)| = %.3e\n', max(abs(y)));
fprintf('eta(5) numerical = %.9f\n', x(end,1));
fprintf('eta(5) exact     = %.9f\n', eta0*exp(5));

figure;
plot(t, x(:,1), 'LineWidth', 1.5);
grid on;
xlabel('time');
ylabel('eta = x_1');
title('Zero Dynamics: eta_dot = eta');

figure;
plot(t, y, 'LineWidth', 1.5);
grid on;
xlabel('time');
ylabel('y');
title('Output Held at Zero');

% Simulink construction idea:
% 1. Add a State-Space block with A, B, C, D.
% 2. Use an internal state output or an equivalent integrator model.
% 3. Feed u = 24*x1 into the State-Space input.
% 4. Scope y and x1.  The output y remains zero while x1 grows.
      

Mathematica can represent the same state-space system, convert it to a transfer function, and solve the zero dynamics symbolically.

(* Chapter21_Lesson3.nb *)
ClearAll["Global`*"];

A = { {0, 1, 0}, {0, 0, 1}, {-6, -11, -6} };
B = { {0}, {0}, {1} };
Cmat = { {-1, 1, 0} };
Dmat = { {0} };

sys = StateSpaceModel[{A, B, Cmat, Dmat}];
tf = TransferFunctionModel[sys, s] // FullSimplify

zeros = s /. Solve[Numerator[tf[[1, 1, 1]]] == 0, s]
poles = s /. Solve[Denominator[tf[[1, 1, 1]]] == 0, s]

eta0 = 0.02;
solEta = DSolveValue[{eta'[t] == eta[t], eta[0] == eta0}, eta[t], t]
etaAt5 = solEta /. t -> 5

u[t_] := 24 solEta;
x[t_] := {solEta, solEta, solEta};
y[t_] := Cmat.x[t] // Simplify

Plot[Evaluate[solEta], {t, 0, 5},
     PlotLabel -> "Zero Dynamics eta_dot = eta",
     AxesLabel -> {"t", "eta"}]

Plot[Evaluate[y[t][[1]]], {t, 0, 5},
     PlotLabel -> "Output y(t) held at zero",
     AxesLabel -> {"t", "y"}]
      

Problem 1: Consider the SISO system \( \dot{x}=Ax+bu,\ y=cx \) with \( D=0 \) and relative degree \( r=1 \). Derive the input that maintains \( y(t)=0 \).

Solution: Since \( r=1 \), \( cb\ne 0 \). The output derivative is

\[ \dot{y}=c\dot{x}=cAx+cbu. \]

To keep \( y=0 \), impose \( \dot{y}=0 \). Therefore

\[ u=-\frac{cAx}{cb}. \]

Substitution gives the constrained internal motion

\[ \dot{x} = \left(A-b\frac{cA}{cb}\right)x, \qquad cx=0. \]

Problem 2: Show that if \( z \) is a transmission zero, then there exists an output-nulling exponential trajectory.

Solution: By definition, a finite transmission zero satisfies

\[ \begin{bmatrix} zI-A & -B\\ C & D \end{bmatrix} \begin{bmatrix} x_z\\ u_z \end{bmatrix} = \begin{bmatrix} 0\\ 0 \end{bmatrix}, \qquad x_z\ne 0. \]

Therefore

\[ (zI-A)x_z-Bu_z=0, \qquad Cx_z+Du_z=0. \]

Let \( x(t)=e^{zt}x_z \) and \( u(t)=e^{zt}u_z \). Then

\[ \dot{x}(t)=ze^{zt}x_z = e^{zt}(Ax_z+Bu_z) = Ax(t)+Bu(t), \]

and

\[ y(t)=Cx(t)+Du(t)=e^{zt}(Cx_z+Du_z)=0. \]

Thus every invariant zero corresponds to a zero-output internal exponential motion.

Problem 3: For \( G(s)=\frac{s+2}{(s+1)(s+3)} \), determine whether the continuous-time zero dynamics are stable.

Solution: The only finite zero is \( z=-2 \). Since

\[ \operatorname{Re}(-2)<0, \]

the zero dynamics are asymptotically stable. For a minimal SISO realization, the system is minimum phase.

Problem 4: For \( G(s)=\frac{s-4}{(s+2)(s+5)} \), explain why exact output-nulling motion is internally unstable.

Solution: The zero is \( z=+4 \). The output-nulling internal mode is proportional to \( e^{4t} \). Therefore, even if the measured output is zero, the internal state grows exponentially. The zero dynamics are unstable, and the system is non-minimum phase.

Problem 5: Explain why unstable zero dynamics place a fundamental limitation on inversion-based tracking.

Solution: Exact inversion attempts to choose the input so that the output follows a desired trajectory. If the system has unstable zero dynamics, then forcing the output to behave exactly can excite internal modes that grow without appearing immediately in the output. Hence stable causal exact inversion is impossible for non-minimum-phase systems. Practical controllers must accept slower tracking, undershoot, preview, approximate inversion, or other performance compromises.

Zero dynamics describe the internal motion compatible with \( y(t)=0 \). Transmission zeros are not merely numerator roots; they are internal output-nulling modes. In minimal SISO systems, the eigenvalues of the zero dynamics coincide with the finite transmission zeros. If all zeros lie in the open left half-plane, the system is minimum phase. If any zero lies in the right half-plane, the system has unstable zero dynamics and exact inversion or aggressive tracking becomes fundamentally limited.