Consider the continuous-time LTI state-space model introduced in earlier chapters: state \( x(t)\in\mathbb{R}^n \), input \( u(t)\in\mathbb{R}^m \), output \( y(t)\in\mathbb{R}^p \), and constant matrices \( A\in\mathbb{R}^{n\times n} \), \( B\in\mathbb{R}^{n\times m} \), \( C\in\mathbb{R}^{p\times n} \), \( D\in\mathbb{R}^{p\times m} \).

\[ \dot{x}(t) = A x(t) + B u(t), \qquad y(t) = C x(t) + D u(t), \qquad x(0)=x_0. \]

From Lesson 2, the state solution (variation of constants) is:

\[ x(t) = e^{At}x_0 + \int_{0}^{t} e^{A(t-\tau)}\,B\,u(\tau)\,d\tau. \]

flowchart LR
  X0["initial state x0"] --> PHI["state transition Phi(t)=exp(A t)"]
  U["input u(t)"] --> BBLK["input map B"]
  BBLK --> INT["dynamic path via integral"]
  PHI --> SUMX["x(t)"]
  INT --> SUMX
  SUMX --> CBLK["output map C"]
  U --> DBLK["direct path D"]
  CBLK --> SUMY["y(t)"]
  DBLK --> SUMY
        

The diagram expresses the key idea: the output inherits both a dynamic dependence on the input via the internal state, and a possible instantaneous dependence through \( D \).

Substitute the known formula for \( x(t) \) into \( y(t)=Cx(t)+Du(t) \). This yields the fundamental output decomposition:

\[ \begin{aligned} y(t) &= C\left(e^{At}x_0 + \int_{0}^{t} e^{A(t-\tau)}\,B\,u(\tau)\,d\tau\right) + D u(t) \\ &= \underbrace{C e^{At}x_0}_{\text{initial-condition output} } + \underbrace{\int_{0}^{t} C e^{A(t-\tau)}B\,u(\tau)\,d\tau}_{\text{dynamic forced output} } + \underbrace{D u(t)}_{\text{direct feedthrough} }. \end{aligned} \]

Define the \( p\times m \) matrix-valued function \( h(t) \) for \( t \ge 0 \): \( h(t) \equiv C e^{At}B \). Then the forced dynamic part becomes a convolution.

\[ \int_{0}^{t} C e^{A(t-\tau)}B\,u(\tau)\,d\tau = \int_{0}^{t} h(t-\tau)\,u(\tau)\,d\tau \equiv (h * u)(t). \]

Therefore, the output can be written compactly as:

\[ y(t) = C e^{At}x_0 + (h*u)(t) + D u(t), \qquad h(t)=C e^{At}B,\; t \ge 0. \]

This section formalizes properties that are frequently used later in modern control (especially when connecting state-space to input-output representations).

3.1 Convolution Representation for Zero Initial Condition

Proposition 1 (forced output as convolution). If \( x_0 = 0 \) and \( u(\cdot) \) is piecewise continuous on \( [0,T] \), then for \( t\in[0,T] \):

\[ y(t) = \int_{0}^{t} C e^{A(t-\tau)}B\,u(\tau)\,d\tau + D u(t). \]

Proof. With \( x_0=0 \), the state solution is

\[ x(t)=\int_{0}^{t} e^{A(t-\tau)}B\,u(\tau)\,d\tau. \]

Multiply by \( C \) and add \( D u(t) \):

\[ y(t)=C\int_{0}^{t} e^{A(t-\tau)}B\,u(\tau)\,d\tau + D u(t) = \int_{0}^{t} C e^{A(t-\tau)}B\,u(\tau)\,d\tau + D u(t). \]

The interchange of multiplication by \( C \) and integration is valid because the integrand is piecewise continuous and matrix multiplication is linear. ■

3.2 Impulse Response and the Direct Term

In classical linear systems, the impulse response is the output produced by \( u(t)=\delta(t) \) with \( x_0=0 \). In state-space, there are two parts: the dynamic impulse response \( C e^{At}B \) for \( t \gt 0 \), and the instantaneous component due to \( D \) at \( t=0 \). To represent both in a single distributional object, define:

\[ h_{\mathrm{io} }(t) \equiv C e^{At}B + D\,\delta(t). \]

Then, for \( x_0=0 \), the complete input-output relation may be written (in the sense of distributions):

\[ y(t) = (h_{\mathrm{io} } * u)(t). \]

