System Identification for a Digital Twin: Koopman Operator and SINDy

Most digital twins are only as good as the equations inside them, and most plants do not come with trustworthy equations. The pump curve in the datasheet is for a new pump, the thermal model was fitted to a commissioning test three years ago, and the friction term in the robot joint was a guess. System identification for a digital twin is the discipline of learning the dynamics directly from the data the asset already produces, so that the twin’s model tracks the machine you have rather than the machine you bought.

This matters now because two families of methods have matured into usable open-source tooling: sparse regression (SINDy) that recovers a readable differential equation, and Koopman-based methods (dynamic mode decomposition and its extensions) that lift a nonlinear system into a space where it evolves linearly. They solve different problems, fail in different ways, and are often confused. You will leave this post knowing which one to pick for a given twin, how the pipeline is wired from historian to runtime, what the numbers look like on a worked example, and how to detect when the identified model has gone stale.

What this covers: the identification problem and incumbents, the SINDy and Koopman mechanisms, a reference architecture with drift monitoring, working code, a decision matrix, failure modes, and a rollout checklist.

Context and Background

Classical system identification, as codified in Ljung’s textbook and the MATLAB System Identification Toolbox, fits parametric models such as ARX, ARMAX, state-space and Box-Jenkins structures to input-output data. These are linear in the dynamics, well understood statistically, and good enough for a great deal of process control. They struggle when the plant is strongly nonlinear: a compressor near surge, a robot arm with Coulomb friction, a battery with a strongly state-dependent resistance, or a reactor with Arrhenius kinetics.

The usual response is physics-based modelling, which the site covers in Modelica vs Simulink for physics-based digital twin models. First-principles models extrapolate well and carry physical meaning, but they cost expert time and drift from reality as hardware ages. The opposite extreme is a black-box neural network, which fits quickly, extrapolates badly, and cannot be audited by a control engineer.

Data-driven dynamics sits between these. Three ideas dominate the research literature:

Sparse identification of nonlinear dynamics (SINDy), introduced by Brunton, Proctor and Kutz in a 2016 PNAS paper, assumes the governing equation is a sparse combination of candidate terms and finds it with sparse regression.
Dynamic mode decomposition (DMD), introduced by Schmid in 2010 for fluid flows, fits the best linear operator advancing snapshots of data one step in time.
The Koopman operator, defined by Bernard Koopman in 1931, is an infinite-dimensional linear operator that advances any measurement function of a nonlinear system. DMD is a finite-dimensional approximation of it, and extended DMD (EDMD, Williams, Kevrekidis and Rowley, 2015) and deep variants widen the approach.

The practical tooling is open source: PySINDy (described in the Journal of Open Source Software, 2020), PyDMD, and PyKoopman (a package for data-driven Koopman approximation, with an arXiv paper in 2023). We deliberately do not quote version numbers here because we could not confirm current releases this run; pin whatever your lockfile resolves and test it.

Identification also has a downstream consumer. A model that is only accurate for one step is useless for model predictive control in a digital twin, where the optimizer rolls the model forward tens of steps. And a model whose parameters shift with wear is exactly the input a remaining useful life estimator wants. Identification is the connective tissue between monitoring, prognostics and control.

What “identification” must deliver to a twin

A twin consumer needs four things from an identified model, and they are not equally served by every method:

Multi-step accuracy, not just one-step fit, because twins simulate forward.
Inputs, because a twin that cannot respond to actuator commands cannot answer what-if questions.
Interpretability or at least auditability, because operators and regulators ask why the twin predicted what it did.
A refresh path, because the asset will change.

The rest of this post is organised around those four requirements.

Reference Architecture: From Historian to Runtime Model

The short answer: a data-driven twin treats the identified model as a versioned, validated artifact. Raw sensor data is cleaned and resampled, an identification algorithm proposes a candidate model, a held-out validation gate accepts or rejects it, the accepted model runs inside a state estimator, and residual statistics feed a drift monitor that triggers refitting.

Figure 1: End-to-end identification pipeline. The validation gate and the drift loop are what turn a research notebook into a twin component.

The figure shows a closed loop rather than a one-off fit. The left half is offline or near-line work (cleaning, identification, validation). The right half is online (state estimation and drift detection). The arrow from drift monitor back to cleaning is the part most projects omit, and it is why many identified models quietly rot.

