In previous lessons, you linearized scalar or single-input models using a Jacobian about an operating point. In control engineering, realistic plants are often coupled: several actuators influence several measured outputs. MIMO linearization produces an LTI approximation suitable for analysis and controller design near an operating point.

The core idea is unchanged: introduce perturbations \( \tilde{\mathbf{x}}(t) \): deviation from an operating state, \( \tilde{\mathbf{u}}(t) \): deviation from an operating input, and \( \tilde{\mathbf{y}}(t) \): deviation from an operating output, then retain only first-order terms.

flowchart TD
  A["Nonlinear model: xdot=f(x,u), y=h(x,u)"] --> B["Choose operating point (x0,u0)"]
  B --> C["Define perturbations: x_tilde=x-x0, u_tilde=u-u0, y_tilde=y-y0"]
  C --> D["Compute Jacobians at (x0,u0): A,B,C,D"]
  D --> E["Linear model: x_tilde_dot = A x_tilde + B u_tilde"]
  E --> F["Output: y_tilde = C x_tilde + D u_tilde"]
  F --> G["Local validation by simulation"]
        

The distinctive MIMO aspect is that \( \mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D} \) are matrices, and their off-diagonal structure reveals cross-coupling between channels.

Consider a smooth nonlinear state-space model with \( \mathbf{x}(t)\in\mathbb{R}^n \) (states), \( \mathbf{u}(t)\in\mathbb{R}^m \) (inputs), \( \mathbf{y}(t)\in\mathbb{R}^p \) (outputs):

\[ \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x},\mathbf{u}), \qquad \mathbf{y} = \mathbf{h}(\mathbf{x},\mathbf{u}), \qquad \mathbf{f}:\mathbb{R}^n\times\mathbb{R}^m → \mathbb{R}^n,\; \mathbf{h}:\mathbb{R}^n\times\mathbb{R}^m → \mathbb{R}^p. \]

An operating point is a constant pair \( (\mathbf{x}_0,\mathbf{u}_0) \):. If it additionally satisfies \( \mathbf{f}(\mathbf{x}_0,\mathbf{u}_0)=\mathbf{0} \):, it is an equilibrium (steady state).

Define perturbation variables (also called small-signal variables): \( \tilde{\mathbf{x}} \): \( \mathbf{x}-\mathbf{x}_0 \), \( \tilde{\mathbf{u}} \): \( \mathbf{u}-\mathbf{u}_0 \), \( \tilde{\mathbf{y}} \): \( \mathbf{y}-\mathbf{y}_0 \), where \( \mathbf{y}_0=\mathbf{h}(\mathbf{x}_0,\mathbf{u}_0) \):.

\[ \mathbf{x}=\mathbf{x}_0+\tilde{\mathbf{x}},\qquad \mathbf{u}=\mathbf{u}_0+\tilde{\mathbf{u}},\qquad \mathbf{y}=\mathbf{y}_0+\tilde{\mathbf{y}}. \]

Linearization is the first-order truncation of a multivariable Taylor expansion. It is helpful to combine state and input into a single vector \( \mathbf{z} \): where \( \mathbf{z}=\begin{bmatrix}\mathbf{x}\\\mathbf{u}\end{bmatrix}\in\mathbb{R}^{n+m} \): and similarly define \( \tilde{\mathbf{z}}=\begin{bmatrix}\tilde{\mathbf{x}}\\\tilde{\mathbf{u}}\end{bmatrix} \):.

Theorem (First-order expansion with integral remainder). Let \( \mathbf{F}:\mathbb{R}^{n+m} → \mathbb{R}^q \) be continuously differentiable in a neighborhood of \( \mathbf{z}_0 \). Then for any small \( \tilde{\mathbf{z}} \) in that neighborhood:

\[ \mathbf{F}(\mathbf{z}_0+\tilde{\mathbf{z}}) = \mathbf{F}(\mathbf{z}_0) + \mathbf{J}_F(\mathbf{z}_0)\,\tilde{\mathbf{z}} + \mathbf{r}_F(\tilde{\mathbf{z}}), \]

where \( \mathbf{J}_F(\mathbf{z}_0) \): is the Jacobian at \( \mathbf{z}_0 \), and the remainder is

\[ \mathbf{r}_F(\tilde{\mathbf{z}}) = \int_0^1 \left(\mathbf{J}_F(\mathbf{z}_0+s\tilde{\mathbf{z}})-\mathbf{J}_F(\mathbf{z}_0)\right) \tilde{\mathbf{z}}\; ds. \]