Interpretation. The term \( D\,\delta(t) \) encodes the fact that \( D u(t) \) is “instantaneous”: discontinuities in \( u \) can appear immediately in \( y \).

3.3 Continuity of x(t) and Jump Conditions for y(t)

Proposition 2 (state continuity). If \( u(\cdot) \) is piecewise continuous on \( [0,T] \), then \( x(t) \) given by \( x(t)=e^{At}x_0+\int_{0}^{t}e^{A(t-\tau)}B u(\tau)d\tau \) is continuous on \( [0,T] \).

Proof. The map \( t \mapsto e^{At}x_0 \) is continuous. For the integral term, define \( f(t,\tau)=e^{A(t-\tau)}B u(\tau) \). For fixed \( \tau \), the function \( t \mapsto e^{A(t-\tau)} \) is continuous, hence \( f(t,\tau) \) is continuous in \( t \) for each \( \tau \) where \( u(\tau) \) is finite. Since \( u \) is piecewise continuous, the integral is well-defined and depends continuously on the upper limit \( t \) by standard properties of parameter-dependent integrals. Thus, the sum is continuous. ■

Proposition 3 (output jumps are governed by D). Suppose \( x(t) \) is continuous and the input has a jump at \( t=t_0 \): \( \Delta u \equiv u(t_0^+) - u(t_0^-) \). Then the output jump satisfies

\[ \Delta y \equiv y(t_0^+) - y(t_0^-) = D\,\Delta u. \]

Proof. Since \( x(t) \) is continuous, \( Cx(t_0^+)-Cx(t_0^-)=0 \). Therefore,

\[ \Delta y = \big(Cx(t_0^+)+D u(t_0^+)\big) - \big(Cx(t_0^-)+D u(t_0^-)\big) = D\big(u(t_0^+)-u(t_0^-)\big) = D\,\Delta u. \]

■

When \( u(t)=u_0 \) is constant for \( t \ge 0 \), the forced integral becomes:

\[ \int_{0}^{t} e^{A(t-\tau)}B u_0\,d\tau = \left(\int_{0}^{t} e^{A s}\,ds\right) B u_0 \quad \text{(substitution } s=t-\tau \text{)}. \]

If \( A \) is invertible, the integral admits a closed form:

\[ \int_{0}^{t} e^{A s}\,ds = A^{-1}\left(e^{At}-I\right). \]

Proof. Differentiate the right-hand side:

\[ \frac{d}{dt}\left(A^{-1}(e^{At}-I)\right) = A^{-1}\left(A e^{At}\right) = e^{At}, \qquad \left. A^{-1}(e^{At}-I)\right|_{t=0}=0. \]

Hence it is the unique antiderivative satisfying the correct initial value. ■

Substituting into \( x(t) \) and then into \( y(t) \) gives:

\[ \boxed{ \begin{aligned} x(t) &= e^{At}x_0 + A^{-1}(e^{At}-I)B u_0, \\ y(t) &= C e^{At}x_0 + C A^{-1}(e^{At}-I)B u_0 + D u_0. \end{aligned} } \]

For a step input \( u(t)=u_0\,\mathbf{1}(t) \), the same formula applies for \( t \ge 0 \). If \( D\neq 0 \), then \( y(0^+)-y(0^-)=D u_0 \) (assuming \( u(0^-)=0 \)).

In practice, output computation uses either (i) analytic formulas via \( e^{At} \) and integrals, or (ii) numerical integration of the state equation and then \( y(t)=Cx(t)+Du(t) \). The following workflow clarifies typical engineering choices.

flowchart TD
  IN["Given A,B,C,D; x0; u(t)"] --> CH["Choose approach"]
  CH --> A1["Analytic: compute Phi(t)=exp(A t)"]
  A1 --> A2["Compute x(t) via closed form or integral"]
  A2 --> OUT1["y(t)=C x(t)+D u(t)"]
  CH --> N1["Numerical: integrate xdot=A x + B u"]
  N1 --> N2["Compute y(t)=C x(t)+D u(t) pointwise"]
  OUT1 --> VAL["Validate (compare methods, step size, conditioning)"]
  N2 --> VAL
        

Numerical note. When \( A \) is ill-conditioned or large-scale, computing \( e^{At} \) explicitly for many time points may be expensive; direct ODE integration can be preferable. Conversely, for repeated evaluation (e.g., many candidate inputs), precomputing \( e^{At} \) or using structure-exploiting methods can be advantageous.

