In earlier chapters, we modeled dynamics primarily via ordinary differential equations (ODEs) and transfer functions. The state-space viewpoint represents a dynamical system using a vector of internal variables (the \( \mathbf{x}(t) \)) such that, given the current state and future inputs, the future evolution is uniquely determined.

A general (possibly nonlinear, time-varying) continuous-time system is written as:

\[ \dot{\mathbf{x}}(t) = \mathbf{f}\!\big(\mathbf{x}(t), \mathbf{u}(t), t\big), \qquad \mathbf{y}(t) = \mathbf{g}\!\big(\mathbf{x}(t), \mathbf{u}(t), t\big), \]

where the \( \mathbf{u}(t) \) are inputs and \( \mathbf{y}(t) \) are outputs. Dimensions: \( \mathbf{x}(t)\in\mathbb{R}^{n} \), \( \mathbf{u}(t)\in\mathbb{R}^{m} \), \( \mathbf{y}(t)\in\mathbb{R}^{p} \).

flowchart TD
  A["Physical laws / ODEs"] --> B["Choose state variables x1..xn"]
  B --> C["Build state equation: x_dot = f(x,u,t)"]
  C --> D["Choose output equation: y = g(x,u,t)"]
  D --> E["Analyze / simulate / design later"]
        

Intuitively, a good choice of state variables often corresponds to energy storage in lumped-parameter systems (e.g., mass/inertia, springs, capacitors, inductors). This will be made precise in later chapters using energy-based modeling, but for now it serves as a practical guideline.

Consider a deterministic system with input \( \mathbf{u}(\cdot) \). Fix a time \( t_0 \). The system is said to have a state \( \mathbf{x}(t_0) \) if the following holds: knowing \( \mathbf{x}(t_0) \) and the future input \( \mathbf{u}(t), t\ge t_0 \) is sufficient to determine the future output \( \mathbf{y}(t), t\ge t_0 \).

A precise way to formalize this is via an equivalence relation on “pasts.” Let a past be the pair \( (\mathbf{u}_{(-\infty,t_0]}, \mathbf{y}_{(-\infty,t_0]}) \). Two pasts are equivalent at \( t_0 \) if, for every future input continuation, they generate identical future outputs.

Denote two pasts \( \mathcal{P}_1 \sim \mathcal{P}_2 \) when they are equivalent. A (minimal) state can be defined as an index of the equivalence class: \( \mathbf{x}(t_0) = \Xi(\mathcal{P}) \), where \( \Xi \) assigns a representative coordinate to each class.

flowchart TD
  P["Past (u,y) up to t0"] --> X["State x(t0) as compressed memory"]
  X --> FU["Future input u(t) for t>=t0"]
  FU --> YF["Future output y(t) for t>=t0"]
  X --> YF
        

For systems modeled by ODEs, this “sufficient memory” is captured by initial conditions. The next proposition states the key Markov-like property for state-space models.

Proposition 1 (State is sufficient): Suppose the system is described by \( \dot{\mathbf{x}}(t)=\mathbf{f}(\mathbf{x}(t),\mathbf{u}(t),t) \), \( \mathbf{y}(t)=\mathbf{g}(\mathbf{x}(t),\mathbf{u}(t),t) \). If \( \mathbf{f} \) is locally Lipschitz in \( \mathbf{x} \) (a standard condition guaranteeing uniqueness), then for any \( t_0 \) and any admissible input \( \mathbf{u}(\cdot) \), the pair \( (\mathbf{x}(t_0), \mathbf{u}(t), t\ge t_0) \) uniquely determines \( \mathbf{x}(t) \) and \( \mathbf{y}(t) \) for all future times where the solution exists.

Proof (sketch with standard ODE theorem):

Under local Lipschitz continuity in \( \mathbf{x} \), the initial value problem \( \dot{\mathbf{x}}(t)=\mathbf{f}(\mathbf{x}(t),\mathbf{u}(t),t),\ \mathbf{x}(t_0)=\mathbf{x}_0 \) has a unique solution on some interval about \( t_0 \) (Picard–Lindelöf theorem). Since \( \mathbf{y}(t)=\mathbf{g}(\mathbf{x}(t),\mathbf{u}(t),t) \), the output is uniquely determined as well. Therefore, no additional “history” beyond \( \mathbf{x}(t_0) \) is needed.

This proposition explains why state variables are often chosen as variables whose values at \( t_0 \) fully specify the system’s initial condition.

Many physical models initially appear as an \( n \)-th order ODE in a scalar variable \( y(t) \). State-space modeling systematically converts such models into an equivalent first-order vector ODE.