Proof. Define \( \phi(s)=\mathbf{F}(\mathbf{z}_0+s\tilde{\mathbf{z}}) \). By the chain rule, \( \phi'(s)=\mathbf{J}_F(\mathbf{z}_0+s\tilde{\mathbf{z}})\tilde{\mathbf{z}} \). Then

\[ \mathbf{F}(\mathbf{z}_0+\tilde{\mathbf{z}})-\mathbf{F}(\mathbf{z}_0) = \phi(1)-\phi(0) = \int_0^1 \phi'(s)\,ds = \int_0^1 \mathbf{J}_F(\mathbf{z}_0+s\tilde{\mathbf{z}})\tilde{\mathbf{z}}\,ds. \]

Add and subtract \( \mathbf{J}_F(\mathbf{z}_0)\tilde{\mathbf{z}} \) inside the integral to obtain the stated form. □

Corollary (Quadratic error order). If \( \mathbf{J}_F \) is locally Lipschitz near \( \mathbf{z}_0 \), i.e., there exists \( L > 0 \) such that

\[ \|\mathbf{J}_F(\mathbf{z}_a)-\mathbf{J}_F(\mathbf{z}_b)\| \le L\|\mathbf{z}_a-\mathbf{z}_b\| \quad \text{for all }\mathbf{z}_a,\mathbf{z}_b \text{ near }\mathbf{z}_0, \]

then the remainder satisfies \( \|\mathbf{r}_F(\tilde{\mathbf{z}})\| \le \tfrac{L}{2}\|\tilde{\mathbf{z}}\|^2 \):.

Proof. Using the remainder formula and Lipschitz continuity:

\[ \|\mathbf{r}_F(\tilde{\mathbf{z}})\| \le \int_0^1 \|\mathbf{J}_F(\mathbf{z}_0+s\tilde{\mathbf{z}})-\mathbf{J}_F(\mathbf{z}_0)\|\;\|\tilde{\mathbf{z}}\|\; ds \le \int_0^1 (Ls\|\tilde{\mathbf{z}}\|)\;\|\tilde{\mathbf{z}}\|\; ds = \frac{L}{2}\|\tilde{\mathbf{z}}\|^2. \]

□

Apply Section 3 to both \( \mathbf{f} \) and \( \mathbf{h} \). Writing them as functions of \( \mathbf{z}=[\mathbf{x}^\top\; \mathbf{u}^\top]^\top \), their Jacobians partition into blocks:

\[ \mathbf{J}_f(\mathbf{z}_0) = \begin{bmatrix} \dfrac{\partial \mathbf{f}}{\partial \mathbf{x}} & \dfrac{\partial \mathbf{f}}{\partial \mathbf{u}} \end{bmatrix}_{(\mathbf{x}_0,\mathbf{u}_0)} \triangleq \begin{bmatrix}\mathbf{A} & \mathbf{B}\end{bmatrix}, \qquad \mathbf{J}_h(\mathbf{z}_0) = \begin{bmatrix} \dfrac{\partial \mathbf{h}}{\partial \mathbf{x}} & \dfrac{\partial \mathbf{h}}{\partial \mathbf{u}} \end{bmatrix}_{(\mathbf{x}_0,\mathbf{u}_0)} \triangleq \begin{bmatrix}\mathbf{C} & \mathbf{D}\end{bmatrix}. \]

Using \( \mathbf{x}=\mathbf{x}_0+\tilde{\mathbf{x}} \) and \( \mathbf{u}=\mathbf{u}_0+\tilde{\mathbf{u}} \), we obtain:

\[ \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}_0,\mathbf{u}_0) + \mathbf{A}\tilde{\mathbf{x}}+\mathbf{B}\tilde{\mathbf{u}} + \mathbf{r}_f(\tilde{\mathbf{x}},\tilde{\mathbf{u}}), \]

\[ \mathbf{y} = \mathbf{h}(\mathbf{x}_0,\mathbf{u}_0) + \mathbf{C}\tilde{\mathbf{x}}+\mathbf{D}\tilde{\mathbf{u}} + \mathbf{r}_h(\tilde{\mathbf{x}},\tilde{\mathbf{u}}). \]

Because \( \dot{\mathbf{x}}=\dot{\tilde{\mathbf{x}}} \) (the operating point is constant), and \( \mathbf{y}_0=\mathbf{h}(\mathbf{x}_0,\mathbf{u}_0) \), this becomes

\[ \dot{\tilde{\mathbf{x}}} = \mathbf{f}(\mathbf{x}_0,\mathbf{u}_0) + \mathbf{A}\tilde{\mathbf{x}}+\mathbf{B}\tilde{\mathbf{u}} + \mathbf{r}_f(\tilde{\mathbf{x}},\tilde{\mathbf{u}}), \qquad \tilde{\mathbf{y}} = \mathbf{C}\tilde{\mathbf{x}}+\mathbf{D}\tilde{\mathbf{u}} + \mathbf{r}_h(\tilde{\mathbf{x}},\tilde{\mathbf{u}}). \]

If \( (\mathbf{x}_0,\mathbf{u}_0) \) is an equilibrium, then \( \mathbf{f}(\mathbf{x}_0,\mathbf{u}_0)=\mathbf{0} \):. Neglecting higher-order remainders yields the linearized MIMO model:

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

If \( (\mathbf{x}_0,\mathbf{u}_0) \) is not an equilibrium, the linearized dynamics contain an affine constant term \( \mathbf{f}(\mathbf{x}_0,\mathbf{u}_0) \). In practice, steady-state operating points are preferred because they remove this constant term and produce a standard LTI model.

The Jacobian blocks have fixed dimensions: \( \mathbf{A}\in\mathbb{R}^{n\times n} \): state-to-state sensitivity, \( \mathbf{B}\in\mathbb{R}^{n\times m} \): input-to-state sensitivity, \( \mathbf{C}\in\mathbb{R}^{p\times n} \): state-to-output sensitivity, \( \mathbf{D}\in\mathbb{R}^{p\times m} \): direct input-to-output sensitivity.

For MIMO systems, coupling appears as off-diagonal structure: for example, if \( B_{12}\neq 0 \), then input channel 2 influences state 1 to first order. Similarly, \( C_{21}\neq 0 \) means state 1 influences output 2.

A common local question is: “Does input \( u_j \) immediately affect output \( y_i \)?” The first-order answer is determined by the feedthrough element \( D_{ij} \):. If \( D_{ij}=0 \), then any effect of \( u_j \) on \( y_i \) must pass through the state.

