In Chapter 4, we introduced the continuous-time LTI state-space model

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

In the multi-input, multi-output (MIMO) case, we track dimensions explicitly: \( \mathbf{x}(t)\in\mathbb{R}^n \), \( \mathbf{u}(t)\in\mathbb{R}^m \), \( \mathbf{y}(t)\in\mathbb{R}^p \), with matrices \( \mathbf{A}\in\mathbb{R}^{n\times n} \), \( \mathbf{B}\in\mathbb{R}^{n\times m} \), \( \mathbf{C}\in\mathbb{R}^{p\times n} \), \( \mathbf{D}\in\mathbb{R}^{p\times m} \).

Compared with SISO conversion, the key structural differences are:

• Multiple inputs: the input channel becomes a vector and \( \mathbf{B} \) gains multiple columns—one per input.
• Multiple outputs: outputs become a vector; \( \mathbf{C} \) stacks output-selection/combination rows; \( \mathbf{D} \) captures any direct feedthrough from inputs.
• Coupling is allowed: outputs can depend on the same state components, and each input can drive multiple states.

flowchart TD
  U["Input vector u(t) (m channels)"] --> B["State dynamics: x_dot = A x + B u"]
  B --> X["State vector x(t) (n states)"]
  X --> C["Output map: y = C x + D u"]
  U --> C
  C --> Y["Output vector y(t) (p channels)"]
        

In this chapter, we focus on constructing \( (\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}) \) from given higher-order ODE descriptions. No frequency-domain tools are required.

Consider an n-th order scalar input–output differential equation (single output signal \( y \)), but driven by m inputs collected in \( \mathbf{u}(t)=[u_1(t)\;\;u_2(t)\;\;\dots\;\;u_m(t)]^\top \):

\[ y^{(n)}(t) + a_{n-1}y^{(n-1)}(t) + \cdots + a_1\dot{y}(t) + a_0 y(t) = \mathbf{b}^\top \mathbf{u}(t), \]

where \( \mathbf{b}\in\mathbb{R}^m \) is the vector of input gains. Following Lessons 1–2 (phase-variable/companion construction), define the state:

\[ \mathbf{x}(t) = \begin{bmatrix} x_1(t)\\x_2(t)\\ \vdots\\ x_n(t) \end{bmatrix} := \begin{bmatrix} y(t)\\ \dot{y}(t)\\ \vdots\\ y^{(n-1)}(t) \end{bmatrix}. \]

Then, by definition of the state components,

\[ \dot{x}_1 = x_2,\;\dot{x}_2=x_3,\;\dots,\;\dot{x}_{n-1}=x_n. \]

The last equation comes from the given ODE:

\[ \dot{x}_n = y^{(n)}(t) = -a_0 x_1 - a_1 x_2 - \cdots - a_{n-1} x_n + \mathbf{b}^\top \mathbf{u}(t). \]

Collecting these into matrix form yields the companion \( \mathbf{A} \) and a multi-column \( \mathbf{B} \):

\[ \dot{\mathbf{x}}(t) = \underbrace{\begin{bmatrix} 0 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & 1\\ -a_0 & -a_1 & -a_2 & \cdots & -a_{n-1} \end{bmatrix}}_{\mathbf{A}\in\mathbb{R}^{n\times n}} \mathbf{x}(t) + \underbrace{\begin{bmatrix} 0 & 0 & \cdots & 0\\ 0 & 0 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & 0\\ b_1 & b_2 & \cdots & b_m \end{bmatrix}}_{\mathbf{B}\in\mathbb{R}^{n\times m}} \mathbf{u}(t). \]

Here the last row of \( \mathbf{B} \) is \( \mathbf{b}^\top \), and all other entries are zero, reflecting that inputs enter only the highest derivative equation in this canonical construction.

2.1 Output equation for the scalar-output case

If the measured output is exactly \( y(t) \), then \( y(t)=x_1(t) \) and we choose:

\[ y(t) = \underbrace{\begin{bmatrix}1 & 0 & \cdots & 0\end{bmatrix}}_{\mathbf{C}\in\mathbb{R}^{1\times n}} \mathbf{x}(t) + \underbrace{\begin{bmatrix}0 & 0 & \cdots & 0\end{bmatrix}}_{\mathbf{D}\in\mathbb{R}^{1\times m}} \mathbf{u}(t). \]