Consider the scalar ODE (input enters without derivatives, for simplicity):

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

Theorem 1 (Standard state construction): Define the state variables \( x_1(t)=y(t),\ x_2(t)=\dot{y}(t),\ \ldots,\ x_n(t)=y^{(n-1)}(t) \). Then the ODE is equivalent to the first-order system \( \dot{\mathbf{x}}(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}u(t) \), \( y(t)=\mathbf{C}\mathbf{x}(t) \), with:

\[ \mathbf{A} = \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},\quad \mathbf{B}= \begin{bmatrix} 0\\0\\ \vdots\\0\\ b \end{bmatrix},\quad \mathbf{C}= \begin{bmatrix} 1 & 0 & \cdots & 0 \end{bmatrix},\quad \mathbf{D}=0. \]

Proof:

By definition, \( \dot{x}_1=x_2 \), \( \dot{x}_2=x_3 \), …, \( \dot{x}_{n-1}=x_n \). Also, \( \dot{x}_n = y^{(n)} \). From the original ODE,

\[ \dot{x}_n(t) = -a_{n-1}x_n(t) - a_{n-2}x_{n-1}(t) - \cdots - a_1 x_2(t) - a_0 x_1(t) + b\,u(t). \]

Collecting these equations yields exactly \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}u \). Since \( x_1=y \), the output equation is \( y=\mathbf{C}\mathbf{x} \). This mapping is reversible because knowing \( \mathbf{x}(t) \) recovers \( y \) and its derivatives up to order \( n-1 \).

Interpretation: the state vector dimension equals the order of the ODE. In many engineering models, this dimension equals the number of independent initial conditions required to predict future trajectories.

The most important special case in control engineering is the linear time-invariant (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), \]

where \( \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} \).

The solution can be expressed using the state-transition matrix \( \boldsymbol{\Phi}(t,t_0) \). For LTI systems, \( \boldsymbol{\Phi}(t,t_0)=e^{\mathbf{A}(t-t_0)} \) (matrix exponential).

Theorem 2 (Variation of constants for LTI):

For any initial condition \( \mathbf{x}(t_0)=\mathbf{x}_0 \), the state satisfies:

\[ \mathbf{x}(t)=e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 + \int_{t_0}^{t} e^{\mathbf{A}(t-\sigma)}\,\mathbf{B}\mathbf{u}(\sigma)\, d\sigma. \]

Consequently,

\[ \mathbf{y}(t)=\mathbf{C}e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 + \mathbf{C}\int_{t_0}^{t} e^{\mathbf{A}(t-\sigma)}\,\mathbf{B}\mathbf{u}(\sigma)\, d\sigma + \mathbf{D}\mathbf{u}(t). \]

Proof:

Define \( \mathbf{z}(t)=e^{-\mathbf{A}t}\mathbf{x}(t) \). Differentiate:

\[ \dot{\mathbf{z}}(t) = -\mathbf{A}e^{-\mathbf{A}t}\mathbf{x}(t) + e^{-\mathbf{A}t}\dot{\mathbf{x}}(t) = -\mathbf{A}e^{-\mathbf{A}t}\mathbf{x}(t) + e^{-\mathbf{A}t}\big(\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t)\big) = e^{-\mathbf{A}t}\mathbf{B}\mathbf{u}(t). \]

Integrate from \( t_0 \) to \( t \):

\[ \mathbf{z}(t)-\mathbf{z}(t_0) = \int_{t_0}^{t} e^{-\mathbf{A}\sigma}\mathbf{B}\mathbf{u}(\sigma)\, d\sigma. \]

Multiply by \( e^{\mathbf{A}t} \) and use \( \mathbf{z}(t_0)=e^{-\mathbf{A}t_0}\mathbf{x}_0 \):

\[ \mathbf{x}(t)=e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 + \int_{t_0}^{t} e^{\mathbf{A}(t-\sigma)}\mathbf{B}\mathbf{u}(\sigma)\, d\sigma, \]

which proves the result. Substituting into \( \mathbf{y}(t)=\mathbf{C}\mathbf{x}(t)+\mathbf{D}\mathbf{u}(t) \) gives the output formula.

These solution formulas are foundational for simulation, stability analysis (via eigenvalues of \( \mathbf{A} \)), and control design (later chapters).

Consider the standard mass–spring–damper model (from earlier mechanical modeling):

\[ m\ddot{q}(t) + c\dot{q}(t) + k q(t) = F(t), \]