When \( \mathbf{D}\neq \mathbf{0} \), the linearized model is sometimes called proper with feedthrough. Physically, \( \mathbf{D}\neq \mathbf{0} \) can arise from sensor/actuator algebraic relations or measurement models that depend directly on inputs.

For the LTI perturbation model, and assuming zero initial perturbation \( \tilde{\mathbf{x}}(0)=\mathbf{0} \):, take Laplace transforms:

\[ s\tilde{\mathbf{X}}(s) = \mathbf{A}\tilde{\mathbf{X}}(s)+\mathbf{B}\tilde{\mathbf{U}}(s). \]

Solve for \( \tilde{\mathbf{X}}(s) \) (where \( \mathbf{I} \) is the \( n\times n \) identity):

\[ (s\mathbf{I}-\mathbf{A})\tilde{\mathbf{X}}(s) = \mathbf{B}\tilde{\mathbf{U}}(s) \quad \Rightarrow \quad \tilde{\mathbf{X}}(s) = (s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}\tilde{\mathbf{U}}(s). \]

Substitute into \( \tilde{\mathbf{Y}}(s)=\mathbf{C}\tilde{\mathbf{X}}(s)+\mathbf{D}\tilde{\mathbf{U}}(s) \):

\[ \tilde{\mathbf{Y}}(s) = \left(\mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}+\mathbf{D}\right)\tilde{\mathbf{U}}(s). \]

Hence the MIMO transfer matrix from \( \tilde{\mathbf{u}} \) to \( \tilde{\mathbf{y}} \) is

\[ \boxed{\mathbf{G}(s) = \mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}+\mathbf{D}\in\mathbb{R}^{p\times m}.} \]

Each element \( G_{ij}(s) \) is the (local) transfer function from input \( \tilde{u}_j \) to output \( \tilde{y}_i \), with all other inputs held at zero perturbation.

For models derived from physics, you often have explicit formulas for \( \mathbf{f} \) and \( \mathbf{h} \). Then the Jacobians can be computed:

• Symbolically (exact Jacobians) using CAS tools: SymPy (Python), Symbolic Math Toolbox (MATLAB), Mathematica.
• Numerically (approximate Jacobians) using finite differences in C++/Java or when formulas are complex.

A standard finite-difference approximation for \( \partial \mathbf{f}/\partial x_k \) uses a small step \( \epsilon \):

\[ \frac{\partial \mathbf{f}}{\partial x_k}(\mathbf{x}_0,\mathbf{u}_0) \approx \frac{\mathbf{f}(\mathbf{x}_0+\epsilon\mathbf{e}_k,\mathbf{u}_0)-\mathbf{f}(\mathbf{x}_0-\epsilon\mathbf{e}_k,\mathbf{u}_0)}{2\epsilon}, \]

where \( \mathbf{e}_k \) is the \( k \)-th unit vector. The same approach applies for \( \partial \mathbf{f}/\partial u_j \), \( \partial \mathbf{h}/\partial x_k \), and \( \partial \mathbf{h}/\partial u_j \).

The centered difference has truncation error of order \( O(\epsilon^2) \) (assuming sufficient smoothness), but excessively small \( \epsilon \) amplifies floating-point roundoff. In practice, \( \epsilon \) is chosen heuristically (e.g., \( 10^{-6} \) to \( 10^{-8} \) in normalized units), and results should be sanity-checked.

We demonstrate a 2-state, 2-input, 2-output nonlinear MIMO system:

\[ \begin{aligned} \dot{x}_1 &= -x_1 + \sin(x_2) + u_1,\\ \dot{x}_2 &= x_1^2 - 2x_2 + u_2,\\ y_1 &= x_1 + 0.2u_1,\\ y_2 &= x_2. \end{aligned} \]

At the equilibrium \( \mathbf{x}_0=[0\;\;0]^\top \), \( \mathbf{u}_0=[0\;\;0]^\top \), compute \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \) using SymPy. We then compare nonlinear and linearized responses for small inputs using SciPy.

import numpy as np
import sympy as sp

# --- Symbolic model definition ---
x1, x2, u1, u2 = sp.symbols('x1 x2 u1 u2', real=True)
x = sp.Matrix([x1, x2])
u = sp.Matrix([u1, u2])

f = sp.Matrix([
    -x1 + sp.sin(x2) + u1,
    x1**2 - 2*x2 + u2
])

h = sp.Matrix([
    x1 + 0.2*u1,
    x2
])

# Jacobian blocks
A_sym = f.jacobian(x)
B_sym = f.jacobian(u)
C_sym = h.jacobian(x)
D_sym = h.jacobian(u)

# Operating point (equilibrium)
op = {x1: 0.0, x2: 0.0, u1: 0.0, u2: 0.0}

A = np.array(A_sym.subs(op), dtype=float)
B = np.array(B_sym.subs(op), dtype=float)
C = np.array(C_sym.subs(op), dtype=float)
D = np.array(D_sym.subs(op), dtype=float)

print("A =\n", A)
print("B =\n", B)
print("C =\n", C)
print("D =\n", D)

# --- Build linear state-space model and simulate ---
# Libraries commonly used for system dynamics in Python:
# numpy, scipy.integrate for ODEs, sympy for symbolic Jacobians, and python-control for LTI models.
from scipy.integrate import solve_ivp

