In state-space design, a disturbance is not merely an unspecified input; it has a precise entry channel. We write the plant as

\[ \dot{x}=Ax+Bu+Ed,\qquad y=Cx+Du+Jd. \]

Here \( x\in\mathbb{R}^n \) is the state, \( u\in\mathbb{R}^m \) is the control input, \( d\in\mathbb{R}^q \) is the disturbance, and \( y\in\mathbb{R}^p \) is the controlled output. In many regulation problems \(D=0\) and \(J=0\), but the disturbance channel \(E\) is crucial. If the same physical actuator channel creates both \(Bu\) and \(Ed\), then cancellation may be easy. If \(E\) excites directions poorly affected by \(B\), exact rejection may be impossible with the available actuator.

flowchart TD
  R["reference r"] --> SUM["tracking error / integrator"]
  SUM --> CTRL["state-space controller"]
  X["measured state x"] --> CTRL
  D["disturbance d through E"] --> PLANT["plant: xdot = A x + B u + E d"]
  CTRL -->|"control u"| PLANT
  PLANT -->|"state x"| X
  PLANT -->|"controlled output y"| Y["output y"]
  Y --> SUM
        

A state-feedback law \(u=-Kx\) changes the internal dynamics but does not remove the disturbance channel. With \(A_K=A-BK\), the closed-loop model is

\[ \dot{x}=A_Kx+Ed, \qquad y=Cx. \]

The disturbance-to-output transfer matrix is therefore

\[ G_{yd,K}(s)=C(sI-A_K)^{-1}E. \]

Suppose \(A_K\) is asymptotically stable and the disturbance is constant: \(d(t)=d_0\). If no disturbance feedforward or internal model is used, the equilibrium is

\[ 0=A_Kx_\infty+Ed_0, \qquad x_\infty=-A_K^{-1}Ed_0. \]

The corresponding steady output is

\[ y_\infty=-CA_K^{-1}Ed_0. \]

Thus pole placement can reduce the transient response to disturbances, but it does not generally enforce \(y_\infty=0\). Exact rejection of a persistent disturbance requires either a measured disturbance compensation term, an internal model of the disturbance, or a structural decoupling condition.

Proof by final value theorem. For zero initial state,

\[ Y(s)=C(sI-A_K)^{-1}E\frac{d_0}{s}. \]

Since \(A_K\) is stable, the final value theorem gives

\[ \lim_{t\to\infty}y(t)=\lim_{s\to 0}sY(s) =C(0I-A_K)^{-1}Ed_0=-CA_K^{-1}Ed_0. \]

If \(d(t)\) is measured or accurately estimated by a separate instrument, one can use

\[ u=-Kx+F_dd. \]

For constant disturbances, exact steady-state rejection is formulated by the regulator equations. Find matrices \(\Pi\) and \(\Gamma\) such that

\[ 0=(A-BK)\Pi+B\Gamma+E, \qquad 0=C\Pi. \]

Then choosing \(F_d=\Gamma\) makes \(x_\infty=\Pi d_0\) and \(y_\infty=C\Pi d_0=0\). In block form, the constant disturbance regulator equation is

\[ \begin{bmatrix} A-BK & B \\ C & 0 \end{bmatrix} \begin{bmatrix} \Pi \\ \Gamma \end{bmatrix} =\begin{bmatrix} -E \\ 0 \end{bmatrix}. \]

A particularly transparent special case is \(\operatorname{im}E\subseteq\operatorname{im}B\). Then there exists \(F_d\) such that \(BF_d+E=0\), and the disturbance is canceled directly in the state equation:

\[ \dot{x}=(A-BK)x+(BF_d+E)d=(A-BK)x. \]

This cancellation is nonrobust: if the disturbance channel or measured disturbance is inaccurate, residual disturbance enters the loop.

If the constant disturbance is not measured, the standard state-space design is to add an integral state for the regulated output error. For regulation to zero, define

\[ \dot{\xi}=y=Cx. \]

The augmented system is

\[ \frac{d}{dt}\begin{bmatrix}x\\ \xi\end{bmatrix} =\underbrace{\begin{bmatrix}A&0\\ C&0\end{bmatrix} }_{A_a} \begin{bmatrix}x\\ \xi\end{bmatrix} +\underbrace{\begin{bmatrix}B\\0\end{bmatrix} }_{B_a}u +\underbrace{\begin{bmatrix}E\\0\end{bmatrix} }_{E_a}d. \]

