Modern control is matrix-centered. Even a small physical model produces coupled matrix equations whose theoretical properties must be checked numerically. For a continuous-time LTI model,

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

software is used to compute eigenstructure, simulate trajectories, assess controllability and observability, design feedback gains, and compare physical models against linearized state-space realizations. A correct workflow should distinguish symbolic theory, floating-point linear algebra, and engineering validation.

The main software families used in this lesson are MATLAB/Simulink, Python, C++, Java, and Wolfram Mathematica. MATLAB provides mature control-design functions such as ss, ctrb, obsv, place, lqr, gram, lsim, and Simulink linearization workflows. Python typically combines NumPy, SciPy, python-control, Slycot, CVXPY, and CasADi. C++ often uses Eigen, Armadillo, or SLICOT-linked routines; Java often uses EJML, Apache Commons Math, Hipparchus, or ojAlgo. Mathematica is useful for exact symbolic checks, closed-form matrix exponentials, and notebook derivations.

The transfer matrix associated with a state-space realization is

\[ G(s)=C(sI-A)^{-1}B+D. \]

Internal analysis requires more than the transfer function. For example, controllability and observability are checked by the Kalman matrices

\[ \mathcal{C}(A,B)=\begin{bmatrix}B&AB&\cdots&A^{n-1}B\end{bmatrix}, \qquad \mathcal{O}(A,C)=\begin{bmatrix}C\\CA\\\vdots\\CA^{n-1}\end{bmatrix}. \]

The exact algebraic criteria are

\[ \operatorname{rank}\mathcal{C}(A,B)=n,\qquad \operatorname{rank}\mathcal{O}(A,C)=n. \]

In floating-point software the computed rank depends on a tolerance. A common singular-value rule is