# Try to use python-control if available; otherwise fall back to direct simulation
try:
    import control

    sys = control.ss(A, B, C, D)

    # Define a small input (piecewise constant) for local validation
    def u_t(t):
        return np.array([0.05*(t >= 0.5), -0.04*(t >= 1.0)], dtype=float)

    # Nonlinear simulation
    def f_nl(t, xvec):
        uu = u_t(t)
        return np.array([
            -xvec[0] + np.sin(xvec[1]) + uu[0],
            xvec[0]**2 - 2*xvec[1] + uu[1]
        ], dtype=float)

    t_span = (0.0, 5.0)
    t_eval = np.linspace(t_span[0], t_span[1], 1001)
    x0 = np.zeros(2)

    sol_nl = solve_ivp(f_nl, t_span, x0, t_eval=t_eval, rtol=1e-9, atol=1e-12)

    # Linear simulation: xdot = A x + B u (perturbation variables)
    # We'll simulate via explicit integration to keep dependencies minimal.
    def f_lin(t, xvec):
        uu = u_t(t)
        return (A @ xvec + B @ uu)

    sol_lin = solve_ivp(f_lin, t_span, x0, t_eval=t_eval, rtol=1e-12, atol=1e-14)

    # Outputs
    y_nl = np.vstack([
        sol_nl.y[0] + 0.2*np.array([u_t(t)[0] for t in sol_nl.t]),
        sol_nl.y[1]
    ])

    y_lin = np.vstack([
        (C @ sol_lin.y + D @ np.vstack([u_t(t) for t in sol_lin.t]).T)[0, :],
        (C @ sol_lin.y + D @ np.vstack([u_t(t) for t in sol_lin.t]).T)[1, :]
    ])

    # Report a simple local-error metric
    err = np.max(np.abs(y_nl - y_lin))
    print("Max |y_nonlinear - y_linear| over time =", err)

except Exception as e:
    print("python-control not available or another error occurred:", e)
      

Interpretation at the operating point: \( \mathbf{A} \): contains the coupling term \( \cos(x_2) \) evaluated at \( x_2=0 \), hence \( A_{12}=1 \). This indicates that small changes in \( x_2 \) influence \( \dot{x}_1 \) linearly near the equilibrium. \( \mathbf{D} \): shows a direct path from \( u_1 \) to \( y_1 \).

When symbolic expressions are unavailable (e.g., legacy code, lookup tables), you can compute \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \) numerically using centered finite differences. Below are compact reference implementations:

9.1 C++ (Eigen + RK4)

Common system-dynamics-related libraries in C++ include: Eigen (linear algebra), Boost::odeint (ODE integration), and custom solvers. The code below uses Eigen and a simple RK4 integrator.

#include <iostream>
#include <cmath>
#include <Eigen/Dense>

// Nonlinear MIMO model (same as Python example)
// xdot = f(x,u), y = h(x,u)
static Eigen::VectorXd f_nl(const Eigen::VectorXd& x, const Eigen::VectorXd& u) {
  Eigen::VectorXd dx(2);
  dx(0) = -x(0) + std::sin(x(1)) + u(0);
  dx(1) = x(0)*x(0) - 2.0*x(1) + u(1);
  return dx;
}

static Eigen::VectorXd h_nl(const Eigen::VectorXd& x, const Eigen::VectorXd& u) {
  Eigen::VectorXd y(2);
  y(0) = x(0) + 0.2*u(0);
  y(1) = x(1);
  return y;
}

// Centered finite-difference Jacobian wrt x
static Eigen::MatrixXd jacobian_x(
    const Eigen::VectorXd& x0, const Eigen::VectorXd& u0, double eps,
    Eigen::VectorXd (*F)(const Eigen::VectorXd&, const Eigen::VectorXd&)) {

  const int n = x0.size();
  const int q = F(x0,u0).size();
  Eigen::MatrixXd J(q, n);

  for (int k = 0; k < n; ++k) {
    Eigen::VectorXd e = Eigen::VectorXd::Zero(n);
    e(k) = 1.0;
    Eigen::VectorXd Fp = F(x0 + eps*e, u0);
    Eigen::VectorXd Fm = F(x0 - eps*e, u0);
    J.col(k) = (Fp - Fm) / (2.0*eps);
  }
  return J;
}

// Centered finite-difference Jacobian wrt u
static Eigen::MatrixXd jacobian_u(
    const Eigen::VectorXd& x0, const Eigen::VectorXd& u0, double eps,
    Eigen::VectorXd (*F)(const Eigen::VectorXd&, const Eigen::VectorXd&)) {

  const int m = u0.size();
  const int q = F(x0,u0).size();
  Eigen::MatrixXd J(q, m);

  for (int j = 0; j < m; ++j) {
    Eigen::VectorXd e = Eigen::VectorXd::Zero(m);
    e(j) = 1.0;
    Eigen::VectorXd Fp = F(x0, u0 + eps*e);
    Eigen::VectorXd Fm = F(x0, u0 - eps*e);
    J.col(j) = (Fp - Fm) / (2.0*eps);
  }
  return J;
}

// Simple RK4 integrator step: x_{k+1} = x_k + (dt/6)(k1+2k2+2k3+k4)
static Eigen::VectorXd rk4_step(
    const Eigen::VectorXd& x, const Eigen::VectorXd& u, double dt,
    Eigen::VectorXd (*f)(const Eigen::VectorXd&, const Eigen::VectorXd&)) {

  Eigen::VectorXd k1 = f(x, u);
  Eigen::VectorXd k2 = f(x + 0.5*dt*k1, u);
  Eigen::VectorXd k3 = f(x + 0.5*dt*k2, u);
  Eigen::VectorXd k4 = f(x + dt*k3, u);
  return x + (dt/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4);
}