Apply augmented state feedback \(u=-K_xx-K_i\xi\). The closed-loop matrix is

\[ A_{cl,a}=\begin{bmatrix}A-BK_x&-BK_i\\ C&0\end{bmatrix}. \]

If the pair \((A_a,B_a)\) is stabilizable and the gain is chosen so that all eigenvalues of \(A_{cl,a}\) satisfy \(\operatorname{Re}\lambda_i<0\), then any constant disturbance produces zero steady-state output.

Proof. At a closed-loop equilibrium under constant \(d_0\),

\[ 0=(A-BK_x)x_\infty-BK_i\xi_\infty+Ed_0, \qquad 0=Cx_\infty. \]

The second equation is exactly the desired regulation condition. Since the augmented closed-loop matrix is stable, the equilibrium is unique and globally attractive for the linear system. Hence \(\lim_{t\to\infty}y(t)=Cx_\infty=0\).

Integral disturbance rejection requires that the augmented integrator be controllable or at least stabilizable. For a SISO output with \(D=0\), the critical obstruction is a transmission zero at the origin. The augmented controllability matrix has full rank precisely when the Rosenbrock matrix at \(s=0\) has full row rank:

\[ \operatorname{rank}\begin{bmatrix}A&B\\ C&0\end{bmatrix}=n+p. \]

If this rank condition fails, the integrator introduces an uncontrollable mode at the origin, so no static augmented feedback can move all augmented poles into the stable half-plane. This is the state-space version of a familiar servo limitation: one cannot force exact DC rejection through an input-output channel whose zero blocks DC authority.

Lesson 3 introduced the internal model principle. Here we express the same idea in regulator-equation form. Suppose the disturbance is generated by

\[ \dot{w}=Sw, \qquad d=Fw. \]

A constant disturbance corresponds to \(S=0\). A sinusoidal disturbance of frequency \(\omega\) uses a two-state oscillator

\[ S=\begin{bmatrix}0&\omega\\-\omega&0\end{bmatrix}. \]

If \(w\) is available, steady-state rejection is characterized by

\[ \Pi S=A\Pi+B\Gamma+EF, \qquad 0=C\Pi. \]

If \(w\) is not available, a controller must contain a copy of the unstable or marginally stable modes of \(S\) that appear in the disturbance. This is why an integrator rejects constants and a resonant internal model rejects sinusoids.

Steady-state rejection is weaker than exact decoupling. Exact disturbance decoupling asks for a feedback \(u=Fx+v\) such that the transfer from \(d\) to the controlled output is identically zero:

\[ C(sI-A-BF)^{-1}E\equiv 0. \]

A geometric sufficient and necessary condition can be stated using a controlled-invariant subspace. There must exist a subspace \(\mathcal{V}\) satisfying

\[ \operatorname{im}E\subseteq\mathcal{V}\subseteq\ker C, \qquad A\mathcal{V}\subseteq\mathcal{V}+\operatorname{im}B. \]

The second inclusion means that the control input can keep the disturbance-excited state component inside an output-null subspace. If a feedback \(F\) is selected so that \((A+BF)\mathcal{V}\subseteq\mathcal{V}\), then a trajectory starting in \(\mathcal{V}\) remains in \(\mathcal{V}\), and since \(\mathcal{V}\subseteq\ker C\), its contribution to the controlled output is zero.

This result is powerful but structural. It says exact decoupling is not a matter of placing poles faster; it is a question of whether disturbance directions can be confined to output-invisible controlled-invariant directions.

The programming labs use a second-order plant with a constant force-like disturbance:

\[ A=\begin{bmatrix}0&1\\-2&-0.6\end{bmatrix},\qquad B=\begin{bmatrix}0\\1\end{bmatrix},\qquad E=\begin{bmatrix}0\\1\end{bmatrix},\qquad C=\begin{bmatrix}1&0\end{bmatrix}. \]

Since \(E=B\), measured disturbance feedforward can directly cancel the disturbance. For the unmeasured case, use the augmented matrices