Data preparation decides the outcome

Identification methods differentiate or difference your data, and differentiation amplifies noise. A signal sampled at 100 Hz with measurement noise standard deviation of 0.5 percent of full scale gives a finite-difference derivative whose noise grows roughly as noise divided by the sampling interval. With a 10 ms interval, that is about 100 times the signal noise, in units per second. This is the single largest practical obstacle for SINDy, and it is why people use smoothing, total-variation regularized differentiation, or the weak (integral) formulation described later.

Sampling rate must be matched to the fastest dynamics you want to capture. A rule of thumb from the Nyquist argument is to sample at least ten times the highest dynamic frequency of interest for derivative-based identification, well above the bare factor of two. This is a heuristic, not a theorem. Industrial historians that compress with deadband or swinging-door algorithms are especially dangerous: they deliver irregular, lossy samples that break the uniform-time-step assumption of DMD. Pull raw high-rate data from the edge buffer for identification, not the compressed archive used for dashboards.

Excitation matters as much as rate. A plant held at steady state teaches the algorithm nothing about transients. Identification needs data in which the states and inputs move through the region the twin will operate in. In practice that means scheduled step tests, PRBS-style (pseudo-random binary sequence) perturbations, or harvesting naturally occurring transients such as start-ups and load changes.

Where the model runs

The identified model is normally small: a sparse polynomial with a dozen non-zero coefficients, or a linear matrix of a few hundred entries. It runs comfortably at the edge. The identification job, which may involve a library with thousands of columns or a singular value decomposition of a tall snapshot matrix, belongs in the cloud or on a workstation. The interface between the two is a model file (coefficients, library definition, scaling constants, validity envelope) and a version identifier, which pairs naturally with the artifact-registry practices used for any other deployed model.

The validity envelope

Every identified model should ship with an envelope: the ranges of states and inputs seen during training. Data-driven models are interpolators. Outside the envelope the polynomial terms of a SINDy model can explode and a Koopman model’s eigenvalues no longer describe reality. The twin runtime should check the current operating point against the envelope and flag extrapolation rather than silently returning a number.

How SINDy Works

The short answer: SINDy assumes that the time derivative of the state is a sparse linear combination of nonlinear candidate functions of the state. It builds a library of those functions, then uses sequentially thresholded least squares or a related sparse solver to find which few terms explain the data, returning a human-readable ordinary differential equation.

Figure 2: The SINDy loop. A library matrix is regressed against measured derivatives, a sparsity-promoting solver zeroes out most terms, and the identified ODE is integrated and compared with data.

The formulation

Collect state measurements x(t) at times t1 through tm, forming a matrix X with one row per time and one column per state variable. Estimate the derivative matrix dX/dt. Build a library matrix Theta(X) whose columns are candidate functions: a constant, each state, all pairwise products, squares, cubes, sines, cosines, or domain-specific terms such as 1/x or exp(-1/x) for kinetics. The model is

dX/dt = Theta(X) Xi

where Xi is a coefficient matrix with one column per state equation. Ordinary least squares would give a dense Xi that uses every term. The sparsity prior says most entries should be exactly zero, mirroring the fact that physical laws are written with few terms.

Sequentially thresholded least squares

The original algorithm, STLSQ, is almost embarrassingly simple. Solve least squares for Xi. Set every coefficient with absolute value below a threshold lambda to zero. Re-solve least squares using only the surviving columns. Repeat until the support stops changing. It converges in a handful of iterations and is fast because each solve is on a shrinking set of columns.

The threshold lambda is the key hyperparameter and it is scale-dependent, so columns should be normalized. Alternatives include LASSO (an L1 penalty), sparse relaxed regularized regression (SR3), and ensembling approaches that bootstrap the data and keep terms that appear in most fits. Ensembling is worth the cost in a twin context because it yields inclusion probabilities, which are an honest uncertainty signal.

A worked example with illustrative numbers