int main() {
  Eigen::VectorXd x0(2); x0 << 0.0, 0.0;
  Eigen::VectorXd u0(2); u0 << 0.0, 0.0;
  double eps = 1e-7;

  Eigen::MatrixXd A = jacobian_x(x0, u0, eps, f_nl);
  Eigen::MatrixXd B = jacobian_u(x0, u0, eps, f_nl);
  Eigen::MatrixXd C = jacobian_x(x0, u0, eps, h_nl);
  Eigen::MatrixXd D = jacobian_u(x0, u0, eps, h_nl);

  std::cout << "A=\n" << A << "\n\n";
  std::cout << "B=\n" << B << "\n\n";
  std::cout << "C=\n" << C << "\n\n";
  std::cout << "D=\n" << D << "\n\n";

  // Local validation: simulate nonlinear and linearized (perturbation) models
  auto u_of_t = [](double t) {
    Eigen::VectorXd u(2);
    u(0) = (t >= 0.5) ? 0.05 : 0.0;
    u(1) = (t >= 1.0) ? -0.04 : 0.0;
    return u;
  };

  double dt = 0.001;
  double T  = 5.0;
  int N = static_cast<int>(T/dt);

  Eigen::VectorXd x_nl = x0;
  Eigen::VectorXd x_lin = Eigen::VectorXd::Zero(2);

  double max_err = 0.0;
  for (int k = 0; k <= N; ++k) {
    double t = k*dt;
    Eigen::VectorXd u = u_of_t(t);

    // Nonlinear state (absolute variables)
    x_nl = rk4_step(x_nl, u, dt, f_nl);

    // Linearized perturbation: x_tilde_dot = A x_tilde + B u_tilde (here u_tilde=u)
    auto f_lin = [&](const Eigen::VectorXd& xt) {
      return A*xt + B*u;
    };

    // RK4 for linear model
    Eigen::VectorXd k1 = f_lin(x_lin);
    Eigen::VectorXd k2 = f_lin(x_lin + 0.5*dt*k1);
    Eigen::VectorXd k3 = f_lin(x_lin + 0.5*dt*k2);
    Eigen::VectorXd k4 = f_lin(x_lin + dt*k3);
    x_lin = x_lin + (dt/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4);

    // Outputs
    Eigen::VectorXd y_nl = h_nl(x_nl, u);
    Eigen::VectorXd y_lin = C*x_lin + D*u;

    max_err = std::max(max_err, (y_nl - y_lin).cwiseAbs().maxCoeff());
  }

  std::cout << "Max |y_nonlinear - y_linear| = " << max_err << "\n";
  return 0;
}
      

9.2 Java (EJML for Matrices)

In Java, common choices are EJML (linear algebra) and Apache Commons Math (ODE solvers). Below is a compact EJML-based Jacobian computation and RK4 simulation skeleton.

import org.ejml.data.DMatrixRMaj;
import org.ejml.dense.row.CommonOps_DDRM;

public class MimoLinearization {

  // Nonlinear dynamics f(x,u) for 2-state, 2-input example
  static double[] f_nl(double[] x, double[] u) {
    return new double[] {
      -x[0] + Math.sin(x[1]) + u[0],
      x[0]*x[0] - 2.0*x[1] + u[1]
    };
  }

  // Output map h(x,u)
  static double[] h_nl(double[] x, double[] u) {
    return new double[] {
      x[0] + 0.2*u[0],
      x[1]
    };
  }

  // Centered finite-difference Jacobian wrt x
  static DMatrixRMaj jacobianX(double[] x0, double[] u0, double eps, boolean dynamics) {
    int n = x0.length;
    double[] F0 = dynamics ? f_nl(x0,u0) : h_nl(x0,u0);
    int q = F0.length;
    DMatrixRMaj J = new DMatrixRMaj(q, n);

    for (int k = 0; k < n; k++) {
      double[] xp = x0.clone(); xp[k] += eps;
      double[] xm = x0.clone(); xm[k] -= eps;

      double[] Fp = dynamics ? f_nl(xp,u0) : h_nl(xp,u0);
      double[] Fm = dynamics ? f_nl(xm,u0) : h_nl(xm,u0);

      for (int i = 0; i < q; i++) {
        J.set(i, k, (Fp[i]-Fm[i])/(2.0*eps));
      }
    }
    return J;
  }

  // Centered finite-difference Jacobian wrt u
  static DMatrixRMaj jacobianU(double[] x0, double[] u0, double eps, boolean dynamics) {
    int m = u0.length;
    double[] F0 = dynamics ? f_nl(x0,u0) : h_nl(x0,u0);
    int q = F0.length;
    DMatrixRMaj J = new DMatrixRMaj(q, m);

    for (int j = 0; j < m; j++) {
      double[] up = u0.clone(); up[j] += eps;
      double[] um = u0.clone(); um[j] -= eps;

      double[] Fp = dynamics ? f_nl(x0,up) : h_nl(x0,up);
      double[] Fm = dynamics ? f_nl(x0,um) : h_nl(x0,um);

      for (int i = 0; i < q; i++) {
        J.set(i, j, (Fp[i]-Fm[i])/(2.0*eps));
      }
    }
    return J;
  }