\[ A_a=\begin{bmatrix}0&1&0\\-2&-0.6&0\\1&0&0\end{bmatrix}, \qquad B_a=\begin{bmatrix}0\\1\\0\end{bmatrix}. \]

Placing the augmented poles at \(-2,-2.5,-3\) gives

\[ K_a=\begin{bmatrix}16.5&6.9&15\end{bmatrix}, \qquad u=-16.5x_1-6.9x_2-15\xi. \]

The Python implementation uses numpy, scipy, and matplotlib. It compares state feedback alone, measured disturbance feedforward from the regulator equations, and integral action for an unmeasured constant disturbance.

Chapter27_Lesson4.py

"""
Chapter27_Lesson4.py
State-space design for disturbance rejection.

Benchmark plant:
    x_dot = A x + B u + E d
    y     = C x

The script compares:
1. state feedback only,
2. measured-disturbance feedforward from regulator equations,
3. unmeasured constant-disturbance rejection using integral action.

Required libraries:
    numpy, scipy, matplotlib
"""

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


A = np.array([[0.0, 1.0], [-2.0, -0.6]])
B = np.array([[0.0], [1.0]])
E = np.array([[0.0], [1.0]])
C = np.array([[1.0, 0.0]])

d0 = 0.5
x0 = np.array([0.4, 0.0])
t_span = (0.0, 10.0)
t_eval = np.linspace(t_span[0], t_span[1], 1000)


def simulate_state_feedback_only():
    """Closed loop u = -Kx. Constant disturbance generally leaves offset."""
    K = place_poles(A, B, [-2.0, -3.0]).gain_matrix

    def rhs(_t, x):
        u = (-K @ x.reshape(-1, 1)).item()
        return (A @ x.reshape(-1, 1) + B * u + E * d0).ravel()

    sol = solve_ivp(rhs, t_span, x0, t_eval=t_eval, rtol=1e-9, atol=1e-11)
    y = (C @ sol.y).ravel()
    return sol.t, y, K


def measured_disturbance_feedforward():
    """Solve regulator equations for known constant d."""
    K = place_poles(A, B, [-2.0, -3.0]).gain_matrix

    block = np.block([[A - B @ K, B], [C, np.zeros((1, 1))]])
    rhs_reg = np.vstack([-E, [[0.0]]])
    solution = np.linalg.solve(block, rhs_reg)
    Pi = solution[:2, :]
    Gamma = solution[2:, :]

    def rhs(_t, x):
        # u = -Kx + Gamma d, where Gamma is obtained from regulator equations.
        u = (-K @ x.reshape(-1, 1) + Gamma * d0).item()
        return (A @ x.reshape(-1, 1) + B * u + E * d0).ravel()

    sol = solve_ivp(rhs, t_span, x0, t_eval=t_eval, rtol=1e-9, atol=1e-11)
    y = (C @ sol.y).ravel()
    return sol.t, y, K, Pi, Gamma


def integral_disturbance_rejection():
    """Unmeasured constant disturbance rejection with xi_dot = y."""
    A_aug = np.block([[A, np.zeros((2, 1))], [C, np.zeros((1, 1))]])
    B_aug = np.vstack([B, [[0.0]]])
    E_aug = np.vstack([E, [[0.0]]])

    K_aug = place_poles(A_aug, B_aug, [-2.0, -2.5, -3.0]).gain_matrix
    K = K_aug[:, :2]
    Ki = K_aug[:, 2:]

    def rhs(_t, xa):
        x = xa[:2].reshape(-1, 1)
        xi = np.array([[xa[2]]])
        u = (-K @ x - Ki @ xi).item()
        dx = A @ x + B * u + E * d0
        dxi = C @ x
        return np.vstack([dx, dxi]).ravel()

    sol = solve_ivp(rhs, t_span, np.array([x0[0], x0[1], 0.0]),
                    t_eval=t_eval, rtol=1e-9, atol=1e-11)
    y = (C @ sol.y[:2, :]).ravel()
    return sol.t, y, K, Ki, A_aug, B_aug, E_aug