2.2 Equivalence proof (input–output ODE ↔ state-space)

We now prove that the state-space model above is equivalent to the original ODE, in the precise sense that any sufficiently differentiable solution \( y(t) \) of the ODE generates a state \( \mathbf{x}(t) \) that satisfies the state equations, and conversely any solution of the state equations yields a \( y(t)=x_1(t) \) satisfying the ODE.

Proposition (Equivalence in the M-input scalar-output canonical construction): Assume \( \mathbf{u}(t) \) is continuous and \( y(t) \) is \( n \)-times differentiable. Define \( \mathbf{x}(t)=[y,\dot{y},\dots,y^{(n-1)}]^\top \). Then \( y(t) \) satisfies the ODE if and only if \( (\mathbf{x}(t),\mathbf{u}(t)) \) satisfies \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \) and \( y=\mathbf{C}\mathbf{x}+\mathbf{D}\mathbf{u} \).

Proof:

• (ODE ⇒ state-space) If \( y(t) \) satisfies the ODE, then by the state definition we have \( \dot{x}_i = x_{i+1} \) for \( i=1,\dots,n-1 \). Moreover, the ODE implies \( y^{(n)} = -\sum_{k=0}^{n-1} a_k y^{(k)} + \mathbf{b}^\top \mathbf{u} \), which is exactly \( \dot{x}_n = -a_0 x_1 - \cdots - a_{n-1} x_n + \mathbf{b}^\top \mathbf{u} \). Stacking these gives \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \), and \( y=x_1=\mathbf{C}\mathbf{x} \).
• (state-space ⇒ ODE) If \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \) holds, the first \( n-1 \) rows force \( \dot{x}_1=x_2,\dots,\dot{x}_{n-1}=x_n \). Therefore, \( x_k = y^{(k-1)} \) for all \( k=1,\dots,n \) (by repeated differentiation). The last row then becomes \( \dot{x}_n = y^{(n)} = -a_0 y - a_1 \dot{y} - \cdots - a_{n-1} y^{(n-1)} + \mathbf{b}^\top \mathbf{u} \), which is exactly the original ODE. ∎

The proof is identical in structure to the SISO case, but note that the forcing term is now the scalar \( \mathbf{b}^\top \mathbf{u}(t) \), i.e., a linear combination of the input channels.

Suppose the internal state vector \( \mathbf{x}(t)\in\mathbb{R}^n \) has already been constructed from the underlying ODEs (Lessons 1–2), and we now specify p measured outputs \( \mathbf{y}(t)\in\mathbb{R}^p \).

In this chapter, we restrict to output equations that are algebraic linear combinations of the state and the (possibly direct) inputs:

\[ \mathbf{y}(t) = \mathbf{C}\mathbf{x}(t) + \mathbf{D}\mathbf{u}(t), \qquad \mathbf{C}\in\mathbb{R}^{p\times n},\;\mathbf{D}\in\mathbb{R}^{p\times m}. \]

3.1 Common output patterns in ODE-derived models

• Output is a subset of states: e.g., measure position and velocity, so \( y_1=x_1 \), \( y_2=x_2 \). Then \( \mathbf{C} \) is a row-selector matrix and typically \( \mathbf{D}=\mathbf{0} \).
• Output is a linear combination: e.g., \( y_1 = x_1 + 2x_3 \). Then the corresponding row of \( \mathbf{C} \) contains these coefficients.
• Direct feedthrough (algebraic output dependence on input): e.g., \( y_2 = x_2 + 0.1 u_1 \). Then \( \mathbf{D} \neq \mathbf{0} \).

3.2 Stacking outputs to build \( \mathbf{C} \) and \( \mathbf{D} \)

If each output is given as \( y_i(t) = \mathbf{c}_i^\top \mathbf{x}(t) + \mathbf{d}_i^\top \mathbf{u}(t) \) for \( i=1,\dots,p \), then stacking them yields:

\[ \mathbf{C} = \begin{bmatrix} \mathbf{c}_1^\top\\ \mathbf{c}_2^\top\\ \vdots\\ \mathbf{c}_p^\top \end{bmatrix}, \qquad \mathbf{D} = \begin{bmatrix} \mathbf{d}_1^\top\\ \mathbf{d}_2^\top\\ \vdots\\ \mathbf{d}_p^\top \end{bmatrix}. \]