  // RK4 step for xdot = f(x,u)
  static double[] rk4Step(double[] x, double[] u, double dt) {
    double[] k1 = f_nl(x,u);
    double[] x2 = new double[] { x[0] + 0.5*dt*k1[0], x[1] + 0.5*dt*k1[1] };
    double[] k2 = f_nl(x2,u);
    double[] x3 = new double[] { x[0] + 0.5*dt*k2[0], x[1] + 0.5*dt*k2[1] };
    double[] k3 = f_nl(x3,u);
    double[] x4 = new double[] { x[0] + dt*k3[0], x[1] + dt*k3[1] };
    double[] k4 = f_nl(x4,u);

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

  public static void main(String[] args) {
    double[] x0 = new double[] {0.0, 0.0};
    double[] u0 = new double[] {0.0, 0.0};
    double eps = 1e-7;

    DMatrixRMaj A = jacobianX(x0,u0,eps,true);
    DMatrixRMaj B = jacobianU(x0,u0,eps,true);
    DMatrixRMaj C = jacobianX(x0,u0,eps,false);
    DMatrixRMaj D = jacobianU(x0,u0,eps,false);

    System.out.println("A=\n" + A);
    System.out.println("B=\n" + B);
    System.out.println("C=\n" + C);
    System.out.println("D=\n" + D);

    // Simulate linear perturbation: xdot = A x + B u
    double dt = 0.001;
    double T = 5.0;
    int N = (int)(T/dt);

    double[] xLin = new double[] {0.0, 0.0};
    double[] xNl  = x0.clone();

    double maxErr = 0.0;

    for (int k = 0; k <= N; k++) {
      double t = k*dt;
      double[] u = new double[] {
        (t >= 0.5) ? 0.05 : 0.0,
        (t >= 1.0) ? -0.04 : 0.0
      };

      // Nonlinear step
      xNl = rk4Step(xNl, u, dt);

      // Linear step (RK4 with matrix ops)
      // k1 = A x + B u, etc.
      double[] k1 = linDyn(A,B,xLin,u);
      double[] x2 = new double[] { xLin[0] + 0.5*dt*k1[0], xLin[1] + 0.5*dt*k1[1] };
      double[] k2 = linDyn(A,B,x2,u);
      double[] x3 = new double[] { xLin[0] + 0.5*dt*k2[0], xLin[1] + 0.5*dt*k2[1] };
      double[] k3 = linDyn(A,B,x3,u);
      double[] x4 = new double[] { xLin[0] + dt*k3[0], xLin[1] + dt*k3[1] };
      double[] k4 = linDyn(A,B,x4,u);

      xLin[0] += (dt/6.0)*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0]);
      xLin[1] += (dt/6.0)*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1]);

      double[] yNl = h_nl(xNl, u);
      double[] yLin = linOut(C,D,xLin,u);

      maxErr = Math.max(maxErr, Math.max(Math.abs(yNl[0]-yLin[0]), Math.abs(yNl[1]-yLin[1])));
    }

    System.out.println("Max |y_nonlinear - y_linear| = " + maxErr);
  }

  static double[] linDyn(DMatrixRMaj A, DMatrixRMaj B, double[] x, double[] u) {
    // dx = A x + B u
    double dx0 = A.get(0,0)*x[0] + A.get(0,1)*x[1] + B.get(0,0)*u[0] + B.get(0,1)*u[1];
    double dx1 = A.get(1,0)*x[0] + A.get(1,1)*x[1] + B.get(1,0)*u[0] + B.get(1,1)*u[1];
    return new double[] {dx0, dx1};
  }

  static double[] linOut(DMatrixRMaj C, DMatrixRMaj D, double[] x, double[] u) {
    // y = C x + D u
    double y0 = C.get(0,0)*x[0] + C.get(0,1)*x[1] + D.get(0,0)*u[0] + D.get(0,1)*u[1];
    double y1 = C.get(1,0)*x[0] + C.get(1,1)*x[1] + D.get(1,0)*u[0] + D.get(1,1)*u[1];
    return new double[] {y0, y1};
  }
}
      

MATLAB provides symbolic Jacobians (Symbolic Math Toolbox) and standard state-space tools (Control System Toolbox). Simulink additionally supports extracting linear models around operating points using linmod (legacy) or linearize (modern workflows; may require Simulink Control Design depending on usage).

10.1 MATLAB symbolic Jacobians

% Define symbolic variables
syms x1 x2 u1 u2 real
x = [x1; x2];
u = [u1; u2];

% Nonlinear model
f = [-x1 + sin(x2) + u1;
      x1^2 - 2*x2 + u2];

h = [x1 + 0.2*u1;
     x2];

% Jacobians
A_sym = jacobian(f, x);
B_sym = jacobian(f, u);
C_sym = jacobian(h, x);
D_sym = jacobian(h, u);

% Operating point
x0 = [0; 0];
u0 = [0; 0];

A = double(subs(A_sym, [x1 x2 u1 u2], [x0.' u0.']));
B = double(subs(B_sym, [x1 x2 u1 u2], [x0.' u0.']));
C = double(subs(C_sym, [x1 x2 u1 u2], [x0.' u0.']));
D = double(subs(D_sym, [x1 x2 u1 u2], [x0.' u0.']));

% Linearized perturbation model
sys = ss(A,B,C,D);

disp(A); disp(B); disp(C); disp(D);
      

10.2 Simulink workflow (model-based linearization)

In Simulink, you typically (i) build the nonlinear block model, (ii) specify the operating point, then (iii) extract \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \) and export an ss object for analysis.

flowchart TD
  M1["Build block model of f,h"] --> OP["Set operating point (x0,u0)"]
  OP --> TR["Select linearization I/O points"]
  TR --> LIN["Compute linear model (A,B,C,D)"]
  LIN --> EXP["Export as state-space object"]
  EXP --> ANA["Analyze with LTI tools"]
  ANA --> VER["Validate locally by simulation"]
        

A minimal script pattern (exact commands depend on your Simulink configuration and model name) is:

% Suppose the Simulink model is named: 'mimo_nonlinear_model'
% and is already built and saved.

model = 'mimo_nonlinear_model';
load_system(model);