if __name__ == "__main__":
    t1, y1, K1 = simulate_state_feedback_only()
    t2, y2, K2, Pi, Gamma = measured_disturbance_feedforward()
    t3, y3, K3, Ki, A_aug, B_aug, E_aug = integral_disturbance_rejection()

    print("State feedback K:", K1)
    print("Measured-disturbance regulator Pi:\n", Pi)
    print("Measured-disturbance feedforward Gamma:", Gamma)
    print("Integral controller K:", K3)
    print("Integral gain Ki:", Ki)
    print("Final outputs:")
    print("  state feedback only        =", y1[-1])
    print("  measured feedforward       =", y2[-1])
    print("  integral disturbance reject=", y3[-1])

    plt.figure(figsize=(8, 4.8))
    plt.plot(t1, y1, label="state feedback only")
    plt.plot(t2, y2, label="measured feedforward")
    plt.plot(t3, y3, label="integral action")
    plt.axhline(0.0, linestyle="--", linewidth=1)
    plt.xlabel("time [s]")
    plt.ylabel("controlled output y")
    plt.title("Disturbance rejection in state space")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()
      

The C++ implementation uses Eigen and Ackermann pole placement on the augmented SISO pair.

Chapter27_Lesson4.cpp

/*
Chapter27_Lesson4.cpp
State-space disturbance rejection using integral action.

Dependency:
    Eigen 3
Compile example:
    g++ -std=c++17 Chapter27_Lesson4.cpp -I /path/to/eigen -O2 -o Chapter27_Lesson4
*/

#include <Eigen/Dense>
#include <fstream>
#include <iomanip>
#include <iostream>

using Matrix3 = Eigen::Matrix3d;
using Vector3 = Eigen::Vector3d;
using Row3 = Eigen::RowVector3d;

Matrix3 desiredPolynomialMatrix(const Matrix3& Aaug) {
    // Desired poles: -2, -2.5, -3
    // p(s) = s^3 + 7.5 s^2 + 18.5 s + 15
    return Aaug * Aaug * Aaug + 7.5 * Aaug * Aaug + 18.5 * Aaug + 15.0 * Matrix3::Identity();
}

Row3 ackermannGain(const Matrix3& Aaug, const Vector3& Baug) {
    Matrix3 ctrb;
    ctrb.col(0) = Baug;
    ctrb.col(1) = Aaug * Baug;
    ctrb.col(2) = Aaug * Aaug * Baug;

    Row3 last;
    last << 0.0, 0.0, 1.0;

    Matrix3 phiA = desiredPolynomialMatrix(Aaug);
    return last * ctrb.inverse() * phiA;
}

Vector3 dynamics(const Vector3& xa, const Row3& Kaug, double d0) {
    // Plant: x1_dot = x2, x2_dot = -2 x1 -0.6 x2 + u + d0, xi_dot = x1
    double x1 = xa(0);
    double x2 = xa(1);
    double xi = xa(2);
    double u = -(Kaug(0) * x1 + Kaug(1) * x2 + Kaug(2) * xi);

    Vector3 dx;
    dx << x2,
          -2.0 * x1 - 0.6 * x2 + u + d0,
          x1;
    return dx;
}