Take the Lorenz system, the standard SINDy benchmark: dx/dt = sigma(y – x), dy/dt = x(rho – z) – y, dz/dt = xy – beta z, with sigma = 10, rho = 28 and beta = 8/3. Expanding, the true x equation has two terms (-10x and 10y), the y equation has three (28x, -y, -xz), and the z equation has two (-8/3 z and xy). That is seven non-zero coefficients out of a third-degree polynomial library in three variables, which has 20 columns per equation, so 60 coefficients in all. SINDy’s job is to find those seven among 60.

In the original paper the authors recover the Lorenz equations from simulated data. The following figures are an illustrative estimate, not a benchmark we ran: with clean data at a small time step, a recovered coefficient typically lands within a fraction of a percent of the truth; with additive noise around the percent level on the states, plain finite-difference SINDy often keeps spurious terms or drops small true ones such as the -y term, because its coefficient of magnitude 1 is small next to 28 and 10. This is precisely the situation where the weak formulation and ensembling help. We make the point to calibrate expectations: SINDy is not magic against noise, and small true coefficients are the first casualties.

Beyond the basic algorithm

Several extensions matter for twins:

SINDy with control (SINDYc) adds the actuator signals u to the library, so the model becomes dx/dt = Theta(x, u) Xi. This is the variant a controllable twin needs.
Weak SINDy (Messenger and Bortz, 2021) multiplies the equation by test functions and integrates by parts, removing the need to differentiate noisy data. It is markedly more noise-robust than finite differences.
SINDy-PI handles rational functions, relevant for kinetics and flow, by identifying implicit equations.
Constrained SINDy enforces known physics such as energy conservation or symmetries as linear constraints on Xi, which both improves accuracy and keeps identified models physically plausible.

Choosing the library is the real modelling step

The sparse solver cannot find a term that is not in the library. Library design is therefore where domain knowledge enters. A thermal system gets temperature differences and fourth-power radiation terms. A rotating machine gets trigonometric functions of angle. A chemical reactor gets Arrhenius-type exponentials. Add too much and the library becomes ill-conditioned, because polynomial columns of degree three or four on correlated sensors are nearly collinear. Condition number growth is the practical limit on library size, more than compute.

How the Koopman Operator Approach Works

The short answer: the Koopman operator advances measurement functions of the state rather than the state itself. Although the underlying system is nonlinear, the operator is linear, so finding a good set of measurement functions turns a nonlinear problem into a linear one that can be analysed with eigenvalues and controlled with linear methods.

Figure 3: Koopman lifting. The nonlinear state is mapped to observables, evolves linearly in the lifted space, and is projected back to state space for prediction and control.

The idea in one paragraph

For a discrete-time system x(k+1) = F(x(k)), the Koopman operator K acts on a function g of the state by (Kg)(x) = g(F(x)). Whatever F does, K is linear in g, because applying it to a sum of functions gives the sum of the results. The price is dimension: the space of functions is infinite. Practical methods pick a finite dictionary of observables psi(x) and find a matrix K that advances their values: psi(x(k+1)) is approximately K psi(x(k)).

DMD as the base case

Dynamic mode decomposition is the simplest instance, where the observables are the raw states. Stack snapshots into X (columns x1 to x(m-1)) and X prime (columns x2 to xm), then solve X prime = A X in the least-squares sense, A = X prime times the pseudo-inverse of X. For high-dimensional data, such as a thermal field with thousands of sensor or CFD nodes, A is never formed directly. A truncated SVD of X reduces the problem to a small matrix, and its eigenvalues and eigenvectors are the DMD modes. Each eigenvalue lambda encodes a frequency and a growth or decay rate: the continuous-time rate is log(lambda) divided by the time step.

For a twin this is powerful for two reasons. First, the eigenvalues are physical diagnostics. A mode whose damping shrinks over months is a bearing or structural-health signal, and it connects directly to condition monitoring. Second, DMD scales to fields: a reduced-order model of a temperature or flow field can run in milliseconds where the CFD would run in hours.

DMD’s weakness is that if the state evolution is truly nonlinear, a single matrix A is a poor fit. It can capture a limit cycle’s fundamental frequency but not the amplitude dependence of a nonlinear oscillator.

EDMD and dictionaries