% Option A (legacy): linmod (works for many continuous models)
% [A,B,C,D] = linmod(model, x0, u0);
% sys = ss(A,B,C,D);

% Option B (operating point + linearize):
% op = operpoint(model);           % or create using findop
% sys = linearize(model, op);      % returns ss model

% After obtaining sys:
% step(sys) or lsim(sys, u, t) can be used for LTI analysis.
      

Recommendation: always compare nonlinear vs linearized outputs for small perturbations to confirm the local validity expected from Lesson 3.

Mathematica can compute Jacobians exactly and manipulate state-space models symbolically.

(* Define variables *)
x = {x1, x2};
u = {u1, u2};

(* Nonlinear model *)
f = {
  -x1 + Sin[x2] + u1,
  x1^2 - 2 x2 + u2
};

h = {
  x1 + 0.2 u1,
  x2
};

(* Jacobian blocks *)
A = D[f, {x}];
B = D[f, {u}];
C = D[h, {x}];
Dmat = D[h, {u}];

(* Operating point *)
op = {x1 -> 0, x2 -> 0, u1 -> 0, u2 -> 0};

A0 = A /. op // N;
B0 = B /. op // N;
C0 = C /. op // N;
D0 = Dmat /. op // N;

A0
B0
C0
D0

(* Create linear state-space model and compute transfer matrix *)
sys = StateSpaceModel[{A0, B0, C0, D0}];
G = TransferFunctionModel[sys];
G
      

The resulting transfer matrix \( \mathbf{G}(s) \) is the local MIMO approximation mapping \( \tilde{\mathbf{u}} → \tilde{\mathbf{y}} \) in the Laplace domain.

Problem 1 (Block Jacobians for a MIMO model): Consider

\[ \dot{\mathbf{x}}= \begin{bmatrix} -x_1 + \sin(x_2) + u_1\\ x_1^2 - 2x_2 + u_2 \end{bmatrix}, \qquad \mathbf{y}= \begin{bmatrix} x_1 + 0.2u_1\\ x_2 \end{bmatrix}. \]

Compute \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \) at the equilibrium \( \mathbf{x}_0=\begin{bmatrix}0\\0\end{bmatrix} \), \( \mathbf{u}_0=\begin{bmatrix}0\\0\end{bmatrix} \).

Solution: By definition, \( \mathbf{A}=\left.\dfrac{\partial \mathbf{f}}{\partial \mathbf{x}}\right|_{(\mathbf{x}_0,\mathbf{u}_0)} \), \( \mathbf{B}=\left.\dfrac{\partial \mathbf{f}}{\partial \mathbf{u}}\right|_{(\mathbf{x}_0,\mathbf{u}_0)} \), \( \mathbf{C}=\left.\dfrac{\partial \mathbf{h}}{\partial \mathbf{x}}\right|_{(\mathbf{x}_0,\mathbf{u}_0)} \), \( \mathbf{D}=\left.\dfrac{\partial \mathbf{h}}{\partial \mathbf{u}}\right|_{(\mathbf{x}_0,\mathbf{u}_0)} \).

\[ \frac{\partial \mathbf{f}}{\partial \mathbf{x}} = \begin{bmatrix} \dfrac{\partial}{\partial x_1}(-x_1+\sin x_2+u_1) & \dfrac{\partial}{\partial x_2}(-x_1+\sin x_2+u_1)\\[6pt] \dfrac{\partial}{\partial x_1}(x_1^2-2x_2+u_2) & \dfrac{\partial}{\partial x_2}(x_1^2-2x_2+u_2) \end{bmatrix} = \begin{bmatrix} -1 & \cos(x_2)\\ 2x_1 & -2 \end{bmatrix}. \]

\[ \Rightarrow\quad \mathbf{A} = \begin{bmatrix} -1 & \cos(0)\\ 2\cdot 0 & -2 \end{bmatrix} = \begin{bmatrix} -1 & 1\\ 0 & -2 \end{bmatrix}. \]

\[ \mathbf{B} = \left. \frac{\partial \mathbf{f}}{\partial \mathbf{u}} \right|_{(\mathbf{x}_0,\mathbf{u}_0)} = \begin{bmatrix} \dfrac{\partial}{\partial u_1}(\cdot) & \dfrac{\partial}{\partial u_2}(\cdot)\\ \dfrac{\partial}{\partial u_1}(\cdot) & \dfrac{\partial}{\partial u_2}(\cdot) \end{bmatrix} = \begin{bmatrix} 1 & 0\\ 0 & 1 \end{bmatrix}. \]

\[ \mathbf{C} = \left. \frac{\partial \mathbf{h}}{\partial \mathbf{x}} \right|_{(\mathbf{x}_0,\mathbf{u}_0)} = \begin{bmatrix} 1 & 0\\ 0 & 1 \end{bmatrix}, \qquad \mathbf{D} = \left. \frac{\partial \mathbf{h}}{\partial \mathbf{u}} \right|_{(\mathbf{x}_0,\mathbf{u}_0)} = \begin{bmatrix} 0.2 & 0\\ 0 & 0 \end{bmatrix}. \]

Problem 2 (Why equilibria remove the affine term): Show that if \( \mathbf{f}(\mathbf{x}_0,\mathbf{u}_0)=\mathbf{0} \), then the linearized perturbation dynamics have no constant term.

Solution: From Section 4,

\[ \dot{\tilde{\mathbf{x}}} = \mathbf{f}(\mathbf{x}_0,\mathbf{u}_0)+\mathbf{A}\tilde{\mathbf{x}}+\mathbf{B}\tilde{\mathbf{u}}+\mathbf{r}_f(\tilde{\mathbf{x}},\tilde{\mathbf{u}}). \]

