Linear Systems

Estimation examples for linear systems.

Linear Systems via Dynamic Mode Decomposition

We will start by estimating the underlying dynamical system of a time discrete process based on some measurements via Dynamic Mode Decomposition. We will model a simple linear system of the form $u_{i+1} = A u_i$.

using DataDrivenDiffEq
using LinearAlgebra
using ModelingToolkit
using OrdinaryDiffEq

A = [0.9 -0.2; 0.0 0.2]
u0 = [10.0; -10.0]
tspan = (0.0, 11.0)

f(u,p,t) = A*u

sys = DiscreteProblem(f, u0, tspan)
sol = solve(sys, FunctionMap())

To estimate the underlying operator in the states $u_1, u_2$, we simply define a discrete DataDrivenProblem using the measurements and time, and then solve the estimation problem using the DMDSVD algorithm for approximating the operator.

prob = DiscreteDataDrivenProblem(sol)

res = solve(prob, DMDSVD(), digits = 1)
system = result(res)
Model ##Koopman#451 with 2 equations
States : x[1] x[2]
Parameters : p[1] p[2] p[3]
Independent variable: t
Equations
Differential(t)(x[1]) = p[1]*x[1] + p[2]*x[2]
Differential(t)(x[2]) = p[3]*x[2]

The DataDrivenSolution contains an explicit result which is a Koopman, defining all necessary information, e.g. the associated operator (which corresponds to our matrix A).

Matrix(system)
2×2 Matrix{Float64}:
  0.9          -0.2
 -1.07938e-17   0.2

In general, we can skip the expensive process of deriving a callable symbolic system and return just the basic definitions using the operator_only keyword.

res = solve(prob, DMDSVD(), digits = 1, operator_only = true)
(K = LinearAlgebra.Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}([0.1999999999999999, 0.9000000000000006], [0.05494422557947558 -0.9000000000000002; 0.19230478952816457 1.3877787807814457e-17]), C = [1.0 0.0; 0.0 1.0], B = Matrix{Float64}(undef, 0, 0), Q = [658.3410161960379 -135.16260071909275; -25.406503863502085 20.833333333333325], P = [658.3410161960379 -135.16260071909275; -25.406503863502085 20.833333333333325])

Where K is the associated operator given as its eigendecomposition, B is a possible mapping of inputs onto the states, C is the linear mapping from the lifted observables back onto the original states and Q and P are used for updating the operator.