Simulink note. If \( D \neq 0 \), a direct feedthrough can create algebraic loops in block diagrams, depending on how feedback is wired. This is a modeling issue (not a mathematical inconsistency): it means the output depends on the current input instantaneously.

The script below computes \( y(t) \) for a step input in three ways: (i) closed-form constant-input response using \( e^{At} \), (ii) scipy.signal.lsim (when available), (iii) a from-scratch forward Euler integrator (for learning, not best accuracy).

File: Chapter7_Lesson4.py

# Chapter7_Lesson4.py
"""
Modern Control — Chapter 7, Lesson 4
Output Response Using x(t) and y(t) = Cx(t) + Du(t)

This script demonstrates:
1) Analytic output response for constant (step) input using expm(A t)
2) Numeric simulation via SciPy (lsim) when available
3) "From scratch" forward-Euler simulation of xdot = A x + B u, y = C x + D u

Author: (student)
"""

import numpy as np

try:
    from scipy.linalg import expm
    from scipy import signal
    SCIPY_AVAILABLE = True
except Exception:
    SCIPY_AVAILABLE = False


def is_invertible(A: np.ndarray, tol: float = 1e-12) -> bool:
    return np.linalg.cond(A) < 1.0 / tol


def analytic_response_constant_u(A, B, C, D, x0, u0, t_grid):
    """
    For u(t)=u0 (constant), and invertible A:
        x(t) = e^{At} x0 + A^{-1}(e^{At} - I) B u0
        y(t) = C x(t) + D u0
    If A is not invertible or SciPy is unavailable, we fall back to numerical quadrature.
    """
    n = A.shape[0]
    I = np.eye(n)
    y = np.zeros((len(t_grid), C.shape[0]))
    x = np.zeros((len(t_grid), n))

    if SCIPY_AVAILABLE and is_invertible(A):
        Ainv = np.linalg.inv(A)
        for k, t in enumerate(t_grid):
            Phi = expm(A * t)
            xk = Phi @ x0 + (Ainv @ (Phi - I) @ B) @ u0
            yk = (C @ xk) + (D @ u0)
            x[k, :] = xk.ravel()
            y[k, :] = yk.ravel()
        return x, y

    # Fallback: numerical integration of convolution integral using trapezoidal rule
    # x(t) = e^{At} x0 + ∫_0^t e^{A(t-τ)} B u0 dτ
    if not SCIPY_AVAILABLE:
        raise RuntimeError("SciPy not available; cannot compute expm for fallback integration.")

    for k, t in enumerate(t_grid):
        Phi_t = expm(A * t)
        # integrate over τ in [0,t]
        m = max(2, int(400 * t / (t_grid[-1] + 1e-12)))
        tau = np.linspace(0.0, t, m)
        dtau = tau[1] - tau[0]
        integ = np.zeros((n, 1))
        for j, tj in enumerate(tau):
            weight = 0.5 if (j == 0 or j == m - 1) else 1.0
            integ += weight * (expm(A * (t - tj)) @ B) @ u0
        integ *= dtau
        xk = Phi_t @ x0 + integ
        yk = (C @ xk) + (D @ u0)
        x[k, :] = xk.ravel()
        y[k, :] = yk.ravel()

    return x, y


def euler_simulation(A, B, C, D, x0, u_fun, t_grid):
    """
    Forward-Euler simulation (educational; step size affects accuracy):
        x_{k+1} = x_k + h (A x_k + B u(t_k))
        y_k = C x_k + D u(t_k)
    """
    n = A.shape[0]
    p = C.shape[0]
    x = np.zeros((len(t_grid), n))
    y = np.zeros((len(t_grid), p))

    x[0, :] = x0.ravel()
    for k in range(len(t_grid) - 1):
        h = t_grid[k + 1] - t_grid[k]
        uk = u_fun(t_grid[k])
        y[k, :] = (C @ x[k, :].reshape(-1, 1) + D @ uk).ravel()
        xdot = (A @ x[k, :].reshape(-1, 1)) + (B @ uk)
        x[k + 1, :] = (x[k, :].reshape(-1, 1) + h * xdot).ravel()
    # last output
    y[-1, :] = (C @ x[-1, :].reshape(-1, 1) + D @ u_fun(t_grid[-1])).ravel()
    return x, y