3.3 Output differentiation is not required here

You may be tempted to treat outputs like \( \dot{y} \) as “new equations.” In this chapter, we do not create new states from output definitions. If you want \( \dot{y} \) as an output and \( \dot{y} \) is already a state component (as in phase-variable form), it is simply a row selection in \( \mathbf{C} \).

A typical “multiple-output” ODE description provides multiple differential equations, e.g., two outputs \( y_1(t),y_2(t) \) each governed by its own higher-order equation, possibly driven by the same input vector \( \mathbf{u}(t) \). A principled (and chapter-consistent) construction is to form an augmented state by stacking the phase-variable states of each output equation.

Example structure (two scalar ODEs of order \( n_1 \) and \( n_2 \)):

\[ \begin{aligned} y_1^{(n_1)} + a_{1,n_1-1}y_1^{(n_1-1)} + \cdots + a_{1,0}y_1 &= \mathbf{b}_1^\top \mathbf{u},\\ y_2^{(n_2)} + a_{2,n_2-1}y_2^{(n_2-1)} + \cdots + a_{2,0}y_2 &= \mathbf{b}_2^\top \mathbf{u}. \end{aligned} \]

Define states for each output separately: \( \mathbf{x}^{(1)}=[y_1,\dot{y}_1,\dots,y_1^{(n_1-1)}]^\top \in\mathbb{R}^{n_1} \), \( \mathbf{x}^{(2)}=[y_2,\dot{y}_2,\dots,y_2^{(n_2-1)}]^\top \in\mathbb{R}^{n_2} \), and form the augmented state:

\[ \mathbf{x} = \begin{bmatrix} \mathbf{x}^{(1)}\\ \mathbf{x}^{(2)} \end{bmatrix} \in \mathbb{R}^{n_1+n_2}. \]

In the absence of explicit coupling terms between \( y_1 \) and \( y_2 \) in the ODEs, the resulting state matrix is block diagonal:

\[ \mathbf{A} = \begin{bmatrix} \mathbf{A}_1 & \mathbf{0}\\ \mathbf{0} & \mathbf{A}_2 \end{bmatrix}, \qquad \mathbf{B} = \begin{bmatrix} \mathbf{B}_1\\ \mathbf{B}_2 \end{bmatrix}, \]

where \( (\mathbf{A}_i,\mathbf{B}_i) \) are companion-form blocks like in Section 2, with the last row of \( \mathbf{B}_i \) equal to \( \mathbf{b}_i^\top \).

Outputs \( \mathbf{y}=[y_1\;\;y_2]^\top \) are extracted by selecting the first state in each block:

\[ \mathbf{y} = \underbrace{ \begin{bmatrix} 1 & 0 & \cdots & 0 & 0 & \cdots & 0\\ 0 & \cdots & 0 & 1 & 0 & \cdots & 0 \end{bmatrix}}_{\mathbf{C}\in\mathbb{R}^{2\times(n_1+n_2)}} \mathbf{x} + \underbrace{\mathbf{0}}_{\mathbf{D}\in\mathbb{R}^{2\times m}}\mathbf{u}. \]

Remark (coupled ODEs): if the ODEs contain cross-terms such as \( y_1 \) depending on \( y_2 \) or its derivatives, then block diagonality is lost; cross-terms appear as off-diagonal blocks in \( \mathbf{A} \). The construction is still systematic: every term that is linear in the chosen states becomes an entry in \( \mathbf{A} \), while input terms become entries in \( \mathbf{B} \).

flowchart TD
  S0["Start: given ODEs + chosen outputs"] --> S1["Choose states per ODE: x = [y, y_dot, ..., y^(n-1)]"]
  S1 --> S2["Write x_dot equations; isolate highest derivatives"]
  S2 --> S3["Assemble A from state-to-state coefficients"]
  S2 --> S4["Assemble B from input coefficients (m columns)"]
  S3 --> S5["Define outputs: y = C x + D u (stack p rows)"]
  S4 --> S5
  S5 --> S6["Check dimensions: A(nxn), B(nxm), C(pxn), D(pxm)"]
  S6 --> S7["Validate by differentiating: recover original ODE(s)"]
        