int main() {
    Matrix3 Aaug;
    Aaug << 0.0, 1.0, 0.0,
           -2.0, -0.6, 0.0,
            1.0, 0.0, 0.0;

    Vector3 Baug;
    Baug << 0.0, 1.0, 0.0;

    Row3 Kaug = ackermannGain(Aaug, Baug);
    std::cout << "Kaug = " << Kaug << std::endl;

    const double d0 = 0.5;
    const double h = 0.001;
    const double tf = 10.0;
    Vector3 xa;
    xa << 0.4, 0.0, 0.0;

    std::ofstream file("Chapter27_Lesson4_cpp_response.csv");
    file << "t,x1,x2,xi,u,y\n";

    for (int k = 0; k <= static_cast<int>(tf / h); ++k) {
        double t = k * h;
        double u = -(Kaug(0) * xa(0) + Kaug(1) * xa(1) + Kaug(2) * xa(2));
        double y = xa(0);
        file << std::setprecision(10) << t << "," << xa(0) << "," << xa(1) << ","
             << xa(2) << "," << u << "," << y << "\n";

        Vector3 k1 = dynamics(xa, Kaug, d0);
        Vector3 k2 = dynamics(xa + 0.5 * h * k1, Kaug, d0);
        Vector3 k3 = dynamics(xa + 0.5 * h * k2, Kaug, d0);
        Vector3 k4 = dynamics(xa + h * k3, Kaug, d0);
        xa += (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
    }

    std::cout << "Final y = " << xa(0) << std::endl;
    std::cout << "CSV written to Chapter27_Lesson4_cpp_response.csv" << std::endl;
    return 0;
}
      

The Java implementation is written from scratch. It builds the controllability matrix, applies Ackermann's formula, and simulates the augmented closed loop by fourth-order Runge-Kutta integration.

Chapter27_Lesson4.java

/*
Chapter27_Lesson4.java
State-space disturbance rejection using integral action.

Compile and run:
    javac Chapter27_Lesson4.java
    java Chapter27_Lesson4
*/

import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;
import java.util.Arrays;

public class Chapter27_Lesson4 {
    static double[][] multiply(double[][] A, double[][] B) {
        int m = A.length;
        int n = B[0].length;
        int p = B.length;
        double[][] C = new double[m][n];
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                for (int k = 0; k < p; k++) {
                    C[i][j] += A[i][k] * B[k][j];
                }
            }
        }
        return C;
    }

    static double[][] add(double[][] A, double[][] B) {
        double[][] C = new double[A.length][A[0].length];
        for (int i = 0; i < A.length; i++) {
            for (int j = 0; j < A[0].length; j++) {
                C[i][j] = A[i][j] + B[i][j];
            }
        }
        return C;
    }

    static double[][] scale(double[][] A, double alpha) {
        double[][] C = new double[A.length][A[0].length];
        for (int i = 0; i < A.length; i++) {
            for (int j = 0; j < A[0].length; j++) {
                C[i][j] = alpha * A[i][j];
            }
        }
        return C;
    }

    static double[][] identity(int n) {
        double[][] I = new double[n][n];
        for (int i = 0; i < n; i++) I[i][i] = 1.0;
        return I;
    }

    static double[][] inverse(double[][] A) {
        int n = A.length;
        double[][] aug = new double[n][2 * n];
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) aug[i][j] = A[i][j];
            aug[i][n + i] = 1.0;
        }

        for (int col = 0; col < n; col++) {
            int pivot = col;
            for (int row = col + 1; row < n; row++) {
                if (Math.abs(aug[row][col]) > Math.abs(aug[pivot][col])) {
                    pivot = row;
                }
            }
            double[] temp = aug[col];
            aug[col] = aug[pivot];
            aug[pivot] = temp;

            double div = aug[col][col];
            if (Math.abs(div) < 1e-12) throw new RuntimeException("Singular matrix");
            for (int j = 0; j < 2 * n; j++) aug[col][j] /= div;

            for (int row = 0; row < n; row++) {
                if (row == col) continue;
                double factor = aug[row][col];
                for (int j = 0; j < 2 * n; j++) {
                    aug[row][j] -= factor * aug[col][j];
                }
            }
        }

        double[][] inv = new double[n][n];
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) inv[i][j] = aug[i][n + j];
        }
        return inv;
    }

    static double[][] column(double a, double b, double c) {
        return new double[][]{ {a}, {b}, {c} };
    }

    static double[][] controllability(double[][] A, double[][] B) {
        double[][] AB = multiply(A, B);
        double[][] A2B = multiply(A, AB);
        return new double[][]{
                {B[0][0], AB[0][0], A2B[0][0]},
                {B[1][0], AB[1][0], A2B[1][0]},
                {B[2][0], AB[2][0], A2B[2][0]}
        };
    }

    static double[][] desiredPolynomialMatrix(double[][] A) {
        // Desired poles: -2, -2.5, -3
        // p(s) = s^3 + 7.5 s^2 + 18.5 s + 15
        double[][] A2 = multiply(A, A);
        double[][] A3 = multiply(A2, A);
        return add(add(add(A3, scale(A2, 7.5)), scale(A, 18.5)), scale(identity(3), 15.0));
    }

    static double[] ackermannGain(double[][] A, double[][] B) {
        double[][] ctrb = controllability(A, B);
        double[][] phiA = desiredPolynomialMatrix(A);
        double[][] last = new double[][]{ {0.0, 0.0, 1.0} };
        double[][] K = multiply(multiply(last, inverse(ctrb)), phiA);
        return K[0];
    }

    static double[] dynamics(double[] xa, double[] K, double d0) {
        double x1 = xa[0];
        double x2 = xa[1];
        double xi = xa[2];
        double u = -(K[0] * x1 + K[1] * x2 + K[2] * xi);
        return new double[]{
                x2,
                -2.0 * x1 - 0.6 * x2 + u + d0,
                x1
        };
    }

    static double[] addScaled(double[] x, double[] dx, double h) {
        return new double[]{x[0] + h * dx[0], x[1] + h * dx[1], x[2] + h * dx[2]};
    }

    public static void main(String[] args) throws IOException {
        double[][] Aaug = new double[][]{
                {0.0, 1.0, 0.0},
                {-2.0, -0.6, 0.0},
                {1.0, 0.0, 0.0}
        };
        double[][] Baug = column(0.0, 1.0, 0.0);
        double[] K = ackermannGain(Aaug, Baug);
        System.out.println("Kaug = " + Arrays.toString(K));

        double d0 = 0.5;
        double h = 0.001;
        double tf = 10.0;
        double[] xa = new double[]{0.4, 0.0, 0.0};

        try (PrintWriter out = new PrintWriter(new FileWriter("Chapter27_Lesson4_java_response.csv"))) {
            out.println("t,x1,x2,xi,u,y");
            int steps = (int) (tf / h);
            for (int k = 0; k <= steps; k++) {
                double t = k * h;
                double u = -(K[0] * xa[0] + K[1] * xa[1] + K[2] * xa[2]);
                out.printf("%.10f,%.10f,%.10f,%.10f,%.10f,%.10f%n", t, xa[0], xa[1], xa[2], u, xa[0]);

                double[] k1 = dynamics(xa, K, d0);
                double[] k2 = dynamics(addScaled(xa, k1, 0.5 * h), K, d0);
                double[] k3 = dynamics(addScaled(xa, k2, 0.5 * h), K, d0);
                double[] k4 = dynamics(addScaled(xa, k3, h), K, d0);
                for (int i = 0; i < 3; i++) {
                    xa[i] += (h / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
                }
            }
        }

        System.out.println("Final y = " + xa[0]);
        System.out.println("CSV written to Chapter27_Lesson4_java_response.csv");
    }
}
      