Extended DMD fixes this by choosing a richer dictionary: polynomials, radial basis functions, Hermite polynomials. The matrix K is then fitted by least squares between the lifted snapshots. The state is recovered either by including the state itself among the observables, so projection is a selection of coordinates, or by fitting a separate linear readout. EDMD converges, in the limit of infinite data and dictionary, to the Koopman operator’s projection. In finite practice, the dictionary quality dictates everything, and a poor dictionary produces eigenvalues outside the unit circle that make long-horizon prediction diverge.

Learning the observables

Deep Koopman methods (for example Lusch, Kutz and Brunton, 2018) train an autoencoder so that the encoder produces the observables and a small linear layer advances them, with loss terms for reconstruction, linear prediction and multi-step consistency. The benefit is that you do not hand-design the dictionary. The cost is everything you expect from neural networks: more data, hyperparameter sensitivity, no guarantee of eigenvalue stability, and a loss of the auditability that made Koopman attractive. We treat deep Koopman as a research tool for systems where hand dictionaries fail, not a default.

Koopman with control

For actuated systems, the lifted model becomes z(k+1) = A z(k) + B u(k), a linear system with a possibly nonlinear map back to the state. Korda and Mezić (2018) showed that this lifted linear predictor can be embedded in a linear model predictive controller, which solves a convex quadratic program instead of a nonlinear one. That is a major computational win for fast loops, and it is the bridge to the MPC post linked earlier. The caveat, discussed under failure modes, is that constraints on the original state become constraints on observables only if the state is among them.

HAVOK and time delays

When only part of the state is measured, which is the normal situation in industrial sensing, time-delay embedding recovers the missing information. The Hankel alternative view of Koopman (HAVOK, Brunton and colleagues, 2017) stacks delayed copies of a single measurement into a Hankel matrix and applies the SVD. The dominant singular vectors form coordinates in which chaotic or strongly nonlinear dynamics look like a linear system with an intermittent forcing term. For a twin with a single vibration or pressure signal this is often the only practical route.

Deeper Analysis: Running It, Comparing It, Keeping It Honest

This section turns the mechanisms into a working pipeline: code you can adapt, a method comparison, and the drift-monitoring loop that keeps an identified model trustworthy after deployment.

Runnable SINDy example

The snippet below identifies the Lorenz system with PySINDy. It follows the library’s documented SINDy interface (optimizer, feature_library, fit, print, simulate), but we did not execute it in this run and API details can shift between releases, so treat it as untested and check it against the version you install.

import numpy as np
from scipy.integrate import solve_ivp
import pysindy as ps

def lorenz(t, s, sigma=10.0, rho=28.0, beta=8/3):
    x, y, z = s
    return [sigma*(y - x), x*(rho - z) - y, x*y - beta*z]

dt = 0.002
t = np.arange(0, 20, dt)
sol = solve_ivp(lorenz, (t[0], t[-1]), [-8, 8, 27], t_eval=t, rtol=1e-10, atol=1e-10)
X = sol.y.T
X_noisy = X + np.random.default_rng(0).normal(0, 0.01 * X.std(0), X.shape)

model = ps.SINDy(
    optimizer=ps.STLSQ(threshold=0.1, alpha=0.05),
    feature_library=ps.PolynomialLibrary(degree=3),
    differentiation_method=ps.SmoothedFiniteDifference(),
)
model.fit(X_noisy, t=dt)
model.print()          # inspect which terms survived
print(model.score(X_noisy, t=dt))

# multi-step check: simulate from the first state and compare
X_sim = model.simulate(X_noisy[0], t[:2000])

Three habits separate a twin-grade script from a tutorial. Hold out a time block, not random rows, for validation because adjacent samples are correlated. Compare simulated trajectories, not just one-step derivative scores, because a model can fit derivatives well and still diverge when integrated. And for a chaotic system like Lorenz, compare invariants (attractor shape, Lyapunov-like statistics) rather than point-wise trajectories past a few Lyapunov times, since even a perfect model will diverge from a noisy initial condition.

Runnable DMD example

import numpy as np
from pydmd import DMD

# snapshots: rows are sensors, columns are time samples (assumed uniform dt)
snapshots = np.load("thermal_field.npy")      # shape (n_sensors, n_times)
dmd = DMD(svd_rank=12)                        # truncate to 12 modes
dmd.fit(snapshots)