\[ \widehat{r}=\#\{i:\sigma_i(M)>\epsilon\,\max(m,n)\,\sigma_1(M)\}. \]

Therefore, a system may be theoretically controllable but numerically ill-conditioned. The condition number \( \kappa_2(M)=\sigma_1(M)/\sigma_r(M) \) is often more informative than rank alone.

For stable continuous-time systems, controllability and observability Gramians solve Lyapunov equations:

\[ AW_c+W_cA^T+BB^T=0,\qquad A^TW_o+W_oA+C^TC=0. \]

The finite-time reachability Gramian is

\[ W_c(t_f)=\int_0^{t_f}e^{A\theta}BB^Te^{A^T\theta}\,d\theta. \]

If \( W_c(t_f) \) is nonsingular, the minimum-energy input that transfers \( \mathbf{x}(0)=\mathbf{0} \) to \( \mathbf{x}(t_f)=\mathbf{x}_f \) is

\[ \mathbf{u}^*(t)=B^Te^{A^T(t_f-t)}W_c(t_f)^{-1}\mathbf{x}_f, \qquad J^*=\mathbf{x}_f^TW_c(t_f)^{-1}\mathbf{x}_f. \]

A reliable toolchain starts from a physical or symbolic model, then passes through numerical checks before controller synthesis. The workflow below is independent of whether the implementation is in MATLAB, Python, C++, Java, or Mathematica.

flowchart TD
  A["Physical or symbolic model"] --> B["Build A, B, C, D matrices"]
  B --> C["Check units, signs, scaling"]
  C --> D["Analyze eigenvalues, rank tests, gramians"]
  D --> E["Choose design method: place or LQR"]
  E --> F["Compute K and closed-loop matrix"]
  F --> G["Verify poles, residuals, conditioning"]
  G --> H["Simulate time response and input effort"]
  H --> I["Export report, code, and reproducible data"]
        

Each stage should produce a checkable numerical artifact. For example, a pole-placement routine should not only return \(K\); it should also verify the eigenvalues of \(A-BK\) and measure residuals.

\[ r_i=\|(A-BK)\mathbf{v}_i-\lambda_i\mathbf{v}_i\|_2,\qquad \max_i r_i \;\text{small relative to}\; \|A-BK\|_2. \]

Software frequently changes coordinates for balancing, scaling, realization conversion, or canonical form construction. Let \( \mathbf{x}=T\mathbf{z} \) with nonsingular \(T\). Then

\[ \dot{\mathbf{z}}=T^{-1}AT\mathbf{z}+T^{-1}B\mathbf{u}, \qquad \mathbf{y}=CT\mathbf{z}+D\mathbf{u}. \]

Hence the transformed realization is

\[ \tilde{A}=T^{-1}AT,\qquad \tilde{B}=T^{-1}B, \qquad \tilde{C}=CT, \qquad \tilde{D}=D. \]

Proof of transfer-matrix invariance.

\[ \begin{aligned} \tilde{G}(s) &=\tilde{C}(sI-\tilde{A})^{-1}\tilde{B}+\tilde{D}\\ &=CT\left(sI-T^{-1}AT\right)^{-1}T^{-1}B+D\\ &=CT\left(T^{-1}(sI-A)T\right)^{-1}T^{-1}B+D\\ &=CT\left(T^{-1}(sI-A)^{-1}T\right)T^{-1}B+D\\ &=C(sI-A)^{-1}B+D=G(s). \end{aligned} \]

Thus MATLAB, Python, or Mathematica may internally use a transformed realization without changing input-output behavior. However, \( \mathcal{C}(A,B) \) and \( \mathcal{O}(A,C) \) can have very different condition numbers after a poor scaling, even though their exact ranks are invariant.

Pole placement computes a gain \(K\) such that the closed-loop matrix \(A-BK\) has a prescribed spectrum:

\[ \operatorname{spec}(A-BK)=\{\lambda_1,\ldots,\lambda_n\}, \qquad \operatorname{Re}(\lambda_i)<0. \]

For a controllable SISO system, Ackermann's formula is

\[ K=\mathbf{e}_n^T\mathcal{C}(A,B)^{-1}\phi_d(A),\qquad \phi_d(s)=s^n+\alpha_{n-1}s^{n-1}+\cdots+\alpha_0. \]

where \( \mathbf{e}_n^T=\begin{bmatrix}0&\cdots&0&1\end{bmatrix} \) and

\[ \phi_d(A)=A^n+\alpha_{n-1}A^{n-1}+\cdots+\alpha_1A+ \alpha_0I. \]

For quadratic optimal state feedback, continuous-time LQR solves

\[ \min_{\mathbf{u}}\int_0^\infty \left(\mathbf{x}^TQ\mathbf{x}+\mathbf{u}^TR\mathbf{u}\right)dt, \qquad Q\succeq0,\quad R\succ0. \]

The stabilizing solution \(P=P^T\succeq0\) of the continuous algebraic Riccati equation satisfies

\[ A^TP+PA-PBR^{-1}B^TP+Q=0,\qquad K=R^{-1}B^TP. \]

For simulation, both MATLAB and Python rely on the state-transition formula

\[ \mathbf{x}(t)=e^{At}\mathbf{x}(0)+\int_0^t e^{A(t-\sigma)}B\mathbf{u}(\sigma)\,d\sigma. \]

For digital implementation with zero-order-hold sampling period \(h\), the exact discrete model is

\[ A_d=e^{Ah},\qquad B_d=\int_0^h e^{A\sigma}B\,d\sigma. \]

Software tools can return a gain even when the design is numerically fragile. Before accepting a design, check scaling, ranks, conditioning, closed-loop residuals, simulation consistency, and input magnitude.

flowchart TD
  S["Scale states and inputs"] --> R["Rank tests with singular values"]
  R --> C["Condition numbers of Ctrb and Obsv"]
  C --> P["Design gain K"]
  P --> V["Verify poles and eigenvector residuals"]
  V --> T["Time-domain simulation"]
  T --> U["Check actuator effort and saturation"]
  U --> REP["Save scripts, random seeds, versions"]
        

A diagonal scaling \( \mathbf{x}=S\mathbf{z} \) produces

\[ A_s=S^{-1}AS, \qquad B_s=S^{-1}B, \qquad C_s=CS. \]

If the controller is designed in scaled coordinates as \( \mathbf{u}=-K_s\mathbf{z} \), then the physical feedback is

\[ \mathbf{u}=-K_sS^{-1}\mathbf{x}, \qquad K=K_sS^{-1}. \]

A useful acceptance test is to compare the achieved characteristic polynomial against the desired polynomial:

\[ \rho=\frac{\|\operatorname{coeff}(\det(sI-A+BK))-\operatorname{coeff}(\phi_d(s))\|_2} {\|\operatorname{coeff}(\phi_d(s))\|_2}. \]

In a robust educational workflow, students should report \( \operatorname{rank}\mathcal{C} \), \( \kappa_2(\mathcal{C}) \), closed-loop poles, simulation plots, and software version numbers.

The Python file below uses NumPy for matrix construction, SciPy for pole placement, Riccati equations, matrix exponentials, and simulation, and Matplotlib for response plots. It also implements controllability, observability, and a finite-horizon reachability Gramian directly.

Chapter30_Lesson3.py

"""
Chapter30_Lesson3.py
Software tools for state-space analysis and design in Python.

Dependencies:
    pip install numpy scipy matplotlib control
Optional:
    pip install slycot
"""

import numpy as np
from numpy.linalg import matrix_rank, eigvals, cond
from scipy.linalg import expm, solve_continuous_are
from scipy.signal import StateSpace, lsim, place_poles
import matplotlib.pyplot as plt


def controllability_matrix(A: np.ndarray, B: np.ndarray) -> np.ndarray:
    """Return C = [B, AB, ..., A^(n-1)B]."""
    n = A.shape[0]
    blocks = [B]
    for k in range(1, n):
        blocks.append(np.linalg.matrix_power(A, k) @ B)
    return np.hstack(blocks)


def observability_matrix(A: np.ndarray, C: np.ndarray) -> np.ndarray:
    """Return O = [C; CA; ...; CA^(n-1)]."""
    n = A.shape[0]
    blocks = [C]
    for k in range(1, n):
        blocks.append(C @ np.linalg.matrix_power(A, k))
    return np.vstack(blocks)


def lqr_gain(A: np.ndarray, B: np.ndarray, Q: np.ndarray, R: np.ndarray) -> np.ndarray:
    """Continuous-time LQR gain K = R^(-1) B' P."""
    P = solve_continuous_are(A, B, Q, R)
    return np.linalg.solve(R, B.T @ P)


def finite_horizon_reachability_gramian(A: np.ndarray, B: np.ndarray, tf: float, steps: int = 4000) -> np.ndarray:
    """Numerically compute Wc(tf) = integral_0^tf exp(A t) B B' exp(A' t) dt."""
    n = A.shape[0]
    W = np.zeros((n, n))
    h = tf / steps
    for i in range(steps + 1):
        t = i * h
        Phi = expm(A * t)
        integrand = Phi @ B @ B.T @ Phi.T
        weight = 0.5 if i in (0, steps) else 1.0
        W += weight * integrand
    return h * W


def main() -> None:
    # Mass-spring-damper example: x1 = position, x2 = velocity
    m, c, k = 1.0, 0.35, 2.0
    A = np.array([[0.0, 1.0], [-k / m, -c / m]])
    B = np.array([[0.0], [1.0 / m]])
    C = np.array([[1.0, 0.0]])
    D = np.array([[0.0]])

    Ctrb = controllability_matrix(A, B)
    Obsv = observability_matrix(A, C)
    print("A eigenvalues:", eigvals(A))
    print("rank(Ctrb):", matrix_rank(Ctrb), "condition(Ctrb):", cond(Ctrb))
    print("rank(Obsv):", matrix_rank(Obsv), "condition(Obsv):", cond(Obsv))

    desired_poles = np.array([-2.0 + 1.5j, -2.0 - 1.5j])
    K_place = place_poles(A, B, desired_poles).gain_matrix
    print("Pole-placement K:", K_place)
    print("closed-loop poles:", eigvals(A - B @ K_place))

    Q = np.diag([20.0, 2.0])
    R = np.array([[0.5]])
    K_lqr = lqr_gain(A, B, Q, R)
    print("LQR K:", K_lqr)
    print("LQR poles:", eigvals(A - B @ K_lqr))

    Wc = finite_horizon_reachability_gramian(A, B, tf=5.0)
    print("Wc(5 s):\n", Wc)
    print("energy to reach xf=[1,0] approximately:", np.array([[1.0, 0.0]]) @ np.linalg.inv(Wc) @ np.array([[1.0], [0.0]]))

    # Closed-loop free response from an initial displacement.
    Acl = A - B @ K_lqr
    sys_cl = StateSpace(Acl, B, C, D)
    t = np.linspace(0.0, 8.0, 600)
    u = np.zeros_like(t)
    tout, y, x = lsim(sys_cl, U=u, T=t, X0=np.array([1.0, 0.0]))

    plt.figure()
    plt.plot(tout, x[:, 0], label="position")
    plt.plot(tout, x[:, 1], label="velocity")
    plt.xlabel("time [s]")
    plt.ylabel("state")
    plt.title("Chapter 30 Lesson 3: LQR closed-loop response")
    plt.grid(True)
    plt.legend()
    plt.show()


if __name__ == "__main__":
    main()

      

MATLAB is usually the fastest environment for teaching state-space workflows because its control functions directly match the notation in the theory. Simulink should be used to validate block-level models, nonlinear models, saturation effects, and implementation timing.

Chapter30_Lesson3.m

% Chapter30_Lesson3.m
% MATLAB/Simulink workflow for state-space analysis and design.
% Toolboxes commonly used: Control System Toolbox, Robust Control Toolbox,
% Simulink Control Design, Symbolic Math Toolbox, Model Predictive Control Toolbox.

clear; clc; close all;

% Mass-spring-damper plant: x1 = position, x2 = velocity
m = 1.0; c = 0.35; k = 2.0;
A = [0 1; -k/m -c/m];
B = [0; 1/m];
C = [1 0];
D = 0;
sys = ss(A, B, C, D);

fprintf('Open-loop poles:\n'); disp(eig(A));

Co = ctrb(A, B);
Ob = obsv(A, C);
fprintf('rank(ctrb) = %d, cond(ctrb) = %.3e\n', rank(Co), cond(Co));
fprintf('rank(obsv) = %d, cond(obsv) = %.3e\n', rank(Ob), cond(Ob));

% Pole placement design
p_des = [-2+1.5i, -2-1.5i];
K_place = place(A, B, p_des);
Acl_place = A - B*K_place;
fprintf('K_place =\n'); disp(K_place);
fprintf('Closed-loop poles by place:\n'); disp(eig(Acl_place));

% Continuous-time LQR design
Q = diag([20, 2]);
R = 0.5;
K_lqr = lqr(A, B, Q, R);
Acl_lqr = A - B*K_lqr;
fprintf('K_lqr =\n'); disp(K_lqr);
fprintf('Closed-loop poles by LQR:\n'); disp(eig(Acl_lqr));

% Gramian analysis
Wc = gram(sys, 'c');
Wo = gram(sys, 'o');
fprintf('Controllability Gramian:\n'); disp(Wc);
fprintf('Observability Gramian:\n'); disp(Wo);

% Closed-loop simulation from initial condition
sys_cl = ss(Acl_lqr, B, C, D);
t = linspace(0, 8, 600);
x0 = [1; 0];
[y, t, x] = initial(sys_cl, x0, t);
figure;
plot(t, x(:,1), 'LineWidth', 1.5); hold on;
plot(t, x(:,2), 'LineWidth', 1.5);
grid on; xlabel('time [s]'); ylabel('state');
legend('position', 'velocity');
title('Chapter 30 Lesson 3: LQR closed-loop response');

% Simulink note:
% 1. Insert a State-Space block with matrices A, B, C, D.
% 2. Add a Gain block with value -K_lqr in feedback from state output.
% 3. If the physical model is built from blocks, use linearize or linmod to obtain ss(A,B,C,D).
% 4. Compare Simulink output against initial(sys_cl,x0,t) for verification.

      

C++ is appropriate when the controller must be embedded, compiled, or integrated into a real-time simulator. The example uses Eigen and implements controllability, observability, SISO Ackermann pole placement, and RK4 simulation.

Chapter30_Lesson3.cpp

/*
Chapter30_Lesson3.cpp
State-space analysis with Eigen: controllability, observability, Ackermann pole placement, RK4 simulation.

Build example:
    g++ -std=c++17 Chapter30_Lesson3.cpp -I /path/to/eigen -O2 -o Chapter30_Lesson3
*/

#include <Eigen/Dense>
#include <complex>
#include <iostream>
#include <vector>

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

MatrixXd controllabilityMatrix(const MatrixXd& A, const MatrixXd& B) {
    const int n = A.rows();
    const int m = B.cols();
    MatrixXd C(n, n * m);
    MatrixXd Ak = MatrixXd::Identity(n, n);
    for (int k = 0; k < n; ++k) {
        C.block(0, k * m, n, m) = Ak * B;
        Ak = A * Ak;
    }
    return C;
}

MatrixXd observabilityMatrix(const MatrixXd& A, const MatrixXd& Cmat) {
    const int n = A.rows();
    const int p = Cmat.rows();
    MatrixXd O(n * p, n);
    MatrixXd Ak = MatrixXd::Identity(n, n);
    for (int k = 0; k < n; ++k) {
        O.block(k * p, 0, p, n) = Cmat * Ak;
        Ak = Ak * A;
    }
    return O;
}

VectorXd characteristicFromDesiredPoles(const std::vector<std::complex<double>>& poles) {
    // Returns coefficients [a0, a1, ..., a_{n-1}] for s^n + a_{n-1}s^{n-1}+...+a0.
    std::vector<std::complex<double>> coeff{1.0};
    for (const auto& pole : poles) {
        std::vector<std::complex<double>> next(coeff.size() + 1, 0.0);
        for (std::size_t i = 0; i < coeff.size(); ++i) {
            next[i] += -pole * coeff[i];
            next[i + 1] += coeff[i];
        }
        coeff = next;
    }
    VectorXd realCoeff(static_cast<int>(poles.size()));
    for (int i = 0; i < static_cast<int>(poles.size()); ++i) {
        realCoeff(i) = coeff[i].real();
    }
    return realCoeff;
}

MatrixXd matrixPolynomial(const MatrixXd& A, const VectorXd& coeff) {
    // phi(A) = A^n + a_{n-1}A^{n-1}+...+a0 I, with coeff=[a0,...,a_{n-1}].
    const int n = A.rows();
    MatrixXd phi = MatrixXd::Zero(n, n);
    MatrixXd Ak = MatrixXd::Identity(n, n);
    for (int i = 0; i < n; ++i) {
        phi += coeff(i) * Ak;
        Ak = A * Ak;
    }
    phi += Ak;
    return phi;
}

MatrixXd ackermannGainSISO(const MatrixXd& A, const MatrixXd& B, const std::vector<std::complex<double>>& poles) {
    const int n = A.rows();
    MatrixXd Ctrb = controllabilityMatrix(A, B);
    Eigen::FullPivLU<MatrixXd> lu(Ctrb);
    if (lu.rank() != n) {
        throw std::runtime_error("System is not controllable; Ackermann formula is not valid.");
    }
    VectorXd coeff = characteristicFromDesiredPoles(poles);
    MatrixXd phiA = matrixPolynomial(A, coeff);
    MatrixXd eT = MatrixXd::Zero(1, n);
    eT(0, n - 1) = 1.0;
    return eT * Ctrb.inverse() * phiA;
}

VectorXd fClosedLoop(const MatrixXd& A, const MatrixXd& B, const MatrixXd& K, const VectorXd& x) {
    return (A - B * K) * x;
}

VectorXd rk4Step(const MatrixXd& A, const MatrixXd& B, const MatrixXd& K, const VectorXd& x, double h) {
    VectorXd k1 = fClosedLoop(A, B, K, x);
    VectorXd k2 = fClosedLoop(A, B, K, x + 0.5 * h * k1);
    VectorXd k3 = fClosedLoop(A, B, K, x + 0.5 * h * k2);
    VectorXd k4 = fClosedLoop(A, B, K, x + h * k3);
    return x + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
}

int main() {
    MatrixXd A(2, 2), B(2, 1), C(1, 2);
    A << 0.0, 1.0,
        -2.0, -0.35;
    B << 0.0,
        1.0;
    C << 1.0, 0.0;

    MatrixXd Ctrb = controllabilityMatrix(A, B);
    MatrixXd Obsv = observabilityMatrix(A, C);
    std::cout << "Controllability matrix:\n" << Ctrb << "\nrank = " << Eigen::FullPivLU<MatrixXd>(Ctrb).rank() << "\n\n";
    std::cout << "Observability matrix:\n" << Obsv << "\nrank = " << Eigen::FullPivLU<MatrixXd>(Obsv).rank() << "\n\n";

    std::vector<std::complex<double>> poles{{-2.0, 1.5}, {-2.0, -1.5}};
    MatrixXd K = ackermannGainSISO(A, B, poles);
    std::cout << "Ackermann K = " << K << "\n";
    std::cout << "A - B*K:\n" << A - B * K << "\n\n";

    VectorXd x(2);
    x << 1.0, 0.0;
    double h = 0.01;
    for (int i = 0; i <= 500; ++i) {
        if (i % 100 == 0) {
            std::cout << "t=" << i * h << ", x=[" << x.transpose() << "]\n";
        }
        x = rk4Step(A, B, K, x, h);
    }
    return 0;
}

      

Java is less common in academic control courses than MATLAB or Python, but it is useful for server-side simulation tools, educational web backends, and industrial middleware. The example below is deliberately written without external dependencies so students can inspect every operation. For larger systems, use EJML or Hipparchus.

Chapter30_Lesson3.java

/*
Chapter30_Lesson3.java
Pure Java state-space utilities for small educational models.

Production libraries to consider:
    EJML, Apache Commons Math, Hipparchus, ojAlgo
*/

import java.util.Arrays;

public class Chapter30_Lesson3 {
    static double[][] matMul(double[][] A, double[][] B) {
        int n = A.length, m = B[0].length, r = B.length;
        double[][] C = new double[n][m];
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                double s = 0.0;
                for (int k = 0; k < r; k++) s += A[i][k] * B[k][j];
                C[i][j] = s;
            }
        }
        return C;
    }

    static double[][] matAdd(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[][] A) {
        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 * A[i][j];
        return C;
    }

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

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

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

    static double[][] controllability(double[][] A, double[][] B) {
        int n = A.length;
        double[][] Ctrb = B;
        double[][] Ak = eye(n);
        for (int k = 1; k < n; k++) {
            Ak = matMul(A, Ak);
            Ctrb = hcat(Ctrb, matMul(Ak, B));
        }
        return Ctrb;
    }

    static double[][] observability(double[][] A, double[][] C) {
        int n = A.length;
        double[][] Obsv = C;
        double[][] Ak = eye(n);
        for (int k = 1; k < n; k++) {
            Ak = matMul(Ak, A);
            Obsv = vcat(Obsv, matMul(C, Ak));
        }
        return Obsv;
    }

    static double det2(double[][] M) {
        return M[0][0] * M[1][1] - M[0][1] * M[1][0];
    }

    static double[][] inv2(double[][] M) {
        double det = det2(M);
        if (Math.abs(det) < 1e-12) throw new RuntimeException("Singular 2x2 matrix");
        return new double[][]{{M[1][1] / det, -M[0][1] / det}, {-M[1][0] / det, M[0][0] / det}};
    }

    static double[][] ackermann2x1(double[][] A, double[][] B, double pole1, double pole2) {
        // Desired polynomial: s^2 + a1 s + a0.
        double a1 = -(pole1 + pole2);
        double a0 = pole1 * pole2;
        double[][] A2 = matMul(A, A);
        double[][] phiA = matAdd(matAdd(A2, scale(a1, A)), scale(a0, eye(2)));
        double[][] CtrbInv = inv2(controllability(A, B));
        double[][] eT = new double[][]{{0.0, 1.0}};
        return matMul(matMul(eT, CtrbInv), phiA);
    }

    static double[] closedLoopDerivative(double[][] A, double[][] B, double[][] K, double[] x) {
        double[][] Acl = matAdd(A, scale(-1.0, matMul(B, K)));
        return new double[]{Acl[0][0] * x[0] + Acl[0][1] * x[1], Acl[1][0] * x[0] + Acl[1][1] * x[1]};
    }

    static double[] rk4(double[][] A, double[][] B, double[][] K, double[] x, double h) {
        double[] k1 = closedLoopDerivative(A, B, K, x);
        double[] k2 = closedLoopDerivative(A, B, K, new double[]{x[0] + 0.5 * h * k1[0], x[1] + 0.5 * h * k1[1]});
        double[] k3 = closedLoopDerivative(A, B, K, new double[]{x[0] + 0.5 * h * k2[0], x[1] + 0.5 * h * k2[1]});
        double[] k4 = closedLoopDerivative(A, B, K, new double[]{x[0] + h * k3[0], x[1] + h * k3[1]});
        return new double[]{x[0] + h * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6.0,
                x[1] + h * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]) / 6.0};
    }

    static void print(String title, double[][] A) {
        System.out.println(title);
        for (double[] row : A) System.out.println(Arrays.toString(row));
    }

    public static void main(String[] args) {
        double[][] A = {{0.0, 1.0}, {-2.0, -0.35}};
        double[][] B = {{0.0}, {1.0}};
        double[][] C = {{1.0, 0.0}};

        print("Controllability matrix", controllability(A, B));
        print("Observability matrix", observability(A, C));

        double[][] K = ackermann2x1(A, B, -2.0, -3.0);
        print("K from Ackermann", K);

        double[] x = {1.0, 0.0};
        double h = 0.01;
        for (int i = 0; i <= 500; i++) {
            if (i % 100 == 0) System.out.printf("t=%.2f, x=[%.6f, %.6f]%n", i * h, x[0], x[1]);
            x = rk4(A, B, K, x, h);
        }
    }
}

      