Consider the scalar 2nd-order ODE driven by two inputs \( u_1,u_2 \):

\[ \ddot{y}(t) + 3\dot{y}(t) + 2y(t) = 4u_1(t) - u_2(t). \]

Define \( x_1=y \), \( x_2=\dot{y} \). Then \( \dot{x}_1=x_2 \) and

\[ \dot{x}_2 = \ddot{y} = -2x_1 - 3x_2 + 4u_1 - u_2. \]

Hence

\[ \mathbf{A}=\begin{bmatrix}0 & 1\\ -2 & -3\end{bmatrix},\quad \mathbf{B}=\begin{bmatrix}0 & 0\\ 4 & -1\end{bmatrix}. \]

Now define two outputs: \( y_1=y=x_1 \) and \( y_2=\dot{y}+0.5u_2=x_2+0.5u_2 \). Then

\[ \mathbf{C}=\begin{bmatrix}1 & 0\\ 0 & 1\end{bmatrix},\quad \mathbf{D}=\begin{bmatrix}0 & 0\\ 0 & 0.5\end{bmatrix}. \]

This example illustrates the central “MIMO bookkeeping” principle: \( \mathbf{B} \) has one column per input channel, while \( \mathbf{C} \) has one row per output channel.

The control (python-control) package provides MIMO state-space objects and simulation utilities. We implement the worked example from Section 5 and simulate a response to piecewise inputs.

import numpy as np

# Example from Section 5
A = np.array([[0.0, 1.0],
              [-2.0, -3.0]])
B = np.array([[0.0, 0.0],
              [4.0, -1.0]])
C = np.array([[1.0, 0.0],
              [0.0, 1.0]])
D = np.array([[0.0, 0.0],
              [0.0, 0.5]])

# Build MIMO state-space model
import control as ctrl
sys = ctrl.ss(A, B, C, D)

# Time grid
t = np.linspace(0.0, 10.0, 2001)

# Inputs: u1 is a step, u2 is a decaying exponential
u1 = np.ones_like(t)
u2 = np.exp(-0.7 * t)

# Stack inputs as shape (m, len(t)) expected by python-control
U = np.vstack([u1, u2])

# Simulate with zero initial state
t_out, y_out, x_out = ctrl.forced_response(sys, T=t, U=U, X0=np.zeros(2), return_x=True)

print("y_out shape:", y_out.shape)  # (p, N)
print("x_out shape:", x_out.shape)  # (n, N)

# Optional: verify that y2 = x2 + 0.5 u2 numerically
y2_check = x_out[1, :] + 0.5 * u2
max_err = np.max(np.abs(y_out[1, :] - y2_check))
print("max|y2 - (x2+0.5u2)| =", max_err)
      

Notes for students:

• In MIMO simulation, the input signal is a matrix whose rows are input channels.
• The output signal is likewise a matrix whose rows are output channels.
• The matrices \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \) must satisfy dimension consistency: \( \mathbf{A}\mathbf{x} \) and \( \mathbf{B}\mathbf{u} \) are both \( n\times 1 \), and \( \mathbf{C}\mathbf{x} \), \( \mathbf{D}\mathbf{u} \) are \( p\times 1 \).

Below is a minimal MIMO state-space simulator using the forward Euler method: \( \mathbf{x}_{k+1}=\mathbf{x}_k + \Delta t(\mathbf{A}\mathbf{x}_k+\mathbf{B}\mathbf{u}_k) \), \( \mathbf{y}_k=\mathbf{C}\mathbf{x}_k+\mathbf{D}\mathbf{u}_k \). This is purely instructional; more accurate integrators will be used later.

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

