Computationally Efficient Data-Driven Discovery and Linear Representation of Nonlinear Systems For Control

Koopman operator theory defines a way to transform any nonlinear dynamical system into an infinite-dimensional linear system [1]. Suppose we have an uncontrolled discrete-time nonlinear dynamical system defined as

where $x_{k}\in\mathcal{M}\subset\mathbb{R}^{n}$ is the system state, $k$ is the time index, and $\boldsymbol{f}$ is the function that evolves the states through state space. We then define observables (unrelated to observability from control theory), which here are real-valued functions of the system state: $g:\mathbb{R}^{n}\rightarrow\mathbb{R}$. The Koopman operator, $\mathcal{K}$, is defined such that, for any observable function $g$,

where $\circ$ is the composition operator. We can now apply this operator to the discrete-time system defined previously to arrive at:

Equation (3) shows that the Koopman operator propagates an observable function of any state, $g(x_{k})$, to the next time step.

It is not practical to apply the infinite-dimensional operator, so several methods have been developed to select a useful finite-dimensional Koopman observable space [13, 3, 4, 10]. Since it is non-trivial to select a set of observable functions that span a Koopman invariant subspace, we use EDMD [14] to find finite-dimensional approximations of the Koopman operator. EDMD involves lifting time-series data $x_{1},\dots,x_{M}$ into higher-dimensional space by selecting a set of observable functions and applying them to each $x_{k}$. Then, a linear matrix is fit to the higher-dimensional system. It is common to consider observable functions such as higher-order polynomials, radial basis functions, trigonometric functions, etc. [13, 15]. Note: in this work, we include the identity function $g(x)=x$ as an observable function (see section II C). We refer readers to Brunton et al. [16] for an in-depth study of Koopman operator theory.

In this work, deep neural networks (DNNs) are employed in conjunction with EDMD to approximate the Koopman operator in finite dimensions. First, rather than hand-selecting a set of observable functions, the DNN defines a set of observable functions. Then, the finite-dimensional approximation to the Koopman operator is calculated using least-squares regression (see Algorithm 1).

For a controlled system, the nonlinear dynamics are propagated with the following mathematical expression:

where $x_{k}\in\mathbb{R}^{n}$ and $u_{k}\in\mathbb{R}^{p}$. After lifting the states to higher dimensions with $\boldsymbol{\Phi}$, we wish to find a linear input-output system representation

where the matrix $\mathbf{K}$ approximates the Koopman operator. This is analogous to $\mathbf{A}$ and $\mathbf{B}$, the state transition and input matrices, respectively, in linear control systems. The goal of the deep learning framework (Fig. 1) is to learn $N$ observables/lifting functions and calculate the $\mathbf{K}$ and $\mathbf{B}$ matrices. We choose $N$ such that $N>>n$ and define

where $\mathbf{\phi}_{i}:\mathbb{R}^{n}\rightarrow\mathbb{R},i=1,...,N$ are the observable functions, which are defined by the DNN. Note that we concatenate the original states $\mathbf{x}_{k}$ with the observables from the DNN. This helps in two ways: (1) the original states can be easily extracted from the observables for control (see below), and (2) there is no need for a decoder to extract the original, thus saving valuable computational resources. Currently, there is no method to dictate the size of $N$ that would guarantee the best balance between simplicity and accurate approximation of the system; therefore, in most cases, $N$ is chosen empirically through trial and error. Some work has been done to study the size $N$ to find controllable systems [12, 4]. In this work, DNNs are applied with EDMD to find the Koopman operator.

To calculate the approximate Koopman operator $\mathbf{K}$ and the input matrix $\mathbf{B}$, the time history of measurement data for $M$ steps is arranged into snapshot matrices. The first matrix, $X$, is the state history from time $k=1$ to $k=M-1$, whilst the second matrix, $X^{\prime}$ is the same state history, right shifted by one-time step:

Mapping the measured states with observable functions leads to

Given the dataset, the matrices $\mathbf{K}$ and $\mathbf{B}$ can be found using least-sqaures

Applying the snapshot matrices of real data yields

therefore,

where $\dagger$ denotes the Moore-Penrose inverse of the matrix [17].

Because the Koopman operator calculated with least-squares is an approximation, as the observable functions do not span a Koopman invariant subspace, the predicted state is an approximation, which we denote with the $\hat{}$ symbol:

Now, we can extract the original states from the observables using a projection matrix $\mathbf{P}$ [5] yielding

where $\mathbf{I}_{n}$ is the $n$ x $n$ identity matrix and $\mathbf{0}_{n\text{x}N}$ is the $n\;\text{x}\;N$ zero matrix. As shown in [5] and [2], the observable functions not spanning a Koopman invariant subspace leads to error accumulation over time, which can lead to misleading predictions. However, if the prediction is corrected at each time step, then this error can be mitigated. This correction is applied by extracting the estimated state variable $\hat{x}_{k+1}$ at each time step with Equation 15 and then reapplying the observables to the extracted state variable with Equation 6.

One of the key improvements introduced in this paper is the use of a neural network with a recursive learning framework for learning the Koopman matrix to improve the training speed and reduce the time to predict the next step. Ultimately, more complex dynamical systems would inherently require more complex, larger Koopman operators in a higher dimensional lifted space, and the proposed architecture could facilitate that case.

