DMD with Control (DMDc)

Identifying forced linear systems x₁₊₁ = A xₜ + B uₜ from snapshot pairs

Why DMDc

Standard DMD assumes the system evolves autonomously. When the data comes from a system driven by an external input -- an actuated mechanical system, a circuit with an applied voltage, a neuron receiving synaptic current -- fitting plain DMD folds the forcing into the identified operator, producing a biased A that confuses the intrinsic dynamics with the input response. DMDc (Proctor, Brunton & Kutz 2016) separates the two by identifying the forced linear system:

where A (n × n) is the state-transition matrix and B (n × q) maps the q-dimensional control input into the state.

Mathematical Background

Snapshot pairs, not trajectories

Unlike core DMD, which takes one contiguous trajectory and splits it internally, dmdc takes explicit pair matrices: X₁ holds states at time t, X₂ the states one step later, and U the control input applied during each transition. Columns may therefore come from many concatenated trajectories, and individual pairs may be freely masked out by the caller.

Joint identification (unknown B)

Stack the state and input into a single regression matrix:

and solve the least-squares problem for both operators at once via the truncated SVD pseudo-inverse of Ω:

The first n columns of the solution are A, the remaining q are B.

Known-B identification

When the input coupling is known by construction, subtract the input response and solve only for A:

This uses fewer degrees of freedom and is immune to the closed-loop caveat below.

Optional output projection

For high-dimensional states, rank_output requests a second projection through the leading left singular vectors Û of X₂, yielding reduced operators in an r-dimensional basis:

DmdcConfig Options

Basic Usage

Simulate a forced two-state system, then recover A and B jointly from the snapshot pairs:

use koopman_dmd::{dmdc, DmdcConfig};

// Simulate x_{t+1} = A0 x_t + B0 u_t with
// A0 = [[0.9, 0.1], [0.0, 0.8]], B0 = [0.5, 1.0].
let m = 120;
let mut x1 = faer::Mat::<f64>::zeros(2, m);
let mut x2 = faer::Mat::<f64>::zeros(2, m);
let mut u = faer::Mat::<f64>::zeros(1, m);
let mut x = [1.0, -0.5];
for t in 0..m {
    // A persistently exciting, exogenous probe input
    let ut = (0.7 * t as f64).sin() + 0.5 * (2.3 * t as f64 + 1.0).cos();
    x1[(0, t)] = x[0];
    x1[(1, t)] = x[1];
    u[(0, t)] = ut;
    x = [0.9 * x[0] + 0.1 * x[1] + 0.5 * ut, 0.8 * x[1] + ut];
    x2[(0, t)] = x[0];
    x2[(1, t)] = x[1];
}

// Identify A and B jointly
let config = DmdcConfig { rank_input: Some(3), ..Default::default() };
let result = dmdc(&x1, &x2, &u, &config).unwrap();

// A and B are recovered to machine precision
println!("A = {:?}", result.a);   // [[0.9, 0.1], [0.0, 0.8]]
println!("B = {:?}", result.b);   // [[0.5], [1.0]]

All matrices in the result are real faer::Mat<f64>, ready for use in stepping loops without conversion.

Known B

When the input matrix is known by construction -- for example, a physical model tells you exactly how the input enters the state -- pin it and estimate only A:

let mut b_known = faer::Mat::<f64>::zeros(2, 1);
b_known[(0, 0)] = 0.5;
b_known[(1, 0)] = 1.0;

let config = DmdcConfig {
    rank_input: Some(2),
    known_b: Some(b_known),
    ..Default::default()
};
let result = dmdc(&x1, &x2, &u, &config).unwrap();
// result.b is a copy of known_b; only A was estimated

Autonomous Multi-Trajectory Fits

Passing a control matrix with zero rows (Mat::zeros(0, m)) turns dmdc into an autonomous identifier over explicit pairs -- something dmd() cannot do, since it requires one contiguous trajectory. This is the natural tool when your data is many short trajectories from different initial conditions:

// Concatenate pairs from any number of separate trajectories
let u = faer::Mat::<f64>::zeros(0, m);
let result = dmdc(&x1, &x2, &u, &DmdcConfig::default()).unwrap();
// result.b has zero columns; result.a is the multi-trajectory fit

DmdcResult Fields

Analyzing the Identified System

Because DmdcResult exposes raw eigenvalues rather than a full modal decomposition, 0.2.0 adds eigenvalue-slice variants of the analysis tools: stability_from_eigenvalues and spectrum_from_eigenvalues generalize dmd_stability / dmd_spectrum to any &[C64]:

use koopman_dmd::{stability_from_eigenvalues, spectrum_from_eigenvalues};

// Stability of the unforced dynamics
let stab = stability_from_eigenvalues(&result.eigenvalues, 1e-6);
println!("Spectral radius: {:.3}", stab.spectral_radius);

// Continuous-time frequencies and growth rates
let modes = spectrum_from_eigenvalues(&result.eigenvalues, None, result.dt);
for m in &modes {
    println!("freq={:.3} Hz, growth={:.3}", m.frequency, m.growth_rate);
}

Python and R

Both bindings expose the same identification modes and add a predict that steps the identified system xt+1 = A xt + B ut under a control input sequence, plus stability and spectrum analysis of the identified operator.

Python

import koopman_dmd

# Joint identification; known_b=B pins the input matrix instead
d = koopman_dmd.DMDc(x1, x2, u, rank_input=3)
d.a, d.b                      # identified matrices
d.eigenvalues                 # (rank_output, 2) re/im pairs
d.stability()                 # (is_stable, is_unstable, is_marginal, spectral_radius)
d.spectrum()                  # per-mode frequency / growth rate / stability

# Simulate under a new input sequence (columns set the horizon),
# or the zero-input free response with n_ahead
pred = d.predict(u=u_future)
free = d.predict(x0=x0, n_ahead=10)

# u=None fits an autonomous multi-trajectory model from the pairs
d_free = koopman_dmd.DMDc(x1, x2)

R

library(koopman.dmd)

# Joint identification; known_B = B pins the input matrix instead
fit <- dmdc(X1, X2, U, rank_input = 3)
fit$a; fit$b                  # identified matrices
dmdc_stability(fit)           # spectral radius and classification
dmdc_spectrum(fit)            # data frame of per-mode information

# Simulate under a new input sequence, or the free response
pred <- predict(fit, U = U_future)
free <- predict(fit, x0 = c(1, 1), n_ahead = 10)

# U = NULL fits an autonomous multi-trajectory model from the pairs
fit_free <- dmdc(X1, X2)

DMDc vs. the Other DMD Forms

• Core DMD -- one contiguous trajectory of an autonomous (unforced) system.
• DMDc -- a forced system with a measured input, or an autonomous system observed as many separate trajectories (q = 0 mode).
• Extended DMD -- nonlinear autonomous systems, via dictionary lifting; composable with core DMD's pipeline.
• Hankel-DMD -- scalar or low-dimensional measurements of a higher-dimensional system, via delay embedding.