int main() {
  using Eigen::MatrixXd;
  using Eigen::VectorXd;

  // Matrices from Section 5
  MatrixXd A(2,2), B(2,2), C(2,2), D(2,2);
  A << 0.0, 1.0,
      -2.0, -3.0;
  B << 0.0, 0.0,
       4.0, -1.0;
  C << 1.0, 0.0,
       0.0, 1.0;
  D << 0.0, 0.0,
       0.0, 0.5;

  double T = 10.0;
  double dt = 0.005;
  int N = static_cast<int>(T / dt) + 1;

  VectorXd x = VectorXd::Zero(2);

  for (int k = 0; k < N; ++k) {
    double t = k * dt;

    // Inputs u1(t)=1, u2(t)=exp(-0.7 t)
    VectorXd u(2);
    u(0) = 1.0;
    u(1) = std::exp(-0.7 * t);

    // Output
    VectorXd y = C * x + D * u;

    // Print a few samples
    if (k % 400 == 0) {
      std::cout << "t=" << t
                << "  y1=" << y(0)
                << "  y2=" << y(1)
                << "  x1=" << x(0)
                << "  x2=" << x(1)
                << std::endl;
    }

    // Euler step
    VectorXd xdot = A * x + B * u;
    x = x + dt * xdot;
  }

  return 0;
}
      

The important MIMO detail is that \( \mathbf{u} \) is a vector and \( \mathbf{B} \) is a matrix with one column per input channel.

In Java, EJML provides efficient dense linear algebra. The structure mirrors the C++ example.

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

public class MimoStateSpaceEuler {
  public static void main(String[] args) {
    // Matrices from Section 5
    DMatrixRMaj A = new DMatrixRMaj(new double[][]{
      {0.0, 1.0},
      {-2.0, -3.0}
    });
    DMatrixRMaj B = new DMatrixRMaj(new double[][]{
      {0.0, 0.0},
      {4.0, -1.0}
    });
    DMatrixRMaj C = new DMatrixRMaj(new double[][]{
      {1.0, 0.0},
      {0.0, 1.0}
    });
    DMatrixRMaj D = new DMatrixRMaj(new double[][]{
      {0.0, 0.0},
      {0.0, 0.5}
    });

    double T = 10.0;
    double dt = 0.005;
    int N = (int)(T / dt) + 1;

    DMatrixRMaj x = new DMatrixRMaj(2, 1); // zero init

    DMatrixRMaj u = new DMatrixRMaj(2, 1);
    DMatrixRMaj Ax = new DMatrixRMaj(2, 1);
    DMatrixRMaj Bu = new DMatrixRMaj(2, 1);
    DMatrixRMaj xdot = new DMatrixRMaj(2, 1);

    DMatrixRMaj y = new DMatrixRMaj(2, 1);
    DMatrixRMaj Cx = new DMatrixRMaj(2, 1);
    DMatrixRMaj Du = new DMatrixRMaj(2, 1);

    for (int k = 0; k < N; k++) {
      double t = k * dt;

      // Inputs
      u.set(0, 0, 1.0);
      u.set(1, 0, Math.exp(-0.7 * t));

      // y = Cx + Du
      CommonOps_DDRM.mult(C, x, Cx);
      CommonOps_DDRM.mult(D, u, Du);
      CommonOps_DDRM.add(Cx, Du, y);

      if (k % 400 == 0) {
        System.out.printf("t=%.3f  y1=%.6f  y2=%.6f  x1=%.6f  x2=%.6f%n",
          t, y.get(0,0), y.get(1,0), x.get(0,0), x.get(1,0));
      }

      // xdot = Ax + Bu
      CommonOps_DDRM.mult(A, x, Ax);
      CommonOps_DDRM.mult(B, u, Bu);
      CommonOps_DDRM.add(Ax, Bu, xdot);

      // Euler update: x = x + dt * xdot
      x.add(0, 0, dt * xdot.get(0,0));
      x.add(1, 0, dt * xdot.get(1,0));
    }
  }
}
      

As in C++, the MIMO structure is reflected in the matrix dimensions and channel stacking.

MATLAB’s Control System Toolbox directly supports MIMO state-space via ss.

% Matrices from Section 5
A = [0 1; -2 -3];
B = [0 0; 4 -1];
C = [1 0; 0 1];
D = [0 0; 0 0.5];

sys = ss(A,B,C,D);

t = linspace(0,10,2001);
u1 = ones(size(t));
u2 = exp(-0.7*t);
U = [u1(:) u2(:)];   % N-by-m

x0 = [0; 0];
[y, tOut, x] = lsim(sys, U, t, x0);

% y is N-by-p, x is N-by-n
disp(size(y));
disp(size(x));

% Verify y2 = x2 + 0.5*u2
y2check = x(:,2) + 0.5*u2(:);
fprintf("max|y2 - (x2+0.5u2)| = %.3e\n", max(abs(y(:,2) - y2check)));
      