As a baseline, we compare to an autoencoder, which is a type of NN architecture in which the output is the same as the input with an intermediate state in the middle that is “encoded.” Some of the most common uses of an autoencoder are in image compression, anomaly detection, and data denoising [18], where the encoder transforms the original data into the (generally compressed) encoded space, and the decoder reconstructs the original data from that encoded space. To use the autoencoder architecture with the Koopman operator, we encode the states into a higher dimensional lifted space with the encoder. One purpose of an autoencoder for Koopman operator theory is to recover the original states for control uses, which, in our case, is bypassed by concatenating the original states, as shown earlier.

Figure 1 shows the proposed architecture for the trained network with prediction and control in continuous time. Here, the given initial states pass through the DNN block to form the lifted states. The lifted states are then concatenated with original states $\mathbf{x}_{k}$ to form the new set of observables $\mathbf{\Phi}$. The states are propagated through the control block where $\mathbf{A}$ and $\mathbf{B}_{con}$ are the continuous form of Koopman matrix $\mathbf{K}$ and $\mathbf{B}$ (see section III B). The predicted states then pass through the correction block, resulting in the original states’ extraction. The entire loop runs until the states are regulated.

The complete algorithm can be seen in Algorithm 1, and the full PyTorch code is on Github¹¹1https://github.com/tiwari-research-group/Koopman-Control-No-Decoder.

The NN is trained using traditional batch training with a custom loss function. The loss is calculated recursively by applying the mean squared error (MSE) and EDMD. The prediction loss function is given as

Here $\hat{\boldsymbol{\Phi}}({x_{k+1}})$ is the predicted observable defined in Equation 14, where $\mathbf{K}$ and $\mathbf{B}$ are calculated recursively using the least-squares method as in Equation 13, and $\mathbf{\Phi}({x_{k+1}})$ are the observables calculated using the true data.

The data generation pipeline for the batch training can be seen in Algorithm 2. It can be noted that to collect data, the system is excited using random control input (+1 or -1). Also, random Gaussian noise is added to the training data to test the effectiveness of the proposed framework.

An LQR controller is employed with the learned dynamics. The objective of the LQR controller is to regulate the states of a system through the use of a quadratic loss function. The loss function $\mathcal{J}$ is given by

where $\mathbf{Q}=\begin{bmatrix}I_{n\times n}&.&.&0\\ .&.&.&0\\ .&.&.&0\\ 0&0&0&0\end{bmatrix}\in R^{N\times N}$ and $\mathbf{R}=1$ are the tuning matrices for states and control respectively. It can be noted that since the first $n$ rows of the observable $\mathbf{\Phi}(x)$ are the states themselves, it is straightforward to weight the controllable states. The LQR control effort is given as

This constant gain is found in MATLAB by calling the LQR function with the appropriate input variables. The trajectory of the LTI system can then easily be found by integrating it with the RK4 function. For purposes of control, we converted the Koopman matrix $\mathbf{K}$ and input matrix $\mathbf{B}$ in continuous form such that

where $\Delta t$ is the RK4 time step of $0.01$ seconds.

To demonstrate the applicability of our model, we have chosen to use the pendulum problem. It is a common example that is widely used due to its relatively simple yet nonlinear dynamics. The nonlinear state dynamics are given as follows:

To generate training data, the pendulum dynamics were integrated for 2 seconds using a Runge-Kutta integrator, using random initial conditions within $-2\leq\theta_{0}\leq 2$. For each initial condition, two data sets were created to provide the inputs and labels for the custom loss function. The first set, $X$ was the state history from time $k=0$ to $k=M-1$, whilst the second set, $X^{\prime}$, was the same state history, right shifted by one time step:

The training data is comprised of 8000 different random initial conditions. From this, 80% of the data is used in the training set, whilst 10% is used for the validation set and 10% is used as the test set. The duration for each initial condition is 2 seconds with a $\Delta t=0.01$.

The training data includes a random control effort of -1 or 1 to simulate the control effort that would be normally applied when deployed. This randomly excites the system and enables the extraction of the control input matrix $\mathbf{B}$. Additionally, to check the robustness of the model, random Gaussian noise is added to each state to simulate sensor measurement noise. This is achieved by adding noise with the same standard deviation to each state at every integration time step with $\mu=0.0158$ and $\sigma=0.00913$. See Algorithm 2.

In this subsection, we discuss the results of the proposed architecture. We also compare the performance of the proposed architecture with a traditional learning framework with an autoencoder. Figure 2 shows the prediction for only one step compared to the true dynamics. Figure 3 shows the prediction for 10 seconds, given only the initial states. It can be seen that learned dynamics closely follow true nonlinear dynamics. However, beyond 10 seconds, there is a significant divergence from the true dynamics, highlighting the limitations of the finite representation of the Koopman operator. The system response when subjected to LQR control is shown in Figure 4. It can be noted that other controllers such as a model predictive control can be employed in this case. Evident in Figure 5, the data that is used for the training of the DNN includes noise to mimic the inputs from a real-world sensor. As shown in the full prediction in Figure 3, the DNN learns to predict smoother dynamics than the noisy training data, shown in Figure 5. It can be seen in Figure 6 that the loss of both models observes the anticipated downward trend.

Figure 7 shows the 10 second comparison between the true dynamics, the proposed method, and the traditional autoencoder. It can be seen that the proposed model more closely follows the true dynamics. Figure 8 shows the absolute error between the true dynamics, the proposed method, and the traditional autoencoder solution. It can be noted that the proposed method is more accurate than the autoencoder framework. Additionally, the training time of our method was 30 seconds, and autoencoder training took 54 seconds. Therefore, the proposed method trains more efficiently than the traditional framework.