def main():
    # Example: 2-state SISO system
    A = np.array([[0.0, 1.0],
                  [-2.0, -3.0]])
    B = np.array([[0.0],
                  [1.0]])
    C = np.array([[1.0, 0.0]])
    D = np.array([[0.25]])  # include direct feedthrough to show output jump at step changes

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

    # Step input u(t) = 1 for t >= 0
    u0 = np.array([[1.0]])

    t = np.linspace(0.0, 6.0, 601)

    # Analytic response for constant input
    if SCIPY_AVAILABLE:
        x_a, y_a = analytic_response_constant_u(A, B, C, D, x0, u0, t)
        print("Analytic (constant-input) y(t) at a few times:")
        for ti in [0.0, 0.5, 1.0, 2.0, 6.0]:
            k = int(round(ti / (t[1] - t[0])))
            print(f"t={t[k]:.2f}, y={y_a[k,0]:.6f}")
    else:
        print("SciPy not available: skipping analytic expm-based response.")

    # Numeric simulation with SciPy lsim (if available)
    if SCIPY_AVAILABLE:
        sys = signal.StateSpace(A, B, C, D)
        u = np.ones_like(t)  # step
        tout, y_lsim, x_lsim = signal.lsim(sys, U=u, T=t, X0=x0.ravel())
        print("\nSciPy lsim y(t) at a few times:")
        for ti in [0.0, 0.5, 1.0, 2.0, 6.0]:
            k = int(round(ti / (t[1] - t[0])))
            print(f"t={tout[k]:.2f}, y={y_lsim[k]:.6f}")
    else:
        print("SciPy not available: skipping lsim.")

    # From-scratch Euler
    def u_fun(_t):
        return u0

    x_e, y_e = euler_simulation(A, B, C, D, x0, u_fun, t)
    print("\nEuler (from scratch) y(t) at a few times:")
    for ti in [0.0, 0.5, 1.0, 2.0, 6.0]:
        k = int(round(ti / (t[1] - t[0])))
        print(f"t={t[k]:.2f}, y={y_e[k,0]:.6f}")

    # Optional: save to CSV for plotting elsewhere
    data = np.column_stack([t, y_e[:, 0]])
    np.savetxt("Chapter7_Lesson4_output_euler.csv", data, delimiter=",", header="t,y", comments="")
    print("\nSaved Euler output to Chapter7_Lesson4_output_euler.csv")


if __name__ == "__main__":
    main()

Recommended Python packages for this lesson’s workflow: \( \) numpy, scipy (matrix exponential, simulation), and (optionally) control (Control Systems Library) for higher-level state-space utilities.

The following C++ example uses Eigen’s matrix exponential (unsupported module) to compute \( \Phi(t)=e^{At} \), then evaluates the constant-input closed-form response and compares it with a forward Euler simulation.

File: Chapter7_Lesson4.cpp

// Chapter7_Lesson4.cpp
/*
Modern Control — Chapter 7, Lesson 4
Output Response Using x(t) and y(t) = Cx(t) + Du(t)

This C++ example uses Eigen to:
1) Compute Phi(t)=exp(A t) using Eigen's MatrixFunctions (unsupported module)
2) Compute analytic response for constant input u(t)=u0 when A is invertible:
       x(t) = exp(A t) x0 + A^{-1}(exp(A t) - I) B u0
       y(t) = C x(t) + D u0
3) Provide an educational forward-Euler simulator as a fallback/validation

Build (example):
  g++ -O2 -std=c++17 Chapter7_Lesson4.cpp -I /path/to/eigen -o Chapter7_Lesson4

Note: Eigen's matrix exponential is in the unsupported module.
*/

#include <iostream>
#include <vector>
#include <cmath>

#include <Eigen/Dense>
#include <unsupported/Eigen/MatrixFunctions>

using Eigen::MatrixXd;
using Eigen::VectorXd;

static bool isInvertible(const MatrixXd& A, double tol = 1e-12) {
    Eigen::FullPivLU<MatrixXd> lu(A);
    return lu.isInvertible() && (1.0 / A.fullPivLu().rcond() < 1.0 / tol);
}