9.1 Simulink model construction (programmatic)

A state-space block in Simulink can realize the same model. The script below creates a simple model with: (i) two input sources, (ii) a MIMO State-Space block, (iii) scopes for outputs.

modelName = 'mimo_state_space_demo';
new_system(modelName); open_system(modelName);

% Add input sources
add_block('simulink/Sources/Step', [modelName '/u1_step'], 'Position', [50 50 80 80]);
add_block('simulink/Sources/Clock', [modelName '/clock'], 'Position', [50 130 80 160]);
add_block('simulink/Math Operations/Gain', [modelName '/gain_minus0p7'], ...
  'Gain', '-0.7', 'Position', [120 130 170 160]);
add_block('simulink/Math Operations/Math Function', [modelName '/exp'], ...
  'Operator', 'exp', 'Position', [200 130 250 160]);

% Mux inputs to 2-channel vector
add_block('simulink/Signal Routing/Mux', [modelName '/Mux'], 'Inputs', '2', 'Position', [300 70 320 150]);

% Add State-Space block
add_block('simulink/Continuous/State-Space', [modelName '/StateSpace'], 'Position', [380 80 480 140]);

A = [0 1; -2 -3];
B = [0 0; 4 -1];
C = [1 0; 0 1];
D = [0 0; 0 0.5];

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

% Add Scope
add_block('simulink/Sinks/Scope', [modelName '/Scope'], 'Position', [540 90 570 120]);

% Wiring
add_line(modelName, 'u1_step/1', 'Mux/1');
add_line(modelName, 'clock/1', 'gain_minus0p7/1');
add_line(modelName, 'gain_minus0p7/1', 'exp/1');
add_line(modelName, 'exp/1', 'Mux/2');
add_line(modelName, 'Mux/1', 'StateSpace/1');
add_line(modelName, 'StateSpace/1', 'Scope/1');

save_system(modelName);
      

This realizes the same mathematical model: \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \), \( \mathbf{y}=\mathbf{C}\mathbf{x}+\mathbf{D}\mathbf{u} \), with \( m=2 \) input channels and \( p=2 \) output channels.

Mathematica supports state-space models natively and can compute time responses for MIMO systems.

A = { {0, 1}, {-2, -3} };
B = { {0, 0}, {4, -1} };
C = { {1, 0}, {0, 1} };
D = { {0, 0}, {0, 1/2} };

sys = StateSpaceModel[{A, B, C, D}];

tmax = 10;
u1[t_] := 1;
u2[t_] := Exp[-0.7 t];

u[t_] := {u1[t], u2[t]};

x0 = {0, 0};

(* State response and output response *)
xResp = StateResponse[sys, u[t], {t, 0, tmax}, x0];
yResp = OutputResponse[sys, u[t], {t, 0, tmax}, x0];

(* Check identity: y2(t) == x2(t) + 0.5 u2(t) *)
check = Simplify[
  yResp[[2]] - (xResp[[2]] + (1/2) u2[t])
];

check
      

The final check should simplify to 0 (or numerically evaluate near zero), confirming that the stacked output equation is being applied consistently.

Problem 1 (Build \( \mathbf{B} \) for multiple inputs): Consider the 3rd-order ODE \[ y^{(3)}(t) + 2\ddot{y}(t) + 5\dot{y}(t) + y(t) = 3u_1(t) + 4u_2(t) - 2u_3(t). \] Using phase-variable states \( x_1=y, x_2=\dot{y}, x_3=\ddot{y} \), construct \( \mathbf{A}\in\mathbb{R}^{3\times 3} \) and \( \mathbf{B}\in\mathbb{R}^{3\times 3} \).

Solution: By definition, \( \dot{x}_1=x_2 \), \( \dot{x}_2=x_3 \). The ODE gives \[ \dot{x}_3 = y^{(3)} = -y - 5\dot{y} - 2\ddot{y} + 3u_1 + 4u_2 - 2u_3 = -x_1 - 5x_2 - 2x_3 + 3u_1 + 4u_2 - 2u_3. \] Therefore