Mathematica is useful when one wants exact characteristic polynomials, symbolic matrix exponentials, and notebook-based derivations. The downloadable file is a notebook expression that can be opened and evaluated in Mathematica.

Chapter30_Lesson3.nb

Notebook[{
Cell["Chapter30 Lesson3: State-Space Analysis and Design", "Title"],
Cell["This notebook computes controllability, observability, pole placement data, LQR equations, and a closed-loop simulation for a mass-spring-damper system.", "Text"],
Cell[BoxData[RowBox[{RowBox[{"A", "=", RowBox[{"{{", RowBox[{RowBox[{"0", ",", "1"}], ",", RowBox[{RowBox[{"-", "2"}], ",", RowBox[{"-", "0.35"}]}]}], "}}"}]}], ";"}]], "Input"],
Cell[BoxData[RowBox[{RowBox[{"B", "=", RowBox[{"{{", RowBox[{"0"}], "},{", RowBox[{"1"}], "}}"}]}], ";"}]], "Input"],
Cell[BoxData[RowBox[{RowBox[{"Cmat", "=", RowBox[{"{{", RowBox[{"1", ",", "0"}], "}}"}]}], ";"}]], "Input"],
Cell[BoxData[RowBox[{"Eigenvalues", "[", "A", "]"}]], "Input"],
Cell[BoxData[RowBox[{RowBox[{"Ctrb", "=", RowBox[{"ArrayFlatten", "[", RowBox[{"{", RowBox[{"{", RowBox[{"B", ",", RowBox[{"A", ".", "B"}]}], "}"}], "}"}], "]"}]}], ";", RowBox[{"MatrixRank", "[", "Ctrb", "]"}]}]], "Input"],
Cell[BoxData[RowBox[{RowBox[{"Obsv", "=", RowBox[{"Join", "[", RowBox[{"Cmat", ",", RowBox[{"Cmat", ".", "A"}]}], "]"}]}], ";", RowBox[{"MatrixRank", "[", "Obsv", "]"}]}]], "Input"],
Cell[BoxData[RowBox[{RowBox[{"desiredPoly", "=", RowBox[{"Expand", "[", RowBox[{RowBox[{"(", RowBox[{"s", "+", "2", "-", RowBox[{"1.5", "I"}]}], ")"}], RowBox[{"(", RowBox[{"s", "+", "2", "+", RowBox[{"1.5", "I"}]}], ")"}]}], "]"}]}], ";", "desiredPoly"}]], "Input"],
Cell[BoxData[RowBox[{RowBox[{"K", "=", RowBox[{"{{", RowBox[{"0", ",", "1"}], "}}", ".", RowBox[{"Inverse", "[", "Ctrb", "]"}], ".", RowBox[{"(", RowBox[{RowBox[{"MatrixPower", "[", RowBox[{"A", ",", "2"}], "]"}], "+", RowBox[{"4", "A"}], "+", RowBox[{"6.25", RowBox[{"IdentityMatrix", "[", "2", "]"}]}]}], ")"}]}]}], ";", "K"}]], "Input"],
Cell[BoxData[RowBox[{"Eigenvalues", "[", RowBox[{"A", "-", RowBox[{"B", ".", "K"}]}], "]"}]], "Input"],
Cell[BoxData[RowBox[{"NDSolveValue", "[", RowBox[{RowBox[{"{", RowBox[{RowBox[{RowBox[{"x1", "'"}], "[", "t", "]"}], "==", RowBox[{RowBox[{"x2", "[", "t", "]"}]}], ",", RowBox[{RowBox[{RowBox[{"x2", "'"}], "[", "t", "]"}], "==", RowBox[{RowBox[{"-", RowBox[{"(", RowBox[{"2", "+", RowBox[{"K", "[[", RowBox[{"1", ",", "1"}], "]]"}]}], ")"}]}], RowBox[{"x1", "[", "t", "]"}], "-", RowBox[{RowBox[{"(", RowBox[{"0.35", "+", RowBox[{"K", "[[", RowBox[{"1", ",", "2"}], "]]"}]}], ")"}], RowBox[{"x2", "[", "t", "]"}]}]}], ",", RowBox[{RowBox[{"x1", "[", "0", "]"}], "==", "1"}], ",", RowBox[{RowBox[{"x2", "[", "0", "]"}], "==", "0"}]}], "}"}], ",", RowBox[{"{", RowBox[{"x1", ",", "x2"}], "}"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "8"}], "}"}]}], "]"}]], "Input"]
}]

      