int main() {
    // Example: 2-state SISO system
    MatrixXd A(2,2);
    A << 0.0,  1.0,
        -2.0, -3.0;

    MatrixXd B(2,1);
    B << 0.0,
         1.0;

    MatrixXd C(1,2);
    C << 1.0, 0.0;

    MatrixXd D(1,1);
    D << 0.25; // direct feedthrough to show output jump at a step

    VectorXd x0(2);
    x0 << 1.0, 0.0;

    VectorXd u0(1);
    u0 << 1.0;

    const double t_end = 6.0;
    const double h = 0.01;
    const int N = static_cast<int>(std::round(t_end / h)) + 1;

    std::vector<double> t(N);
    for (int k = 0; k < N; ++k) t[k] = k * h;

    // Analytic constant-input response (requires invertible A)
    if (!isInvertible(A)) {
        std::cerr << "A is not invertible; analytic constant-input formula not used.\n";
    } else {
        MatrixXd Ainv = A.inverse();
        std::cout << "Analytic y(t) (constant input) at a few times:\n";
        for (double ti : {0.0, 0.5, 1.0, 2.0, 6.0}) {
            MatrixXd Phi = (A * ti).exp();                    // expm
            VectorXd xt = Phi * x0 + (Ainv * (Phi - MatrixXd::Identity(2,2)) * B) * u0(0);
            VectorXd yt = C * xt + D * u0;
            std::cout << "t=" << ti << ", y=" << yt(0) << "\n";
        }
    }

    // From-scratch Euler simulation
    std::vector<VectorXd> x(N, VectorXd::Zero(2));
    std::vector<VectorXd> y(N, VectorXd::Zero(1));

    x[0] = x0;
    for (int k = 0; k < N-1; ++k) {
        VectorXd uk = u0;
        y[k] = C * x[k] + D * uk;

        VectorXd xdot = A * x[k] + B * uk(0);
        x[k+1] = x[k] + h * xdot;
    }
    y[N-1] = C * x[N-1] + D * u0;

    std::cout << "\nEuler y(t) at a few times:\n";
    for (double ti : {0.0, 0.5, 1.0, 2.0, 6.0}) {
        int k = static_cast<int>(std::round(ti / h));
        std::cout << "t=" << t << "\n";
    }

    return 0;
}

Recommended C++ libraries for modern-control computations: Eigen (linear algebra), and optionally SLICOT bindings (advanced state-space routines) in specialized environments.

Java does not ship with matrix exponentials; EJML provides eigen-decompositions and basic matrix operations. The code below computes \( e^{At} \) using diagonalization when possible, with a series fallback.

File: Chapter7_Lesson4.java

// Chapter7_Lesson4.java
/*
Modern Control — Chapter 7, Lesson 4
Output Response Using x(t) and y(t) = Cx(t) + Du(t)

This Java example uses EJML (Efficient Java Matrix Library) to compute
an approximation of exp(A t) via eigen-decomposition (diagonalization):
    A = V Λ V^{-1}  (when A is diagonalizable)
    exp(A t) = V exp(Λ t) V^{-1}

Then, for constant input u(t)=u0 and invertible A:
    x(t) = exp(A t) x0 + A^{-1}(exp(A t) - I) B u0
    y(t) = C x(t) + D u0

If eigen-decomposition fails (e.g., defective A), the code falls back to
a truncated series exp(A t) ≈ Σ_{k=0}^K (A t)^k / k! (educational).

Dependency (Gradle/Maven):
  org.ejml:ejml-all:0.43+

Run:
  javac -cp ejml-all-0.43.jar Chapter7_Lesson4.java
  java  -cp .:ejml-all-0.43.jar Chapter7_Lesson4
*/

import org.ejml.data.DMatrixRMaj;
import org.ejml.dense.row.CommonOps_DDRM;
import org.ejml.interfaces.decomposition.EigenDecomposition_F64;
import org.ejml.dense.row.factory.DecompositionFactory_DDRM;

public class Chapter7_Lesson4 {

    static DMatrixRMaj expmSeries(DMatrixRMaj A, double t, int K) {
        int n = A.numRows;
        DMatrixRMaj At = A.copy();
        CommonOps_DDRM.scale(t, At);

        DMatrixRMaj X = CommonOps_DDRM.identity(n); // term 0
        DMatrixRMaj term = CommonOps_DDRM.identity(n);

        for (int k = 1; k <= K; k++) {
            DMatrixRMaj tmp = new DMatrixRMaj(n, n);
            CommonOps_DDRM.mult(term, At, tmp);
            term = tmp;
            CommonOps_DDRM.scale(1.0 / k, term);     // divide by k iteratively => /k!
            CommonOps_DDRM.addEquals(X, term);
        }
        return X;
    }

