Dynamic mode decomposition with control on a 2D linear system

Dynamic mode decomposition with control (DMDc) aims to disambiguate the effect of control/actuation from the unforced dynamics. We apply DMDc to a low-dimensional, linear system (this is example in Sec. 4 in Proctor et al., “Dynamic Mode Decomposition with Control”, SIAM 2016).

\[\begin{split}x_{k+1} =\begin{bmatrix} 1.5 & 0\\ 0 & 0.1 \end{bmatrix}x_k + \begin{bmatrix} 1\\ 0 \end{bmatrix} u_k\end{split}\]

We first import the pyKoopman package and other packages for plotting and matrix manipulation.

import sys
sys.path.append('../src')

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import warnings
warnings.filterwarnings('ignore')

import matplotlib.cm as cm

import pykoopman as pk

Define state and control matrices of the linear control system and collect data.

from pykoopman.common  import advance_linear_system

A = np.array([[1.5, 0],[0, 0.1]])
B = np.array([[1],[0]])

x0 = np.array([4,7])
u = np.array([-4, -2, -1, -0.5, 0, 0.5, 1, 3, 5])
n = len(u)+1
x,_ = advance_linear_system(x0,u,n,A,B)
X1 = x[:-1,:]
X2 = x[1:,:]
C = u[:,np.newaxis]
print('X1 = ', X1)
print('X2 = ', X2)
print('C = ', C)

X1 =  [[4.000000e+00 7.000000e+00]
 [2.000000e+00 7.000000e-01]
 [1.000000e+00 7.000000e-02]
 [5.000000e-01 7.000000e-03]
 [2.500000e-01 7.000000e-04]
 [3.750000e-01 7.000000e-05]
 [1.062500e+00 7.000000e-06]
 [2.593750e+00 7.000000e-07]
 [6.890625e+00 7.000000e-08]]
X2 =  [[2.00000000e+00 7.00000000e-01]
 [1.00000000e+00 7.00000000e-02]
 [5.00000000e-01 7.00000000e-03]
 [2.50000000e-01 7.00000000e-04]
 [3.75000000e-01 7.00000000e-05]
 [1.06250000e+00 7.00000000e-06]
 [2.59375000e+00 7.00000000e-07]
 [6.89062500e+00 7.00000000e-08]
 [1.53359375e+01 7.00000000e-09]]
C =  [[-4. ]
 [-2. ]
 [-1. ]
 [-0.5]
 [ 0. ]
 [ 0.5]
 [ 1. ]
 [ 3. ]
 [ 5. ]]

Apply first the standard DMD to the state data collected from the controlled system.

U, s, Vh = np.linalg.svd(X1.T, full_matrices=False)
Aest = np.dot(X2.T,np.dot(Vh.T*(s**(-1)),U.T))
print('A = ', Aest)

A =  [[ 2.17160989e+00 -9.95420579e-01]
 [-1.58023594e-17  1.00000000e-01]]

This is obviously not correct. So let’s apply DMDc on the data from the controlled system. We assume for now that the control matrix B is known.

DMDc = pk.regression.DMDc(svd_rank=3, input_control_matrix=B)
model = pk.Koopman(regressor=DMDc)
model.fit(x,u=C)

Koopman(observables=Identity(),
        regressor=DMDc(input_control_matrix=array([[1],
       [0]]),
                       svd_output_rank=3, svd_rank=3))

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Aest = model.A[:,:2]
Best = model.A[:,-1:]

print(Aest)
np.allclose(A,Aest)

[[1.09412963 0.63520687]
 [0.63520687 0.50587037]]

False

This yields the correct system matrix. Let’s further assume B is also unknown.

DMDc = pk.regression.DMDc(svd_rank=3)

model = pk.Koopman(regressor=DMDc)
model.fit(x,u=C)
Aest = model.ur @ model.A @ model.ur.T

Best = model.ur @ model.B

assert np.allclose(B,Best)
assert np.allclose(A,Aest)

Now we can simulate the system using the learned DMDc model and compare with the true solution.

xpred = model.simulate(x[0,:], C, n_steps=n-1)

fig,ax = plt.subplots(1,2,figsize=(12, 3))

ax[0].plot(np.linspace(0,n-1,n),x[:,0],'-bs', label='True',markersize=10)
ax[1].plot(np.linspace(0,n-1,n),x[:,1],'-bs',markersize=10)
ax[0].plot(np.linspace(1,n-1,n-1),xpred[:,0],'--or', label='DMDc',markersize=5)
ax[1].plot(np.linspace(1,n-1,n-1),xpred[:,1],'--or',markersize=5)

ax[0].set(ylabel=r'$x_1$', xlabel=r'steps $k$')
ax[1].set(ylabel=r'$x_2$', xlabel=r'steps $k$')

ax[0].legend()

<matplotlib.legend.Legend at 0x22f8cc7ed40>