Problem 1 (Similarity Invariance): Let \( \tilde{A}=T^{-1}AT \), \( \tilde{B}=T^{-1}B \), and \( \tilde{C}=CT \). Prove that controllability is invariant under this coordinate change.

Solution: The transformed controllability matrix is

\[ \begin{aligned} \mathcal{C}(\tilde{A},\tilde{B}) &=\begin{bmatrix}T^{-1}B&T^{-1}AB&\cdots&T^{-1}A^{n-1}B\end{bmatrix}\\ &=T^{-1}\mathcal{C}(A,B). \end{aligned} \]

Since \(T^{-1}\) is nonsingular, left multiplication by \(T^{-1}\) does not change rank. Therefore \( \operatorname{rank}\mathcal{C}(\tilde{A},\tilde{B})= \operatorname{rank}\mathcal{C}(A,B) \).

Problem 2 (Ackermann Calculation): For \( A=\begin{bmatrix}0&1\\-2&-0.35\end{bmatrix} \) and \( B=\begin{bmatrix}0\\1\end{bmatrix} \), find the SISO pole-placement gain for desired poles \( -2+j1.5 \) and \( -2-j1.5 \).

Solution: The desired polynomial is

\[ \phi_d(s)=(s+2-j1.5)(s+2+j1.5)=s^2+4s+6.25. \]