    static DMatrixRMaj expmEigen(DMatrixRMaj A, double t) {
        int n = A.numRows;
        EigenDecomposition_F64<DMatrixRMaj> eig = DecompositionFactory_DDRM.eig(n, true);
        if (!eig.decompose(A)) {
            return null;
        }

        DMatrixRMaj V = new DMatrixRMaj(n, n);
        DMatrixRMaj Vinv = new DMatrixRMaj(n, n);
        DMatrixRMaj LambdaExp = new DMatrixRMaj(n, n);

        // Build V and exp(Lambda t). This handles only real eigenpairs cleanly.
        for (int i = 0; i < n; i++) {
            if (eig.getEigenvalue(i).isComplex()) {
                return null;
            }
            double lam = eig.getEigenvalue(i).getReal();
            DMatrixRMaj vi = eig.getEigenVector(i);
            if (vi == null) return null;
            for (int r = 0; r < n; r++) {
                V.set(r, i, vi.get(r, 0));
            }
            LambdaExp.set(i, i, Math.exp(lam * t));
        }

        // Vinv = inv(V)
        if (!CommonOps_DDRM.invert(V, Vinv)) {
            return null;
        }

        DMatrixRMaj tmp = new DMatrixRMaj(n, n);
        DMatrixRMaj Phi = new DMatrixRMaj(n, n);
        CommonOps_DDRM.mult(V, LambdaExp, tmp);
        CommonOps_DDRM.mult(tmp, Vinv, Phi);
        return Phi;
    }

    public static void main(String[] args) {
        // Example: 2-state SISO system
        DMatrixRMaj A = new DMatrixRMaj(new double[][]{
                {0.0, 1.0},
                {-2.0, -3.0}
        });
        DMatrixRMaj B = new DMatrixRMaj(new double[][]{
                {0.0},
                {1.0}
        });
        DMatrixRMaj C = new DMatrixRMaj(new double[][]{
                {1.0, 0.0}
        });
        DMatrixRMaj D = new DMatrixRMaj(new double[][]{
                {0.25}
        });

        DMatrixRMaj x0 = new DMatrixRMaj(new double[][]{
                {1.0},
                {0.0}
        });
        DMatrixRMaj u0 = new DMatrixRMaj(new double[][]{
                {1.0}
        });

        // Invert A (for constant-input formula)
        DMatrixRMaj Ainv = A.copy();
        if (!CommonOps_DDRM.invert(Ainv)) {
            System.out.println("A is not invertible; cannot use constant-input closed form.");
            return;
        }

        double[] times = new double[]{0.0, 0.5, 1.0, 2.0, 6.0};

        System.out.println("y(t) for constant u(t)=1:");
        for (double t : times) {
            DMatrixRMaj Phi = expmEigen(A, t);
            if (Phi == null) {
                Phi = expmSeries(A, t, 30);
            }

            // x(t) = Phi x0 + A^{-1}(Phi - I) B u0
            DMatrixRMaj I = CommonOps_DDRM.identity(2);
            DMatrixRMaj PhiMinusI = new DMatrixRMaj(2, 2);
            CommonOps_DDRM.subtract(Phi, I, PhiMinusI);

            DMatrixRMaj tmp1 = new DMatrixRMaj(2, 1);
            CommonOps_DDRM.mult(Phi, x0, tmp1);

            DMatrixRMaj tmp2 = new DMatrixRMaj(2, 1);
            DMatrixRMaj tmp3 = new DMatrixRMaj(2, 1);
            CommonOps_DDRM.mult(PhiMinusI, B, tmp2);
            CommonOps_DDRM.mult(Ainv, tmp2, tmp3);
            CommonOps_DDRM.multAdd(tmp3, u0, tmp1); // tmp1 += tmp3*u0

            DMatrixRMaj y = new DMatrixRMaj(1, 1);
            CommonOps_DDRM.mult(C, tmp1, y);
            CommonOps_DDRM.multAdd(D, u0, y);

            System.out.printf("t=%.2f, y=%.8f%n", t, y.get(0,0));
        }
    }
}

Recommended Java libraries for modern-control computations: EJML (linear algebra), and (when doing simulation-heavy work) Apache Commons Math ODE integrators.

MATLAB’s Control System Toolbox supports direct state-space simulation via ss and lsim. The same script also includes an optional (commented) programmatic Simulink model construction using a State-Space block and a Step source.

File: Chapter7_Lesson4.m

% Chapter7_Lesson4.m
% Modern Control — Chapter 7, Lesson 4
% Output Response Using x(t) and y(t) = Cx(t) + Du(t)
%
% Demonstrates:
% 1) Closed-form y(t) for constant u(t)=u0 when A is invertible
% 2) Simulation using lsim (state-space object)
% 3) Optional Simulink model construction (commented section)

clear; clc;

% Example: 2-state SISO system
A = [0 1; -2 -3];
B = [0; 1];
C = [1 0];
D = 0.25;   % direct feedthrough to show output jump for step changes