The MATLAB script uses place from the Control System Toolbox and optionally creates a simple Simulink model for the augmented closed-loop disturbance-to-output simulation.

Chapter27_Lesson4.m

% Chapter27_Lesson4.m
% State-space disturbance rejection using measured feedforward and integral action.
% Required toolbox for place(): Control System Toolbox.

clear; clc; close all;

A = [0 1; -2 -0.6];
B = [0; 1];
E = [0; 1];
C = [1 0];
d0 = 0.5;
x0 = [0.4; 0];

%% 1) State feedback only: u = -Kx
K = place(A, B, [-2 -3]);
f1 = @(t, x) A*x + B*(-K*x) + E*d0;
[t1, x1] = ode45(f1, [0 10], x0);
y1 = (C*x1')';

%% 2) Measured disturbance feedforward from regulator equations
M = [A - B*K, B; C, 0];
rhs = [-E; 0];
sol = M \ rhs;
Pi = sol(1:2);
Gamma = sol(3);

f2 = @(t, x) A*x + B*(-K*x + Gamma*d0) + E*d0;
[t2, x2] = ode45(f2, [0 10], x0);
y2 = (C*x2')';

%% 3) Unmeasured constant disturbance rejection with xi_dot = y
Aaug = [A, zeros(2,1); C, 0];
Baug = [B; 0];
Kaug = place(Aaug, Baug, [-2 -2.5 -3]);
Kx = Kaug(1:2);
Ki = Kaug(3);

f3 = @(t, xa) [A*xa(1:2) + B*(-Kx*xa(1:2) - Ki*xa(3)) + E*d0;
               C*xa(1:2)];