The controllability matrix is

\[ \mathcal{C}=\begin{bmatrix}0&1\\1&-0.35\end{bmatrix}. \]

Using \(K=\mathbf{e}_2^T\mathcal{C}^{-1}\phi_d(A)\) gives

\[ K=\begin{bmatrix}4.25&3.65\end{bmatrix}. \]

Then \(A-BK=\begin{bmatrix}0&1\\-6.25&-4\end{bmatrix}\), whose characteristic polynomial is exactly \(s^2+4s+6.25\).

Problem 3 (Continuous-Time LQR): State the Riccati equation solved by MATLAB lqr and Python solve_continuous_are, and explain why \(R\succ0\) is required.

Solution: The equation is

\[ A^TP+PA-PBR^{-1}B^TP+Q=0. \]

The gain is \(K=R^{-1}B^TP\). The matrix \(R\) must be positive definite so that the control penalty is strictly convex and \(R^{-1}\) exists. If \(R\) is singular, the standard finite-dimensional unconstrained LQR formula is not directly valid.

Problem 4 (Rank Tolerance): Suppose \(\sigma(\mathcal{C})=\{10,10^{-10}\}\). Explain why different software packages may disagree on controllability.

Solution: Exact rank is two if both singular values are nonzero. Numerically, however, a rank algorithm compares singular values against a tolerance. If the tolerance is larger than \(10^{-10}\), the computed rank is one; if the tolerance is smaller, the computed rank is two. This is not a contradiction: it indicates numerical near-uncontrollability and large steering energy in at least one state direction.

Problem 5 (Discrete-Time Conversion): Show the exact zero-order-hold formulas used to convert \(\dot{\mathbf{x}}=A\mathbf{x}+B\mathbf{u}\) to \(\mathbf{x}_{k+1}=A_d\mathbf{x}_k+B_d\mathbf{u}_k\).

Solution: With constant input over one sample period \(h\), integration gives

\[ \mathbf{x}((k+1)h)=e^{Ah}\mathbf{x}(kh)+ \int_0^h e^{A\sigma}B\,d\sigma\;\mathbf{u}_k. \]

Therefore

\[ A_d=e^{Ah},\qquad B_d=\int_0^h e^{A\sigma}B\,d\sigma. \]

State-space software should be used as a mathematical instrument, not as a black box. MATLAB and Python provide the strongest teaching workflows for state-space analysis, while C++, Java, and Mathematica support embedded implementation, backend simulation, and symbolic verification. A professional workflow checks rank, condition number, closed-loop poles, residuals, time-domain response, input effort, coordinate scaling, and reproducibility.