x0 = [1; 0];
u0 = 1;

t = linspace(0, 6, 601).';

% Closed-form for constant input (requires invertible A):
% x(t) = expm(A t) x0 + A^{-1}(expm(A t) - I) B u0
% y(t) = C x(t) + D u0
I = eye(size(A));
if rcond(A) < 1e-12
    error('A is not invertible; constant-input closed form not applicable.');
end
Ainv = inv(A);

y_closed = zeros(length(t), 1);
x_closed = zeros(length(t), 2);

for k = 1:length(t)
    Phi = expm(A*t(k));
    xk = Phi*x0 + (Ainv*(Phi - I)*B)*u0;
    yk = C*xk + D*u0;
    x_closed(k,:) = xk.';
    y_closed(k) = yk;
end

disp('Closed-form y(t) at a few times:');
for ti = [0 0.5 1 2 6]
    [~,k] = min(abs(t-ti));
    fprintf('t=%.2f, y=%.8f\n', t(k), y_closed(k));
end

% Simulation using lsim
sys = ss(A,B,C,D);
u = ones(size(t)); % step
[y_lsim, t_out, x_lsim] = lsim(sys, u, t, x0);

disp('lsim y(t) at a few times:');
for ti = [0 0.5 1 2 6]
    [~,k] = min(abs(t_out-ti));
    fprintf('t=%.2f, y=%.8f\n', t_out(k), y_lsim(k));
end

% Plot
figure; plot(t, y_closed, t_out, y_lsim, '--'); grid on;
xlabel('t (s)'); ylabel('y(t)');
legend('Closed-form', 'lsim');

% ------------------------------------------------------------
% Optional: Programmatically build a Simulink model for y(t)
% ------------------------------------------------------------
%{
model = 'Chapter7_Lesson4_Simulink';
new_system(model); open_system(model);

add_block('simulink/Sources/Step', [model '/Step']);
add_block('simulink/Continuous/State-Space', [model '/StateSpace']);
add_block('simulink/Sinks/Scope', [model '/Scope']);

set_param([model '/StateSpace'], 'A', mat2str(A), 'B', mat2str(B), 'C', mat2str(C), 'D', mat2str(D));

add_line(model, 'Step/1', 'StateSpace/1');
add_line(model, 'StateSpace/1', 'Scope/1');

set_param(model, 'StopTime', '6');
save_system(model);

sim(model);
%}

Mathematica directly supports matrix exponentials via MatrixExp. The notebook below implements the constant-input closed form and plots \( y(t) \).

File: Chapter7_Lesson4.nb

(* Chapter7_Lesson4.nb
   Modern Control — Chapter 7, Lesson 4
   Output Response Using x(t) and y(t) = C x(t) + D u(t)

   This is a plain-text Mathematica Notebook expression. Open it in Mathematica.
*)

Notebook[{
  Cell["Modern Control — Chapter 7, Lesson 4: Output Response", "Title"],
  Cell["Define a continuous-time LTI state-space model and compute y(t) for a step input.", "Text"],

  Cell[BoxData@ToBoxes[
    Row[{
      "A = { {0, 1}, {-2, -3} };\n",
      "B = { {0}, {1} };\n",
      "C = { {1, 0} };\n",
      "D = { {0.25} };\n",
      "x0 = { {1}, {0} };\n",
      "u0 = { {1} };\n",
      "tmax = 6;\n"
    }]
  ], "Input"],

  Cell["Closed-form (constant input) if A is invertible:", "Text"],
  Cell[BoxData@ToBoxes[
    Row[{
      "Phi[t_] := MatrixExp[A t];\n",
      "x[t_] := Phi[t].x0 + Inverse[A].(Phi[t] - IdentityMatrix[2]).B.u0;\n",
      "y[t_] := C.x[t] + D.u0;\n"
    }]
  ], "Input"],

  Cell["Evaluate y(t) at selected times:", "Text"],
  Cell[BoxData@ToBoxes[
    Row[{
      "Table[{tt, y[tt][[1,1]]}, {tt, {0, 0.5, 1, 2, 6} }]\n"
    }]
  ], "Input"],

  Cell["Plot y(t) over time:", "Text"],
  Cell[BoxData@ToBoxes[
    Row[{
      "Plot[Evaluate[y[t][[1,1]]], {t, 0, tmax}, PlotRange -> All, GridLines -> Automatic, AxesLabel -> {\"t\", \"y(t)\"}]\n"
    }]
  ], "Input"]
}]

