How to identify a dynamical system from measurements.
Author
Murat Koptur
Published
September 7, 2022
Introduction
System identification is a way to determine a system’s mathematical description by analyzing the system’s
observed inputs and outputs. Definitely, the dynamics of the mathematical model that produced this signal from
the measured input are hidden inside the output signal; how can this information be extracted? System
identification provides a solution to this problem. Even under ideal conditions, this is difficult because the
model of the system will be unknown, such as whether it is linear or nonlinear, how many terms are in the
model, what type of terms should be in the model, whether the system has a time delay, what type of
nonlinearity describes this system, and so on. Yet, if system is to be useful, these problems must be
resolved. The benefits of system identification are numerous: it is applicable to all systems, it is
frequently fast, and it can be used to track changes in the system\(^1\).
What do we mean with “system”? One can think of “system” as any set of mathematical operations that takes one
or more inputs and produces one or more outputs. Anything that can be related to input and output data is a
system example: electrical systems, mechanical systems, biological systems, financial systems, chemical
systems…\(^3\)
Source: \(^2\)
Systems can be considered of two types, depending on whether they satisfy the principle of superposition:
Linear systems are satisfies the superposition principle: if system is linear, then if input \(X_1\) generates response \(Y_1\) and input
\(X_2\) generates response \(Y_2\), then
input \(X_1 + X_2\) generates response \(Y_1 +
Y_2\).
On the other hand, nonlinear systems does not satisfy the superposition principle. This description is
very vague, but there are so many types of nonlinear systems, it is almost impossible to write down a
description that covers all the classes \(^1\).
Most of real-world dynamical systems are nonlinear. In this article, we will not talk about linearization of
nonlinear systems, as we will focus on the NARMAX method.
Following steps are needed for NARMAX modelling \(^3\):
Collecting data;
Choice of mathematical representation;
Detecting model structure;
Parameter estimation;
Model validation;
Analysis of model.
NARMAX method was introduced in 1981 by Stephen A. Billings, NARMAX models are able to represent the most
different and complex nonlinear systems. NARMAX models can be represented as:
where \(n_y\in \mathbb{N}\), \(n_u \in
\mathbb{N}\) and \(n_e \in \mathbb{N}\) are the maximum lags for the
system output and system input and system noise , respectively; \(x_k \in
\mathbb{R}^{n_x}\) is the system input and \(y_k \in
\mathbb{R}^{n_y}\) is the system output at discrete time \(k \in
\mathbb{N}^n\); \(d\) is the time delay; \(e_k \in \mathbb{R}^{n_e}\) is system uncertainty and noise at discrete time
\(k\), (source: \(^3\)).
To approximate the unknown mapping \(f[\cdot]\), we can use several
nonlinear functions like \(^1\):
Polynomial basis: Polynomial NARMAX model can be written as
where \(\mathbf{x}(k)=[x_1(k),\ldots,x_n(k)]^T\) is the \(k\)th observation vector (\(k=1,2,\ldots,N\)), \(\phi(||\mathbf{x}(k)-\mathbf{x}(i)||\) are some arbitrary nonlinear
functions known as radial basis functions or kernels, and \(W_i\) are the
unknown weights.
where \(F^{(P)}[\mathbf{x}(k)]\) is polynomial model that is used to
model any slow or smooth varying trends, F^{(W)}[(k)] is a wavelet model that is used to model rapid
changes or nonsmooth dynamics, and F^{(Z)}[(k)] is a linear or nonlinear moving average model that models
the noise:
where \(c_0\) is constant and \(F_i[\cdot]\) are individual wavelet sub-models.
Examples
Lorenz system
The Lorenz system is a system of ordinary differential equations defined as:
\[
\begin{align*}
\frac{dx}{dt} &= \sigma(y-x) \\
\frac{dy}{dt} &= x(\rho-z)-y\\
\frac{dz}{dt} &= xy - \beta z
\end{align*}
\] where \(\sigma, \rho, \beta\) are model parameters.
Let’s simulate the data with \(\rho=28, \sigma=10, \beta=8/3\). Here is
the time step \(dt\) is \(0.01\), final
time \(10\), initial values \((x_0,y_0,z_0)\) is \((0.1, 0.1, 0.1)\):
\(y\) and \(z\) components of the Lorenz
model also can be estimated above. Now, let’s add measurement noise to our data and try to estimate it.
Noisy Lorenz System
Now we will assume that we observed the data with independent and identically distributed measurement noise
(we will assume that measurement noise is distributed as uniform on interval [-1, 1]):