where \( q(t) \) is position, \( F(t) \) is applied force, and \( m,c,k \) are constants. Choose the state variables: \( x_1(t)=q(t) \), \( x_2(t)=\dot{q}(t) \), and input \( u(t)=F(t) \).

Then \( \dot{x}_1 = x_2 \) and \( \dot{x}_2 = \ddot{q} = -\frac{k}{m}q - \frac{c}{m}\dot{q} + \frac{1}{m}F \). Therefore:

\[ \dot{\mathbf{x}}(t)= \begin{bmatrix} 0 & 1\\ -\frac{k}{m} & -\frac{c}{m} \end{bmatrix}\mathbf{x}(t) + \begin{bmatrix} 0\\ \frac{1}{m} \end{bmatrix}u(t). \]

If the measured output is position \( y(t)=q(t)=x_1(t) \), then:

\[ y(t)= \begin{bmatrix} 1 & 0 \end{bmatrix}\mathbf{x}(t) + 0\cdot u(t). \]

This example illustrates the engineering meaning of state: position and velocity at time \( t_0 \) are exactly the information needed to predict future motion under a known force input.

This section implements the mass–spring–damper state-space model and simulates its response to a step input. The goal is to reinforce (i) the meaning of \( \mathbf{x}(t) \) as the evolving internal state and (ii) the practical use of the state-space form for simulation.

6.1 Python (NumPy/SciPy + Control Systems Library)

Common libraries for dynamical systems in Python include numpy, scipy.integrate, scipy.linalg, control (python-control), and sympy (symbolic derivations).

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

# Mass-spring-damper parameters
m, c, k = 1.0, 0.8, 4.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]])

def u_of_t(t):
    # Step force of magnitude 1
    return 1.0

def f(t, x):
    u = u_of_t(t)
    return (A @ x.reshape(-1,1) + B * u).ravel()

t0, tf = 0.0, 10.0
x0 = np.array([0.0, 0.0])  # initial position and velocity

sol = solve_ivp(f, (t0, tf), x0, max_step=0.01, dense_output=True)
t = np.linspace(t0, tf, 2000)
x = sol.sol(t).T
y = (x @ C.T).ravel()

plt.figure()
plt.plot(t, y)
plt.xlabel("t (s)")
plt.ylabel("y(t) = position")
plt.title("Mass-Spring-Damper Step Response (State-Space Simulation)")
plt.grid(True)
plt.show()

# Optional: state-space object via python-control
try:
    import control
    sys = control.ss(A, B, C, D)
    t2, y2 = control.step_response(sys, T=t)
    # y2 has shape (outputs, time) for MIMO; flatten for SISO
    plt.figure()
    plt.plot(t2, np.asarray(y2).ravel())
    plt.xlabel("t (s)")
    plt.ylabel("position")
    plt.title("Step Response via python-control ss()")
    plt.grid(True)
    plt.show()
except Exception as e:
    print("python-control not available or failed:", e)

6.2 C++ (Eigen + RK4 from scratch)

Common C++ tools for state-space simulation include Eigen (linear algebra) and Boost.Odeint (ODE integration). Below is a minimal RK4 integrator implemented from scratch using Eigen.

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

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

double u_of_t(double t) {
    (void)t;
    return 1.0; // unit step
}

int main() {
    // Parameters
    const double m = 1.0, c = 0.8, k = 4.0;

    MatrixXd A(2,2);
    A << 0.0, 1.0,
        -k/m, -c/m;

    VectorXd B(2);
    B << 0.0, 1.0/m;

    // Simulation settings
    double t0 = 0.0, tf = 10.0, h = 0.001;
    int N = static_cast<int>((tf - t0) / h);

    VectorXd x(2);
    x << 0.0, 0.0; // initial state: [position, velocity]

    auto f = [&](double t, const VectorXd& xvec) {
        double u = u_of_t(t);
        return A * xvec + B * u;
    };

    for (int i = 0; i <= N; ++i) {
        double t = t0 + i * h;

        // Output y = position = x1
        if (i % 1000 == 0) {
            std::cout << t << " " << x(0) << std::endl;
        }

        // RK4 step
        VectorXd k1 = f(t, x);
        VectorXd k2 = f(t + 0.5*h, x + 0.5*h*k1);
        VectorXd k3 = f(t + 0.5*h, x + 0.5*h*k2);
        VectorXd k4 = f(t + h, x + h*k3);
        x = x + (h/6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4);
    }

    return 0;
}

6.3 Java (Apache Commons Math + RK4)

In Java, typical choices include Apache Commons Math for linear algebra and ODE solvers, or EJML. Below is a compact RK4 simulation using arrays (from-scratch approach); you may replace array operations with library matrices if desired.