[t3, xa3] = ode45(f3, [0 10], [x0; 0]);
y3 = (C*xa3(:,1:2)')';

fprintf('K = [%g %g]\n', K(1), K(2));
fprintf('Measured disturbance Gamma = %g\n', Gamma);
fprintf('Kaug = [%g %g %g]\n', Kaug(1), Kaug(2), Kaug(3));
fprintf('Final outputs: feedback only=%g, feedforward=%g, integral=%g\n', y1(end), y2(end), y3(end));

figure;
plot(t1, y1, 'LineWidth', 1.2); hold on;
plot(t2, y2, 'LineWidth', 1.2);
plot(t3, y3, 'LineWidth', 1.2);
yline(0, '--'); grid on;
xlabel('time [s]'); ylabel('controlled output y');
title('State-space disturbance rejection');
legend('state feedback only', 'measured feedforward', 'integral action', 'Location', 'best');

%% Optional Simulink model: augmented closed-loop from disturbance to y
% The model is: d0 -> State-Space(Acl, Bd, Cy, Dy) -> Scope.
if exist('simulink', 'file') == 2
    model = 'Chapter27_Lesson4_Simulink_Model';
    if bdIsLoaded(model)
        close_system(model, 0);
    end
    new_system(model);
    open_system(model);

    Acl = [A - B*Kx, -B*Ki; C, 0];
    Bd = [E; 0];
    Cy = [C, 0];
    Dy = 0;

    add_block('simulink/Sources/Constant', [model '/Constant disturbance'], ...
        'Value', num2str(d0), 'Position', [40 80 130 110]);
    add_block('simulink/Continuous/State-Space', [model '/Closed-loop augmented plant'], ...
        'A', mat2str(Acl), 'B', mat2str(Bd), 'C', mat2str(Cy), 'D', mat2str(Dy), ...
        'Position', [190 60 370 130]);
    add_block('simulink/Sinks/Scope', [model '/Output y'], ...
        'Position', [430 72 500 118]);
    add_line(model, 'Constant disturbance/1', 'Closed-loop augmented plant/1');
    add_line(model, 'Closed-loop augmented plant/1', 'Output y/1');
    set_param(model, 'StopTime', '10');
    save_system(model);
    fprintf('Simulink model saved as %s.slx\n', model);
end
      

The Mathematica implementation computes the augmented gain by Ackermann's formula and solves the closed-loop differential equations using NDSolve.

Chapter27_Lesson4.nb

Notebook[{
Cell["Chapter27_Lesson4.nb", "Title"],
Cell["State-space disturbance rejection with integral action", "Text"],
Cell[BoxData["(* Chapter27_Lesson4.nb *)\n(* State-space disturbance rejection with integral action. *)\n\nClearAll[\"Global`*\"];\n\nA = { {0, 1}, {-2, -0.6} };\nB = { {0}, {1} };\nEchan = { {0}, {1} };\nCmat = { {1, 0} };\nd0 = 0.5;\n\nAaug = ArrayFlatten[{ {A, ConstantArray[0, {2, 1}]}, {Cmat, { {0} } } }];\nBaug = Join[B, { {0} }];\n\nCtrb = {Baug[[All, 1]], Aaug.Baug[[All, 1]], MatrixPower[Aaug, 2].Baug[[All, 1]]} // Transpose;\nMatrixRank[Ctrb]\n\n(* Desired poles: -2, -2.5, -3, so p(s)=s^3+7.5 s^2+18.5 s+15. *)\nphiA = MatrixPower[Aaug, 3] + 7.5 MatrixPower[Aaug, 2] + 18.5 Aaug + 15 IdentityMatrix[3];\nKaug = { {0, 0, 1} }.Inverse[Ctrb].phiA;\nKaug // MatrixForm\n\nsol = NDSolve[\n   {\n    x1'[t] == x2[t],\n    x2'[t] == -2 x1[t] - 0.6 x2[t] - Kaug[[1, 1]] x1[t] - Kaug[[1, 2]] x2[t] - Kaug[[1, 3]] xi[t] + d0,\n    xi'[t] == x1[t],\n    x1[0] == 0.4,\n    x2[0] == 0,\n    xi[0] == 0\n    },\n   {x1, x2, xi}, {t, 0, 10}\n   ];\n\nPlot[Evaluate[x1[t] /. sol], {t, 0, 10},\n PlotLabel -> \"Controlled output y(t) under constant disturbance\",\n AxesLabel -> {\"t\", \"y\"}, PlotRange -> All]\n"], "Input"]
}]
      

Problem 1 (Offset under state feedback): Consider \(\dot{x}=(A-BK)x+Ed_0\), \(y=Cx\), where \(A-BK\) is stable and \(d_0\) is constant. Derive the steady-state output.

Solution: At equilibrium,

\[ 0=(A-BK)x_\infty+Ed_0. \]

Therefore

\[ x_\infty=-(A-BK)^{-1}Ed_0, \qquad y_\infty=-C(A-BK)^{-1}Ed_0. \]

Unless \(C(A-BK)^{-1}E=0\), stabilization alone does not produce exact disturbance rejection.

Problem 2 (Measured disturbance regulator equations): Let \(u=-Kx+\Gamma d\). Show that constant disturbance rejection is achieved if matrices \(\Pi\) and \(\Gamma\) satisfy

\[ 0=(A-BK)\Pi+B\Gamma+E, \qquad C\Pi=0. \]

Solution: For constant \(d_0\), set \(x_\infty=\Pi d_0\). Substitution into the equilibrium equation gives

\[ 0=(A-BK)\Pi d_0+B\Gamma d_0+Ed_0. \]

Since this must hold for any constant \(d_0\), the first regulator equation follows. The output at equilibrium is \(y_\infty=C\Pi d_0=0\), which follows from the second regulator equation.

Problem 3 (Augmented pole placement for the benchmark): For the benchmark plant in Section 9, form the augmented pair \((A_a,B_a)\) and compute an augmented feedback gain that assigns poles \(-2,-2.5,-3\).

Solution: The augmented matrices are

\[ A_a=\begin{bmatrix}0&1&0\\-2&-0.6&0\\1&0&0\end{bmatrix}, \qquad B_a=\begin{bmatrix}0\\1\\0\end{bmatrix}. \]

The controllability matrix is

\[ \mathcal{C}_a=\begin{bmatrix}B_a&A_aB_a&A_a^2B_a\end{bmatrix} =\begin{bmatrix}0&1&-0.6\\1&-0.6&-1.64\\0&0&1\end{bmatrix}, \qquad \det(\mathcal{C}_a)=-1. \]

The desired characteristic polynomial is \((s+2)(s+2.5)(s+3)=s^3+7.5s^2+18.5s+15\). By Ackermann's formula,

\[ K_a=\begin{bmatrix}0&0&1\end{bmatrix}\mathcal{C}_a^{-1} \left(A_a^3+7.5A_a^2+18.5A_a+15I\right) =\begin{bmatrix}16.5&6.9&15\end{bmatrix}. \]

Problem 4 (Transmission zero obstruction): Explain why a zero at the origin prevents integral action from rejecting a constant output disturbance through the chosen input-output channel.

Solution: A zero at the origin means that the Rosenbrock matrix

\[ \begin{bmatrix}A&B\\C&0\end{bmatrix} \]

loses rank. This rank loss implies that the augmented pair \((A_a,B_a)\) has an uncontrollable mode at the integrator eigenvalue \(0\). Since uncontrollable modes cannot be shifted by feedback, the augmented closed-loop cannot be made asymptotically stable. Therefore exact constant-disturbance rejection by integral action fails.

Problem 5 (Exact decoupling condition): Suppose \(D=J=0\). State a geometric condition for exact disturbance decoupling from \(d\) to \(y=Cx\) by state feedback.

Solution: There must exist a controlled-invariant subspace \(\mathcal{V}\) such that

\[ \operatorname{im}E\subseteq\mathcal{V}\subseteq\ker C, \qquad A\mathcal{V}\subseteq\mathcal{V}+\operatorname{im}B. \]

The first inclusion places disturbance directions inside an output-invisible subspace. The second says the control can make this subspace invariant. Therefore the disturbance-excited component remains invisible at the controlled output.

Disturbance rejection in state space begins with identifying the disturbance channel \(E\). Stabilizing state feedback generally reduces transients but leaves steady-state offset under persistent disturbances. Measured disturbances can be rejected with regulator-equation feedforward. Unmeasured constant disturbances require an internal model, implemented here as integral action. Exact disturbance decoupling is stronger and depends on structural subspace conditions, not merely on fast pole placement.