eigs = dmd.eigs                               # discrete-time eigenvalues
dt = 1.0                                      # seconds between samples
growth = np.log(eigs).real / dt               # 1/s, negative is decay
freq = np.log(eigs).imag / (2*np.pi*dt)       # Hz
print(np.abs(eigs).max())                     # > 1 means unstable mode

The file name and shape are placeholders; this too is untested here. The diagnostic that matters is the largest eigenvalue magnitude. For a physically stable process that magnitude should be at or just below one. If it creeps above one, long-horizon predictions blow up, and the usual fixes are lower rank, a forward-backward DMD variant that reduces noise bias, or optimized DMD.

A worked numerical illustration for a twin

Suppose a motor-pump skid reports five signals at 50 Hz: shaft speed, inlet pressure, outlet pressure, flow and winding temperature. You collect 30 minutes of mixed operating data, which is 90,000 samples. A DMD fit with rank 8 reduces the 5-by-90,000 snapshot problem to an 8-by-8 matrix, a trivially small model. A cubic SINDy library over five states has 56 columns per equation (all monomials up to degree three in five variables, including the constant), so a full coefficient matrix holds 280 entries; if the sparse solution keeps something like 20 of them, the model is readable and cheap. These figures are arithmetic on library size, not measured results, and the retained count is an assumption for illustration.

The point of the arithmetic is the data-to-parameter ratio. With 90,000 samples and at most 280 coefficients, estimation variance is not the issue; excitation is. If the pump ran at a single speed for the whole 30 minutes, the speed-dependent terms are unidentifiable regardless of sample count, and the regression will assign them arbitrary values. Check the rank and condition number of Theta before trusting any coefficients.

SINDy versus DMD versus Koopman versus neural networks

The right method depends on what the twin must do, as the matrix below summarises. Ratings are our qualitative judgment from the published literature, not measured benchmark scores.

Criterion	SINDy	DMD	EDMD / Koopman	Neural ODE or LSTM
Output	Readable ODE	Linear matrix and modes	Lifted linear operator	Opaque network
Nonlinearity handled	Yes, if in library	No (single linear map)	Yes, via dictionary	Yes
Data needed	Low to moderate	Low	Moderate	High
Interpretability	Highest	High (modes)	Moderate	Lowest
Noise sensitivity	High unless weak form	Moderate (biased)	Moderate	Moderate
Control-ready	SINDYc, nonlinear MPC	Linear MPC	Linear MPC in lifted space	Nonlinear MPC
Extrapolation	Fair, risky	Poor	Poor to fair	Poor
Best for	Low-dimensional physics-like systems	Fields, oscillations, diagnostics	Fast-loop linear control of nonlinear plants	Messy plants with big data

A defensible default for a new asset: start with DMD for diagnostics and a baseline, try SINDy with a physics-informed library when the state dimension is below roughly ten, and move to EDMD or deep Koopman when you need a linear surrogate for MPC or the library approach cannot capture the behaviour. Use the neural network when you have abundant data and no one needs to read the equation.

Combining with physics: grey-box twins

The strongest practical pattern is not choosing between physics and data but composing them. Keep the known structure, such as conservation of mass and energy, as a skeleton, and use SINDy to identify only the unknown closure terms: the friction law, the heat transfer coefficient as a function of flow, the efficiency curve. This is sometimes called residual or hybrid modelling. The identified piece is small, so the data requirement and the extrapolation risk are both lower, and the result can be exported into a Modelica or FMI-based model as a custom component, consistent with the approach in the Modelica and Simulink post linked above.

Drift monitoring and refitting

An identified model is a snapshot of the machine on the day of the test. Fouling, wear, recalibrated sensors, new feedstock and seasonal ambient conditions all move the true dynamics. The twin therefore needs a runtime loop, shown below.

Figure 4: Runtime loop. The state estimator feeds residuals to a drift monitor, which raises an alarm that triggers refitting on a recent window and a validation report before the model is swapped.

The residual is the difference between the model’s one-step prediction and the next measurement. Under a healthy model, residuals should look like zero-mean noise with roughly the sensor variance. Test that with simple statistics: a CUSUM or EWMA control chart on the residual mean, a variance ratio test on the residual spread, and autocorrelation at a few lags. Autocorrelated residuals mean unmodelled dynamics; a shifted mean means a bias such as a calibration change.