public class MassSpringDamperRK4 {

    static double uOfT(double t) {
        return 1.0; // unit step
    }

    // xdot = A x + B u
    static double[] f(double t, double[] x, double[][] A, double[] B) {
        double u = uOfT(t);
        double[] dx = new double[2];
        dx[0] = A[0][0]*x[0] + A[0][1]*x[1] + B[0]*u;
        dx[1] = A[1][0]*x[0] + A[1][1]*x[1] + B[1]*u;
        return dx;
    }

    static double[] add(double[] a, double[] b) {
        return new double[]{a[0] + b[0], a[1] + b[1]};
    }

    static double[] scale(double s, double[] a) {
        return new double[]{s*a[0], s*a[1]};
    }

    public static void main(String[] args) {
        double m = 1.0, c = 0.8, k = 4.0;

        double[][] A = new double[][]{
            {0.0, 1.0},
            {-k/m, -c/m}
        };

        double[] B = new double[]{0.0, 1.0/m};

        double t0 = 0.0, tf = 10.0, h = 0.001;
        int N = (int)((tf - t0)/h);

        double[] x = new double[]{0.0, 0.0}; // [position, velocity]

        for (int i = 0; i <= N; i++) {
            double t = t0 + i*h;

            if (i % 1000 == 0) {
                double y = x[0];
                System.out.println(t + " " + y);
            }

            double[] k1 = f(t, x, A, B);
            double[] k2 = f(t + 0.5*h, add(x, scale(0.5*h, k1)), A, B);
            double[] k3 = f(t + 0.5*h, add(x, scale(0.5*h, k2)), A, B);
            double[] k4 = f(t + h, add(x, scale(h, k3)), A, B);

            double[] incr = add(add(k1, scale(2.0, k2)), add(scale(2.0, k3), k4));
            x = add(x, scale(h/6.0, incr));
        }
    }
}

6.4 MATLAB / Simulink

MATLAB’s Control System Toolbox provides ss, step, and lsim. Simulink can represent the same model using integrator blocks (state accumulation) and gain/sum blocks (state equation).

% Mass-spring-damper state-space in MATLAB
m = 1.0; c = 0.8; k = 4.0;

A = [0 1; -k/m -c/m];
B = [0; 1/m];
C = [1 0];
D = 0;

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

t = linspace(0,10,2000);
u = ones(size(t));          % unit step force
x0 = [0; 0];

[y, t_out, x] = lsim(sys, u, t, x0);

figure; plot(t_out, y);
grid on; xlabel('t (s)'); ylabel('position');
title('Mass-Spring-Damper Response via lsim');

% Simulink (programmatic skeleton): creates a model with State-Space block
modelName = 'msd_statespace_model';
new_system(modelName); open_system(modelName);
add_block('simulink/Continuous/State-Space', [modelName '/StateSpace']);
set_param([modelName '/StateSpace'], 'A', mat2str(A), 'B', mat2str(B), 'C', mat2str(C), 'D', mat2str(D));

add_block('simulink/Sources/Step', [modelName '/Step']);
add_block('simulink/Sinks/Scope', [modelName '/Scope']);

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

set_param(modelName, 'StopTime', '10');
save_system(modelName);
% Run with: sim(modelName);

6.5 Wolfram Mathematica

Mathematica supports state-space models via StateSpaceModel and simulation via OutputResponse (for LTI) or NDSolve (general nonlinear ODEs).

(* Mass-spring-damper state-space in Mathematica *)
m = 1.0; c = 0.8; k = 4.0;

A = { {0, 1}, {-k/m, -c/m} };
B = { {0}, {1/m} };
C = { {1, 0} };
D = { {0} };

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

tmax = 10;
u[t_] := 1; (* unit step *)

(* Output response with zero initial state *)
y = OutputResponse[sys, u[t], {t, 0, tmax}];
Plot[y, {t, 0, tmax}, GridLines -> Automatic, PlotLabel -> "Step Response (position)"]

(* Direct ODE simulation (shows the meaning of state variables) *)
eqs = {
  x1'[t] == x2[t],
  x2'[t] == -(k/m) x1[t] - (c/m) x2[t] + (1/m) u[t],
  x1[0] == 0, x2[0] == 0
};

sol = NDSolve[eqs, {x1, x2}, {t, 0, tmax}];
Plot[Evaluate[x1[t] /. sol], {t, 0, tmax}, GridLines -> Automatic, PlotLabel -> "x1(t) = position"]