The problems below reinforce the output-response decomposition and its mathematical consequences. Unless stated otherwise, assume \( u(\cdot) \) is piecewise continuous on \( [0,T] \).

Problem 1 (Derive the output response). Starting from \( \dot{x}(t)=A x(t)+B u(t) \) with \( x(0)=x_0 \), and the output equation \( y(t)=C x(t)+D u(t) \), derive an explicit formula for \( y(t) \) in terms of \( e^{At} \), \( x_0 \), and \( u(\cdot) \).

Solution. From the known state solution:

\[ x(t)=e^{At}x_0+\int_{0}^{t} e^{A(t-\tau)}B u(\tau)\,d\tau. \]

Substitute into \( y(t)=C x(t)+D u(t) \):

\[ y(t)=C e^{At}x_0+\int_{0}^{t} C e^{A(t-\tau)}B u(\tau)\,d\tau + D u(t). \]

Problem 2 (Show that jumps in y(t) depend only on D). Assume \( x(t) \) is continuous. Suppose the input has a jump at \( t=t_0 \): \( \Delta u = u(t_0^+) - u(t_0^-) \). Prove that \( \Delta y = D\,\Delta u \).

Solution. Write

\[ \Delta y = \big(Cx(t_0^+)+D u(t_0^+)\big) - \big(Cx(t_0^-)+D u(t_0^-)\big) = C\big(x(t_0^+)-x(t_0^-)\big) + D\Delta u. \]

Since \( x(t) \) is continuous, \( x(t_0^+)-x(t_0^-)=0 \), hence \( \Delta y = D\Delta u \).

Problem 3 (Closed form for constant input with invertible A). Let \( u(t)=u_0 \) for \( t \ge 0 \), and assume \( A \) is invertible. Show that \( \int_{0}^{t} e^{A s}\,ds = A^{-1}(e^{At}-I) \), and derive closed forms for \( x(t) \) and \( y(t) \).

Solution. Consider \( F(t)=A^{-1}(e^{At}-I) \). Then

\[ \dot{F}(t)=A^{-1}(A e^{At})=e^{At}, \qquad F(0)=0. \]

Thus \( F(t)=\int_{0}^{t}e^{A s}ds \). Substitute into the state solution:

\[ x(t)=e^{At}x_0 + \left(A^{-1}(e^{At}-I)\right)B u_0, \qquad y(t)=C x(t) + D u_0. \]

Equivalently:

\[ y(t)=C e^{At}x_0 + C A^{-1}(e^{At}-I)B u_0 + D u_0. \]

Problem 4 (Impulse response identification). Assume \( x_0=0 \). Define \( h(t)=C e^{At}B \) for \( t \ge 0 \). Show that the dynamic forced output equals \( (h*u)(t) \).

Solution. Starting from the forced output term:

\[ \int_{0}^{t} C e^{A(t-\tau)}B\,u(\tau)\,d\tau. \]

By definition of convolution with a causal kernel \( h(\cdot) \):

\[ (h*u)(t) \equiv \int_{0}^{t} h(t-\tau)\,u(\tau)\,d\tau = \int_{0}^{t} C e^{A(t-\tau)}B\,u(\tau)\,d\tau. \]

Problem 5 (Output at t=0 and the role of D). Evaluate \( y(0) \) in terms of \( x_0 \) and \( u(0) \). Then, for an input step from 0 to \( u_0 \) at \( t=0 \), compute \( y(0^+) - y(0^-) \).

Solution. Directly from the output equation:

\[ y(0)=C x(0)+D u(0)=C x_0 + D u(0). \]

For a step from 0 to \( u_0 \), we have \( u(0^-)=0 \) and \( u(0^+)=u_0 \). Since \( x(t) \) is continuous, \( x(0^+)=x(0^-)=x_0 \). Therefore:

\[ y(0^+)-y(0^-) = \big(Cx_0 + D u_0\big) - \big(Cx_0 + D\cdot 0\big) = D u_0. \]

We derived the output response directly from the known state solution, obtaining \( y(t)=C e^{At}x_0 + \int_{0}^{t} C e^{A(t-\tau)}B u(\tau)d\tau + D u(t) \). This decomposition cleanly separates initial-condition effects, dynamic forced effects (a convolution with \( h(t)=C e^{At}B \)), and instantaneous feedthrough. We proved continuity of the state under mild input assumptions and derived the jump condition \( \Delta y = D \Delta u \).