Take care not to retrain on faults. If the drift is a developing failure, refitting absorbs it into the nominal model and the twin stops seeing the fault, which defeats condition monitoring. A robust design separates two outcomes: slow drift consistent with wear is absorbed by scheduled refits with bounded parameter change, while abrupt or large residual change raises a maintenance alert and freezes the model. Bounding how far any coefficient or eigenvalue may move per refit is a cheap, effective guard. Eigenvalue trajectories over time also make a useful health indicator in their own right, feeding the prognostics models discussed in the remaining-useful-life post.

Online and streaming variants

Batch refits are fine for slow assets. For faster adaptation, recursive formulations exist: online DMD updates the operator with each new snapshot using a rank-one update (analogous to recursive least squares with a forgetting factor), and streaming variants maintain an SVD incrementally. The forgetting factor trades tracking speed against noise: a short memory follows change quickly but chases noise. Choose it from the time scale of real drift on your asset, which is usually days to months, not seconds.

Uncertainty quantification

A point model with no error bars is a risky twin. Three practical routes: bootstrap or subsample the data and refit to get coefficient distributions (ensemble SINDy does this natively); use Bayesian formulations with sparsity-inducing priors such as the horseshoe or spike-and-slab; or calibrate empirically by propagating the model over many held-out windows and recording the growth of prediction error with horizon, then publishing that curve as the model’s trusted horizon. The last is the most operationally useful. A statement like “the model is within the sensor noise band for the first N steps, then degrades” lets a controller or operator decide how far to trust a forecast.

Trade-offs, Gotchas, and What Goes Wrong

Data-driven identification fails in predictable ways. Most can be caught by a disciplined validation gate; a few cannot be fixed by any algorithm.

Noise and differentiation. The most common failure is a SINDy model that fits derivative data well and produces garbage when integrated, because noise in the estimated derivatives drove it to keep spurious terms. Remedies: weak-form SINDy, smoothing before differentiation, ensembling, and always validating by integration rather than by derivative score.

Insufficient excitation. If the data never explore the operating region, the identified model is only valid on a thin slice. Closed-loop operation makes this worse, because a good controller deliberately suppresses variation, which leaves inputs correlated with states. Identifying a plant under feedback also biases estimates unless the controller’s influence is handled, either by injecting independent excitation or by using methods designed for closed-loop data.

Library and dictionary choice. SINDy cannot find terms outside its library, and EDMD with a poor dictionary produces spectral pollution, meaning spurious eigenvalues that do not correspond to real dynamics. These are modelling errors, not numerical ones, and no hyperparameter sweep will rescue them.

Spectral instability. A Koopman or DMD matrix with an eigenvalue outside the unit circle will diverge on long rollouts, even if the one-step error is tiny. Check the spectrum, and consider projecting eigenvalues onto the unit disk or fitting with a stability constraint when physics says the plant is stable.

Extrapolation. Polynomial SINDy models are dangerous outside the data range because high-degree terms dominate. Koopman models can be accurate only on the region where the dictionary was trained. The validity envelope from the architecture section is the mitigation.

Chaos and sensitivity. For chaotic or near-chaotic plants, long-horizon point prediction is impossible in principle. Measure success by statistical invariants and a bounded useful horizon, not trajectory match.

Non-stationarity and regime switching. A system with distinct modes, such as a machine running different recipes or a building in heating versus cooling season, is poorly served by a single model. Cluster the data by regime and fit one model per regime with a switching logic, or include the regime indicator in the library.

Sensor faults masquerading as dynamics. A stuck sensor or a calibration step shows up as a sudden change that the regression will happily absorb as a physical term. Run data-quality checks before identification and keep a human-reviewable log of what was excluded.

Governance. An auto-refitting model is a changing safety-relevant artifact. In regulated settings, treat each refit as a release: log the data window, the library, the coefficients, the validation metrics and the approver. The ability to diff two SINDy equations in plain text is a genuine advantage over neural alternatives here.

Anti-patterns

Fitting on the full dataset and reporting training error as accuracy. Tuning the SINDy threshold until the equation looks like the textbook one, which injects the answer. Treating a low one-step error as evidence of a good simulator. Retraining nightly without drift gating. Deploying with no envelope check.