If \( (\mathbf{x}_0,\mathbf{u}_0) \) is an equilibrium, then \( \mathbf{f}(\mathbf{x}_0,\mathbf{u}_0)=\mathbf{0} \). Neglecting higher-order terms yields

\[ \dot{\tilde{\mathbf{x}}}=\mathbf{A}\tilde{\mathbf{x}}+\mathbf{B}\tilde{\mathbf{u}}, \]

which has no affine offset. □

Problem 3 (Transfer matrix derivation): Starting from \( \dot{\tilde{\mathbf{x}}}=\mathbf{A}\tilde{\mathbf{x}}+\mathbf{B}\tilde{\mathbf{u}} \), \( \tilde{\mathbf{y}}=\mathbf{C}\tilde{\mathbf{x}}+\mathbf{D}\tilde{\mathbf{u}} \), and \( \tilde{\mathbf{x}}(0)=\mathbf{0} \), derive \( \tilde{\mathbf{Y}}(s)=\mathbf{G}(s)\tilde{\mathbf{U}}(s) \) and provide \( \mathbf{G}(s) \).

Solution: Taking Laplace transforms gives

\[ s\tilde{\mathbf{X}}(s)=\mathbf{A}\tilde{\mathbf{X}}(s)+\mathbf{B}\tilde{\mathbf{U}}(s) \Rightarrow (s\mathbf{I}-\mathbf{A})\tilde{\mathbf{X}}(s)=\mathbf{B}\tilde{\mathbf{U}}(s). \]

\[ \tilde{\mathbf{X}}(s)=(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}\tilde{\mathbf{U}}(s). \]

Substitute into the output equation:

\[ \tilde{\mathbf{Y}}(s) = \left(\mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}+\mathbf{D}\right)\tilde{\mathbf{U}}(s), \]

hence \( \mathbf{G}(s)=\mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}+\mathbf{D} \). □

Problem 4 (Coordinate change and linearization): Let \( \tilde{\mathbf{x}}=\mathbf{T}\tilde{\mathbf{z}} \) where \( \mathbf{T}\in\mathbb{R}^{n\times n} \) is invertible. Show that the linearized matrices transform as \( \mathbf{A}_z=\mathbf{T}^{-1}\mathbf{A}\mathbf{T} \) and \( \mathbf{B}_z=\mathbf{T}^{-1}\mathbf{B} \), while \( \mathbf{C}_z=\mathbf{C}\mathbf{T} \) and \( \mathbf{D}_z=\mathbf{D} \).

Solution: Substitute \( \tilde{\mathbf{x}}=\mathbf{T}\tilde{\mathbf{z}} \) into \( \dot{\tilde{\mathbf{x}}}=\mathbf{A}\tilde{\mathbf{x}}+\mathbf{B}\tilde{\mathbf{u}} \):

\[ \mathbf{T}\dot{\tilde{\mathbf{z}}}=\mathbf{A}\mathbf{T}\tilde{\mathbf{z}}+\mathbf{B}\tilde{\mathbf{u}} \Rightarrow \dot{\tilde{\mathbf{z}}}=\mathbf{T}^{-1}\mathbf{A}\mathbf{T}\tilde{\mathbf{z}}+\mathbf{T}^{-1}\mathbf{B}\tilde{\mathbf{u}}. \]

So \( \mathbf{A}_z=\mathbf{T}^{-1}\mathbf{A}\mathbf{T} \) and \( \mathbf{B}_z=\mathbf{T}^{-1}\mathbf{B} \). For outputs, \( \tilde{\mathbf{y}}=\mathbf{C}\tilde{\mathbf{x}}+\mathbf{D}\tilde{\mathbf{u}}= \mathbf{C}\mathbf{T}\tilde{\mathbf{z}}+\mathbf{D}\tilde{\mathbf{u}} \), hence \( \mathbf{C}_z=\mathbf{C}\mathbf{T} \), \( \mathbf{D}_z=\mathbf{D} \). □

Problem 5 (Local error order): Using the remainder bound from Section 3, argue that for sufficiently small perturbations, the modeling error between the nonlinear model and the linearized model scales like \( O(\|\tilde{\mathbf{z}}\|^2) \).

Solution: From Section 3,

\[ \mathbf{f}(\mathbf{z}_0+\tilde{\mathbf{z}}) = \mathbf{f}(\mathbf{z}_0)+\mathbf{J}_f(\mathbf{z}_0)\tilde{\mathbf{z}}+\mathbf{r}_f(\tilde{\mathbf{z}}), \]

and if \( \mathbf{J}_f \) is locally Lipschitz with constant \( L \), then \( \|\mathbf{r}_f(\tilde{\mathbf{z}})\|\le \tfrac{L}{2}\|\tilde{\mathbf{z}}\|^2 \). Neglecting \( \mathbf{r}_f \) is therefore a second-order approximation, meaning the first-order linear model error decreases quadratically as perturbations shrink. □

We extended Jacobian linearization to MIMO nonlinear systems by (i) framing the system as maps \( \mathbf{f}(\mathbf{x},\mathbf{u}) \) and \( \mathbf{h}(\mathbf{x},\mathbf{u}) \), (ii) using perturbation variables around an operating point, (iii) deriving the small-signal model with Jacobian blocks \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \), and (iv) connecting the result to the transfer matrix \( \mathbf{G}(s)=\mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}+\mathbf{D} \). We also provided rigorous remainder bounds and practical computation routes (symbolic and finite-difference) across multiple programming environments.