\[ \mathbf{A}= \begin{bmatrix} 0 & 1 & 0\\ 0 & 0 & 1\\ -1 & -5 & -2 \end{bmatrix}, \qquad \mathbf{B}= \begin{bmatrix} 0 & 0 & 0\\ 0 & 0 & 0\\ 3 & 4 & -2 \end{bmatrix}. \]

Problem 2 (Multiple outputs from one state): For the system in Problem 1, define two outputs: \( y_1(t)=y(t) \) and \( y_2(t)=\dot{y}(t) + 0.2u_2(t) \). Construct \( \mathbf{C}\in\mathbb{R}^{2\times 3} \) and \( \mathbf{D}\in\mathbb{R}^{2\times 3} \).

Solution: Since \( x_1=y \) and \( x_2=\dot{y} \), \[ y_1 = x_1,\qquad y_2 = x_2 + 0.2u_2. \] Hence

\[ \mathbf{C}= \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0 \end{bmatrix}, \qquad \mathbf{D}= \begin{bmatrix} 0 & 0 & 0\\ 0 & 0.2 & 0 \end{bmatrix}. \]

Problem 3 (Recover the ODE from a companion-form MIMO model): Let \( \mathbf{x}=[y,\dot{y}]^\top \) and \[ \dot{\mathbf{x}}= \begin{bmatrix} 0 & 1\\ -6 & -5 \end{bmatrix}\mathbf{x} + \begin{bmatrix} 0 & 0\\ 2 & -1 \end{bmatrix}\mathbf{u},\quad y = \begin{bmatrix}1 & 0\end{bmatrix}\mathbf{x}. \] Derive the 2nd-order ODE relating \( y \) to \( u_1,u_2 \).

Solution: From \( x_1=y \) and \( \dot{x}_1=x_2 \), we have \( x_2=\dot{y} \). The second state equation gives \[ \dot{x}_2 = \ddot{y} = -6x_1 - 5x_2 + 2u_1 - u_2 = -6y - 5\dot{y} + 2u_1 - u_2. \] Therefore the equivalent ODE is:

\[ \ddot{y}(t) + 5\dot{y}(t) + 6y(t) = 2u_1(t) - u_2(t). \]

Problem 4 (Dimension consistency check): Suppose \( n=4 \), \( m=3 \), \( p=2 \). If \( \mathbf{x}\in\mathbb{R}^4 \) and \( \mathbf{u}\in\mathbb{R}^3 \), state the required dimensions of \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \) and prove that each term in \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \), \( \mathbf{y}=\mathbf{C}\mathbf{x}+\mathbf{D}\mathbf{u} \) is well-defined.

Solution: The dynamics require \( \dot{\mathbf{x}}\in\mathbb{R}^4 \). Since \( \mathbf{x}\in\mathbb{R}^4 \), \( \mathbf{A}\mathbf{x} \) is defined only if \( \mathbf{A}\in\mathbb{R}^{4\times 4} \), producing a \( 4\times 1 \) vector. Since \( \mathbf{u}\in\mathbb{R}^3 \), \( \mathbf{B}\mathbf{u} \) is defined only if \( \mathbf{B}\in\mathbb{R}^{4\times 3} \), again producing a \( 4\times 1 \) vector. Their sum is therefore in \( \mathbb{R}^4 \). For outputs, \( \mathbf{y}\in\mathbb{R}^2 \). Thus \( \mathbf{C}\mathbf{x} \) requires \( \mathbf{C}\in\mathbb{R}^{2\times 4} \), producing a \( 2\times 1 \) vector, and \( \mathbf{D}\mathbf{u} \) requires \( \mathbf{D}\in\mathbb{R}^{2\times 3} \), also producing a \( 2\times 1 \) vector. Hence \( \mathbf{y}=\mathbf{C}\mathbf{x}+\mathbf{D}\mathbf{u} \) is well-defined. ∎

We extended the ODE-to-state-space conversion to the MIMO setting by: (i) treating inputs as vectors so \( \mathbf{B} \) acquires multiple columns, (ii) stacking outputs so \( \mathbf{C} \) acquires multiple rows, and allowing direct feedthrough via \( \mathbf{D} \), and (iii) proving equivalence to the original higher-order ODEs via direct differentiation and substitution. These constructions prepare us to work systematically with state-space realizations before introducing transfer functions and transfer matrices in Chapter 6.