Dynamic Mode Decomposition (DMD) was presented as a method to produce linear models from data generated through nonlinear dynamical processes by using Singular Value Decomposition (SVD) [1]. Later, this method was developed further to create Extended Dynamic Mode Decomposition (EDMD), which introduced the concept of using observable functions, nonlinear functions of state variables, as a method of augmenting the state space [2]. EDMD referenced the Koopman Operator as justification and a theoretical underpinning for lifting the state space. Decades prior, Bernard Koopman showed the existence of this operator that transforms nonlinear systems into linear systems [3]. Another extension of DMD has shown the viability of using DMD for control on non-autonomous systems [4]. This enabled complex nonlinear Model Predictive Control (MPC) to be converted to linear MPC [5], leading to numerous studies utilizing the methodology to real systems [6]- [9].
To improve the accuracy of the models based on the Koopman Operator, two avenues of research have formed. The first avenue regards the selection of the observable functions used for constructing a lifted state space, as these functions are a key ingredient in creating an accurate linear model. Various methods have been developed, including deep neural networks for learning effective observable functions [10]- [12] and optimization [13]. However, an efficient selection of observables does not solve all the issues that arise when attempting to construct an accurate linear model. The second avenue of research addressed the proper formulation of the linear transition matrix. It is known that unstable modes are involved in Koopman-based DMD models and their extensions although the underlying nonlinear systems are known to be stable [14]. Extensive studies have been done to create stable linear models to remedy the situations where an outright use of DMD would lead to the creation of an unstable linear model [14]- [17]. Recently, an extension of DMD, called Robust Dynamic Mode Decomposition (RDMD), utilizes statistical measures to suppress the effect of outliers on modeling the linear Koopman matrix [18].
A fundamental difficulty in constructing a proper linear model is data dependency. Least Squares Estimation (LSE), involved in all DMD based methods, often produces a significant bias, an over weighting of samples. This bias is due to the distribution, how prevalent similar data points are, of the dataset. To eliminate this dependency on distribution, the current work takes an alternative approach to LSE.
Recently, a new formulation of the Koopman Operator, termed Koopman Direct Encoding (DE), was produced [19]. This method directly encodes the nonlinear dynamics into the lifted linear model. Inner products of observable functions in composition with the nonlinear state transition function are used to construct the state transition matrix without use of LSE. While DE theoretically guarantees the exact linear model, it requires access to the nonlinear state equations, which are often not available in practical applications. The current work aims to fill the gap between DE and data-driven approaches.
There are four significant contributions presented in this work. The first is the conversion of the DE formula of the Koopman Operator to a data-driven formula. The second is a computational algorithm and proof of its convergence to the true inner products that constitute the DE formula. The third is numerical experiments that provide evidence that the proposed method, unlike EDMD, does not exhibit biases to data distribution, but can produce consistently higher accuracy compared to EDMD. Finally, the DDE algorithm is utilized in modeling a high order nonlinear system in combination with deep learning.
In this section, we give a brief overview of the Koopman Operator and dynamic mode decomposition [2], and introduce the direct encoding method for obtaining a Koopman Operator directly from nonlinear dynamics [19].
Consider a discrete-time dynamical system, given by
where x ∈ X ⊂ R n is the independent state variable vector representing the dynamic state of the system, f is a selfmap, nonlinear function f : X → X, and t is the current time step. Also consider a real-valued observable function of the state variables g : X → R. The Koopman Operator K is an infinite-dimensional linear operator acting on the observable function g :
where g • f is the composition of function g with function f :
A common data-driven method for constructing the operator is Extended Dynamic Mode Decomposition (EDMD) [2], where observables that are experimentally obtained or simulated from the governing equation of the system are augmented by including real-valued observable functions of the independent state vector x t . This collection of observables, z t is
. . .
where m is the order of the lifted state corresponding to the number of observable functions. Underpinned by the Koopman Operator theory, EDMD assumes the existence of a linear state transition matrix A relating z t+1 to z t , and determine A by solving a least squares regression that minimizes the Sum of Squared Error (SSE) through
Singular Value Decomposition (SVD) is used for the least squares optimization.
An alternative to the least squares estimate and EDMD is to obtain the exact A matrix by directly encoding the selfmap, nonlinear state transition function f (x) with an independent and complete set of observable functions through inner product computations. This Direct Encoding method is introduced next, while the full proof can be found in [19].
Let us first consider the case where g 1 , g 2 , g 3 , ... are orthonormal basis functions spanning a Hilbert space H. We assume that the self-map nonlinear function f (x) is continuous and that the composition of g j with f is also involved in the Hilbert space.
This implies that the function g j • f can be expanded in [g 1 , g 2 , g 3 , ...].
Concatenating g 1 , g 2 , g 3 , ... and
eq. ( 6) can be written in matrix and vector form.
where Ā is an infinite dimensional matrix consisting of the inner products involved in eq. ( 6),
Eq. ( 8) manifests that the state lifted to the infinite dimensional space zt makes linear state transition with matrix Ā. The observables g 1 , g 2 , g 3 , ... were assumed to be orthonormal basis functions in the above derivation. This assumption can be relaxed to an independent and complete set of basis functions spanning the Hilbert space. Hereafter, let [g 1 , g 2 , g 3 , ...] be an independent and complete set of basis functions spanning the Hilbert space.
It can be shown that the time evolution of lifted state z t is given by a constant matrix A f for the independent and complete set of basis functions [g 1 , g 2 , g 3 , ...].
The matrix A f can be computed directly from the selfmap, state transition function f (x) and an independent and complete set of observables [g 1 , g 2 , g 3 , ...] through inner product computations. Post-multiplying the transpose of z t to both sides of eq. ( 10) and integrating them over X yield:
which can be written as
where
Because the observables are independent, the matrix R is non-singular. Therefore, the matrix A f is given by
This formula for obtaining the matrix A f directly from the governing nonlinear state equation with the function f (x) and the independent observables through inner products, which are guaranteed to exist in Hilbert space H, is the Direct Encoding method.
The prevailing method for construction of the Koopman Operator, EDMD, is based on LSE and SVD. This method, however, cannot provide an unbiased estimate; the result is biased -dependent on the distribution of data within a dataset-as the Koopman Operator is being approximated [20]. This dependency on distribution of a dataset occurs because a core assumption of LSE is that the model structure is correct; when this assumption is violated, LSE is unable to create an unbiased estimator [21]. As the Koopman Operator is truncated in practical use, this assumption does not hold.
Non-uniform data distributions, that is datasets where distances between data points in the state space are not equidistant from their nearest neighbors, inevitably occur in practical applications. For a nonlinear dynamical system with a stable equilibrium, for example, data collected from experiments and/or simulation of the system tend to be dense in the vicinity of the equilibrium, as all trajectories that begin within a region of attraction converge to the equilibrium. Because LSE applies equal weighting to all data points, the model is heavily tuned to the behavior of the densely populated region.
The Direct Encoding method described previously enables us to obtain the exact linear state transition matrix A through inner product computations. As the formulation is based on integration over the entire state space, there is no bias towards particular parts of the domain.
However, the original form of the Direct Encoding method utilizes the nonlinear state equation, i.e. the self-map f (x), to compute the inner products. In practical applications, such a nonlinear function is not always available; only data are available. The objective of this section is to establish a computational algorithm to obtain the A matrix by numerically computing the inner products, ⟨g i , g j ⟩, ⟨g i • f, g j ⟩, from a given set of data.
The method presented consists of three operations.
• The integral of the inner products is reduced in range from the entire state space to the dynamic range encapsulated by the data. • The dynamic range is discretized with data points.
• The inner product integral is reduced to a weighted summation of the integrand evaluated at each data point multiplied by the volume ∆v associated to each point. Naturally, if data are densely populated in a small region, the discretized integral interval is small and thereby the volume also becomes small. Similarly, the volume tends to be larger where the data are sparse. In the summation, the integrand evaluated at individual data points are "weighted" by the size of the volume. This numerical inner product calculation prevents overemphasis of clustered data.
We present the data-driven encoding method (DDE) as an alternative data-driven method to DMD for calculating a finite order approximation of the Koopman Operator. The objective of this method is to compute the matrices R and Q in ( 13) and ( 14) from data. This entails the computation of inner products:
where
are assumed to be Riemann Integrable; the functions are bounded and continuous [22].
There are two data sets used for the inner product computation. The first data set is
Note that all the data values are finite, |x i | < ∞. As such, the integral interval of the inner products is finite in computing them from the data. To define the integral interval, we consider the dynamic range of the system, X D , determined from the data set D N . See Fig. 1. The dynamic range X D is defined to be the minimum domain in the space X that includes all the data points in D N , X D ⊃ D N , and that is convex. Namely, for any two states in
where 0 ≤ α ≤ 1. Each data point x i is mapped to f (x i ),
Fig. 1: Illustration of the dynamic range of a dataset defined by the convex hull containing all points in the set, partitioned using a triangulation method. The data points are in black and the convex hull that encapsulates all data points is in grey.
following the state transition law in eq.( 1). We assume that the transferred state, too, stays within the same dynamic range X D . Collecting all the transferred states yields the second data set.
This implies that the state space of the nonlinear system under consideration is closed within the dynamic range X D . With this dynamic range, we redefine our objective to compute the inner products over X D .
The integrals can be computed by partitioning the domain X D into many segments X 1 , • • • X P , as shown in Fig. 1.
We generate these segments by applying a meshing technique to the data set D N , where the n-dimensional coordinates of individual data points are treated as nodes of a mesh. Delaunay Triangulation, for example, generates a triangular mesh structure with desirable properties [23]. As illustrated in Fig. 1, each triangular element is convex and has no internal node. The volume of the dynamic range V (X D ) is the sum of the volumes of all the elements.
Accordingly, the integral R ij in eq.( 23) can be segmented to
Suppose that the p-th element has K p nodes, as shown in Fig. 2. Renumbering these nodes 1 through K p ,
The integrand G ij within the p-th element can be approximated to the mean of the K p nodes involved in the p-th element.
If Delaunay Triangulation is used, K p = n + 1. See Fig. 2. Substituting this into (23) yields the approximate value of R ij .
where Similarly, each component of the matrix Q can be computed by using the same meshing.
Note that F ij is evaluated by using the data points in both
where f (x[k p ]) ∈ D f N , thus not requiring evaluation of the nonlinear function f .
Consider the center of each partition, xp = Xp xdx/∆v p , and the distance between xp and each point,
See Fig. 2. The maximum distance from the center of the partition to each point that makes up the partition is
Consider a sequence of refining the approximate inner product integral Rij by increasing data points N . We can show that, as the number of partition P tends infinity and the maximum subintervals |∆x p | approach zero, the approximate inner product integral Rij converges to its true integral.
This formulation takes the form of weighted sums, specifically Riemann sums. Given functions that are bounded and continuous over the subdomain of interest, sequences of this form are known to have a common limit and thus converge upon refinement to the Riemann integral value over that subdomain, according to Numerical Integration theory [22, Section 1.5].
In the prior section, integrals (30) and (32) are presented as summations over partitions. This computation can be streamlined by converting the summations over partitions to the one over nodes. Consider node 3 associated to data point x 3 in Fig. 1, for example. This node is an apex of the 5 surrounding triangles. This implies that integrand G ij (x) is calculated or recalled 5 times in computing (30) and (32). This repetition can be eliminated by computing volume ∆v k associated to node k, rather than partition p : ∆v p . Namely, we compute
where I(k, p) is a membership function that takes value 1 when node k is an apex of triangle p, that is, node k is involved in partition p. Using this volume as a new weight we can rewrite (30) and (32) to be
Using this conversion, the computation can be streamlined and cleanly separated into three steps, as shown by pseudocode in Algorithm 1. The steps are: (1) Graph Creation: data are connected to create partitions of the domain using a mesh generator: lines 3 to 8, (2) Weighting Calculation: calculation of the weights for each data point: lines 10 to 17, and (3) Matrix Calculation: the calculation of the R and Q matrices to find the matrix A, lines 19 to 21.
In comparison to EDMD, this algorithm is notably slower. In terms of time complexity, EDMD uses Singular Value Decomposition (SVD) which is O(mn 2 ), while the current method of graph creation for DDE, Delaunay Triangulation, is O(m n/2 ), where n is the dimensions of the space and m is the number of points [24] [25]. In addition, as DDE is based on numerical integration, issues arise when the underlying nonlinear system becomes significantly high order.
IV. EXPERIMENTS In this section, the DDE algorithm is implemented for the sake of evaluating its validity and comparing its modeling accuracy to EDMD. Consider a 2nd order nonlinear system consisting of a pendulum with a nonlinear damper. See Fig. 3. The pendulum also bounces against walls with nonlinear compliance. The state variables for this system are x = [θ, θ] T , and the equation of motion can be written as:
where F k and F c are wall reaction moment and damping moment, respectively, where k = 200 and c = 1. We present a two part numerical experiment for this system. The first experiment regards variations in dataset size and distribution, and the second experiment varies the usage of observable functions.
The datasets tested are of three types: 1) Uniform: These datasets are composed of a rectangular dynamic range which is sampled uniformly, like an evenly divided grid. The range varies from θ = [-0. walls and the damping can vary from 0 to a significant value. See Fig. 3-(b), (c). 2) Gaussian: Data points are sampled with a finitesupport Gaussian distribution. The data are distributed non-uniformly with their highest density at the peak of the Gaussian placed at diverse locations. In addition, each dataset contains 100 data points uniformly distributed along the border of the dynamic range to guarantee the same dynamic range as the uniform datasets. Samples outside the dynamic range are excluded. 3) Trajectories: These datasets are composed of trajectories, beginning from 100 initial conditions that are simulated forward the same number of time steps. The dynamic range of this dataset differs from the two other dataset types.
The models constructed for DDE and EDMD use the same observable functions. The observable functions chosen are two dimensional radial basis functions (RBFs), uniformly distributed between the maximum and minimum values of each state variable in their respective dataset, and the state variables. The total order of the system is 27th order with 25 RBFs and 2 state variables.
A trajectory dataset graph is generated using Delaunay Triangles in DDE, shown in Fig. 4.
The accuracy of the models is tested through calculating sum of squared errors (SSE) for one-step ahead predictions over the dynamic range of the datasets. These error values are calculated for a uniform grid of points, similar to that used in the uniform datasets. A visualization of the SSE is plotted in Fig. 5. The results of these calculations are shown in Table I and II. In the computation, the dynamic range is discretized, and the SSE value of each point is summed.
For the Gaussian datasets, the test is run for eight iterations of each dataset to account for randomness and the average result is noted.
The second experiment varies the number of observable functions selected, thus increasing the order. In this experiment, the number of RBFs is varied through uniformly increasing the density of the centers of the function over the range of the dataset. The results are noted in Table III. In the table, the number detailing number of observables is the number of functions including the state variables.
From these results we can make the following observations.
• All the numerical experiments show that DDE outperforms EDMD in total SSE. • For uniform datasets, both models are nearly equivalent, though DDE has slightly lower SSE in all cases, shown in Table I. This result is expected as all data points are weighted equally in a uniform distribution. • For trajectory datasets, the total SSE converges for DDE with small dataset sizes. This result implies that the elements in the R and Q matrices of DDE, that is, the inner product integral computations, converge as the data size and the data density increase. This convergence is confirmed in Fig. 6, where several elements of the Q matrix are plotted against the data size. • The second experiment, regarding variations in observable function numbers demonstrates the effect of DDE remains even for significant increases in observable functions, referring to Table III. In all cases tested, DDE significantly outperforms EDMD over the dynamic range, as expected for a non-uniform dataset.
The use of deep neural networks for finding effective observable functions and constructing a Koopman linear model has been reported by several groups [10]- [12]. This method, sometimes referred to as Deep Koopman, is effective for approximating the Koopman Operator to a low-order model, compared to the use of locally activated functions, such as RBFs, which scale poorly for high-order nonlinear DDE can be applied to the neural network having an architecture similar to prior works [10]- [12]. See Fig. 7. The input layer receives training data of independent state variables. The hidden layers produce observable functions; these functions feed into the output layer consisting of linear activation units. This linear output layer corresponds to the A matrix that maps the observables of the current time to those of the next time step, i.e. the state transition in the lifted space. In the Deep Koopman approach, the output layer, that is, the A matrix, is trained together with observable functions in the hidden layers. This A matrix can be further improved by replacing it with the A matrix obtained from DDE, captioned in Fig. 7.
This Deep Koopman-DDE method is applied to a simulated cable manipulation system, similar to prior art [6]. The specific system used for the simulation experiment consists of one cable suspending a point mass, where a winch varies the length of one cable. The system is a 6-th order nonlinear, switched system where the cables go slack because of the unidirectional nature of cable tension. Because this system is higher order, it necessitates the application of Deep Koopman for the selection of observable functions to accurately represent the system linearly with a finite order. More details can be found in [26]. The network is constructed using PyTorch with the parameters shown in Table IV. From simulated trajectories starting at diverse initial conditions, 3,000 data points are drawn.
In this work, a new data-driven approach to generating a Koopman linear model based on the direct encoding of Koopman Operator (DDE) is presented as an alternative to dynamic mode decomposition (DMD) and other related methods using least squares estimate (LSE). The major contributions include: 1) The analytical formula of Direct Encoding is converted to a numerical formula for computing the inner product integrals from given data; 2) An efficient algorithm is developed and its convergence conditions to the true results are analyzed; 3) Numerical experiments demonstrate a) greater accuracy compared to EDMD, b) lower sensitivity to data distribution, and c) rapid convergence of inner product computation. Furthermore, the DDE method is incorporated to Deep Koopman, i.e. neural network based methods for construction of the Koopman Operator, for improving prediction accuracy. The current method, however, is for autonomous systems. The extension to systems with control is a challenge for the future and must be addressed rigorously, beyond utilizing an input as an observable.