Problem 1 (State construction for an ODE): Consider the scalar second-order system \( \ddot{y}(t) + 3\dot{y}(t) + 2y(t) = u(t) \). (a) Choose a state vector \( \mathbf{x}(t) \). (b) Write \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}u \) and \( y=\mathbf{C}\mathbf{x} \).

Solution:

Choose \( x_1=y \), \( x_2=\dot{y} \). Then \( \dot{x}_1=x_2 \) and \( \dot{x}_2=\ddot{y}= -3\dot{y}-2y+u = -3x_2 - 2x_1 + u \).

\[ \dot{\mathbf{x}}= \begin{bmatrix} 0 & 1\\ -2 & -3 \end{bmatrix}\mathbf{x} + \begin{bmatrix} 0\\ 1 \end{bmatrix}u,\qquad y= \begin{bmatrix} 1 & 0 \end{bmatrix}\mathbf{x}. \]

Problem 2 (Sufficient memory interpretation): For the mass–spring–damper system \( m\ddot{q}+c\dot{q}+kq=F(t) \), explain why the pair \( (q(t_0),\dot{q}(t_0)) \) is sufficient to predict \( q(t) \) for \( t\ge t_0 \) when \( F(t) \) is known.

Solution:

The equation is second-order, so two initial conditions are required to specify a unique solution. Setting \( q(t_0) \) and \( \dot{q}(t_0) \) fixes those initial conditions; with a known forcing input \( F(t) \), the IVP has a unique solution (under standard smoothness assumptions), so no earlier history is needed.

Problem 3 (Matrix-exponential solution check): Consider the LTI system \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} \) with constant \( \mathbf{A} \). Show that \( \mathbf{x}(t)=e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 \) satisfies the differential equation and the initial condition.

Solution:

By definition of the matrix exponential, \( \frac{d}{dt}e^{\mathbf{A}(t-t_0)} = \mathbf{A}e^{\mathbf{A}(t-t_0)} \). Therefore:

\[ \dot{\mathbf{x}}(t)=\frac{d}{dt}\big(e^{\mathbf{A}(t-t_0)}\mathbf{x}_0\big) = \mathbf{A}e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 = \mathbf{A}\mathbf{x}(t). \]

At \( t=t_0 \), \( e^{\mathbf{A}(t_0-t_0)}=e^{\mathbf{0}}=\mathbf{I} \), hence \( \mathbf{x}(t_0)=\mathbf{x}_0 \).

Problem 4 (Output equation design): For the mass–spring–damper state \( \mathbf{x}=[q,\dot{q}]^\top \), provide matrices \( \mathbf{C},\mathbf{D} \) for each output: (a) output is position only, (b) output is velocity only, (c) output is both position and velocity.

Solution:

In all cases take \( \mathbf{D}=0 \). Then:

• (a) \( \mathbf{C}=\begin{bmatrix}1 & 0\end{bmatrix} \)
• (b) \( \mathbf{C}=\begin{bmatrix}0 & 1\end{bmatrix} \)
• (c) \( \mathbf{C}=\begin{bmatrix}1 & 0\\ 0 & 1\end{bmatrix} \)

Problem 5 (From third-order ODE to state-space): Let \( y^{(3)} + a_2 y'' + a_1 y' + a_0 y = b u \). Construct \( \mathbf{A},\mathbf{B},\mathbf{C} \) using the standard state choice.

Solution:

Set \( x_1=y \), \( x_2=\dot{y} \), \( x_3=y'' \). Then \( \dot{x}_1=x_2 \), \( \dot{x}_2=x_3 \), and \( \dot{x}_3=y^{(3)}=-a_2 x_3 - a_1 x_2 - a_0 x_1 + b u \).

\[ \mathbf{A}= \begin{bmatrix} 0 & 1 & 0\\ 0 & 0 & 1\\ -a_0 & -a_1 & -a_2 \end{bmatrix},\quad \mathbf{B}= \begin{bmatrix} 0\\ 0\\ b \end{bmatrix},\quad \mathbf{C}= \begin{bmatrix} 1 & 0 & 0 \end{bmatrix},\quad \mathbf{D}=0. \]

We defined the system state as the minimal internal memory required to predict future behavior, introduced state variables and state vectors, and showed how higher-order ODEs are systematically converted into first-order vector state equations. We presented the standard nonlinear and LTI state-space forms, derived the state-transition (matrix-exponential) solution for LTI systems, and implemented a canonical example across Python, C++, Java, MATLAB/Simulink, and Mathematica. These foundations enable the derivation of state-space models from physical laws in the next lesson.