Practical Recommendations

Begin with the question the twin must answer, because the method follows from it. If the answer is a diagnosis (is the machine changing?), start with DMD and track eigenvalues. If the answer is a simulation or what-if on a low-dimensional process with physics you partly know, use SINDy with a constrained, physics-informed library. If the answer feeds a fast controller on a nonlinear plant, evaluate EDMD with control and linear MPC, and keep a nonlinear fallback.

Invest first in data, not algorithms. A week spent on excitation design, sample-rate selection and cleaning beats a month tuning thresholds. Pull raw high-rate data, avoid lossy compression for identification windows, and record operating mode as a feature.

Define success as a multi-step, held-out, time-blocked criterion: for example, the simulated trajectory stays inside a stated multiple of sensor noise for a stated horizon on at least a stated fraction of validation windows. The numbers are yours to set from the application; the discipline of setting them before looking at results is what matters.

Rollout checklist

Confirm sampling rate against the fastest dynamics and keep raw data for identification.
Design or harvest excitation covering the full operating envelope; record the envelope.
Normalise features and inspect the condition number of the library before fitting.
Fit with an ensemble or weak form when noise is non-trivial; keep inclusion probabilities.
Validate by time-blocked simulation, not derivative score; check the eigenvalue spectrum.
Ship the model with version, envelope, library definition and trusted horizon.
Monitor residual mean, variance and autocorrelation; gate refits and bound parameter change.
Separate drift (absorb) from fault (alert and freeze).
Log every refit as an auditable release.
Keep a physics or neural fallback for operation outside the envelope.

For most industrial teams the first deployment should be modest: a diagnostic DMD or a SINDy grey-box closure term for one asset class, running in shadow mode next to the existing model. Promote it only after it has beaten the incumbent on held-out data through at least one seasonal or operating-mode change.

Frequently Asked Questions

What is system identification in a digital twin?

It is the process of building or updating the twin’s mathematical model of an asset from measured input and output data, rather than from design documents alone. The result might be a transfer function, a state-space model, a sparse nonlinear equation from SINDy, or a Koopman operator. Because machines wear and operating conditions change, identification is repeated, so the twin follows the real asset instead of its nameplate behaviour.

What is the difference between SINDy and the Koopman operator?

SINDy works in the original state space and returns a sparse nonlinear differential equation you can read and integrate. Koopman methods lift the state into a space of measurement functions where the evolution is linear, returning a matrix whose eigenvalues describe frequencies and decay. SINDy favours interpretability of the physics; Koopman favours linear tools for analysis and control. They can be combined, for example by using SINDy to discover Koopman observables.

How much data does SINDy need?

There is no universal number. The binding constraint is excitation: the data must exercise every term you want to identify. With a library of a few dozen columns, tens of thousands of samples are usually plenty if they span the operating region, while a million samples of steady-state operation can identify almost nothing. Noise level and sampling rate matter too. Validate with held-out time blocks and an ensemble to see whether coefficients are stable.

Can Koopman or SINDy models be used for model predictive control?

Yes. SINDy with control gives a nonlinear model usable in nonlinear MPC. Koopman with control gives a lifted linear model usable in linear MPC, turning an otherwise nonlinear problem into a convex quadratic program, as shown by Korda and Mezić in 2018. In both cases, validate multi-step prediction accuracy over the MPC horizon, enforce constraints carefully, and keep the operating point inside the model’s validity envelope.

How do I know when to retrain an identified model?

Monitor one-step prediction residuals. Healthy residuals are zero-mean with variance near the sensor noise and little autocorrelation. A control chart on the mean, a variance test and an autocorrelation check will flag drift. Retrain on a recent window when drift is slow and consistent with wear, with bounded parameter change. If the shift is abrupt or large, treat it as a possible fault, raise an alert and freeze the model.

Is a neural network better than SINDy for a digital twin?

Not categorically. Neural networks handle high-dimensional, messy data and need no hand-built library, but they require more data, extrapolate poorly, and are hard to audit. SINDy and DMD give compact, inspectable models that run on small edge hardware and are easy to diff between versions. Many teams use a physics or SINDy backbone with a neural residual term, getting interpretability where it matters and flexibility where it is needed.