Closed-Loop Koopman Operator Approximation

This paper proposes a method to identify a Koopman model of a feedback-controlled system given a known controller. The Koopman operator allows a nonlinear system to be rewritten as an infinite-dimensional linear system by viewing it in terms of an infinite set of lifting functions. A finite-dimensional approximation of the Koopman operator can be identified from data by choosing a finite subset of lifting functions and solving a regression problem in the lifted space. Existing methods are designed to identify open-loop systems. However, it is impractical or impossible to run experiments on some systems, such as unstable systems, in an open-loop fashion. The proposed method leverages the linearity of the Koopman operator, along with knowledge of the controller and the structure of the closed-loop system, to simultaneously identify the closed-loop and plant systems. The advantages of the proposed closed-loop Koopman operator approximation method are demonstrated in simulation using a Duffing oscillator and experimentally using a rotary inverted pendulum system. An open-source software implementation of the proposed method is publicly available, along with the experimental dataset generated for this paper.


Introduction
Using Koopman operator theory [1][2][3][4], a finite-dimensional nonlinear system can be rewritten as an infinite-dimensional linear system.The Koopman operator itself serves to advance the infinitedimensional system's state in time, much like the dynamics matrix of a linear system.Since, in practice, working with an infinite-dimensional operator is intractable, the Koopman operator is approximated as a finite-dimensional matrix, often using data-driven methods.While the Koopman operator was originally proposed almost one hundred years ago, recent theoretical results [2][3][4] paired with modern computational resources have led to it becoming a popular tool for identifying dynamical systems from data.
In a typical Koopman workflow, a set of lifting functions is first chosen.In some situations, these lifting functions are inspired by the dynamics of the system [5,6], while in others, they are taken from a standard set of basis functions, like polynomials, sinusoids, or radial basis functions [7,8].Lifting functions can also be chosen to approximate a given kernel [9][10][11].Time delays are often incorporated into the lifted state as well [7,12].Once the chosen lifting functions are applied to the data, linear regression is used to approximate the Koopman matrix, often with some form of regularization [7].
One of the key properties of the Koopman operator is its linearity, making it possible to leverage existing linear systems theory when identifying or analyzing Koopman systems.Moreover, the linear Koopman representation of the nonlinear system allows linear control design methods to be leveraged [5][6][7][12][13][14].Koopman operator approximation methods are also closely related to classical system identification methods [12], most notably subspace identification methods [15, §3].This similarity means that well-established system identification techniques can often be adapted to enrich the Koopman workflow.One area of particular interest is the identification of closed-loop systems [16][17][18].It is often impossible to run an experiment on a plant without a feedback control system, either because the plant is inherently unstable, or because it is unsafe or impractical to do so [18, §1.1, §17.3].In some cases the control loop is simply ignored and control outputs are treated as exogenous input signals to the plant [18, §13.5].However, this simplification, known as the direct approach, is not without its issues.The feedback loop introduces correlations between the plant's output and its input, leading to biased estimates of the plant's dynamics [19, §11.1].In cases where this bias is small, ignoring the control loop may be acceptable.In cases where the bias is not small, closed-loop system identification methods must be employed.These methods identify the closed-loop system as a whole, where the reference signal treated as the exogenous input, and then extract the plant system for use on its own [19, §11.1].Methods that use knowledge of the controller are known as indirect approaches [18, §13.5].In contrast, methods that identify an unknown controller along with the plant are categorized as joint input-output approaches [18, §13.5].While closed-loop system identification can lead to more complex plant models, in many situations it is the appropriate tool to use to avoid biased models.In some situations, closed-loop system identification must be used due to the inability to operate the plant in an open-loop fashion.
The key contribution of this paper is the adaptation of closed-loop system identification methods to Koopman operator approximation.The proposed method, summarized in Figure 1, is a variation of the indirect approach, wherein a Koopman model of the closed-loop system is first identified, and then the Koopman system modelling the plant is extracted given knowledge of the controller.The significance of this work is that a system that cannot be operated in open-loop without a feedback controller can now be identified using Koopman operator theory.The effectiveness of this method is demonstrated on a motor drive suffering from nonlinear gearbox vibrations, which would be difficult to identify in an open-loop setting.
The remainder of this paper is as follows.Section 2 outlines the necessary background information on the Koopman operator and its associated approximation techniques.Section 3 derives the proposed closed-loop Koopman operator approximation method.Section 4 describes the experimental setup used to generate the motor drive training dataset used in this paper, including the control sys- CL Koopman matrix: Plant Koopman matrix: Plant Koopman system: tem used on the drive.Section 5 outlines the important design decisions made in the Koopman operator approximation procedure, namely the lifting function selection rationale and the hyperparameter optimization method.Section 6 analyzes the experimental results of the closed-loop Koopman operator approximation method.The paper's conclusions are presented in Section 7.

Koopman operator theory
Consider the nonlinear difference equation where x k ∈ M evolves on a manifold M ⊆ R m×1 .Also consider an infinite number of scalarvalued lifting functions, ψ : M → R, which span an infinite-dimensional Hilbert space H.The Koopman operator, U : H → H, is a linear operator that composes all lifting functions ψ ∈ H with f(•), thereby advancing them in time by one step.That is [20, §3.2], Evaluating (2) at x k reveals that the dynamics of (1) can be rewritten linearly in terms of ψ as A finite-dimensional approximation of (3) is where ψ : M → R p×1 is the vector-valued lifting function, U ∈ R p×p is the Koopman matrix, and k is the residual error.

Koopman operator theory with inputs
The definition of the Koopman operator can be modified to accommodate nonlinear difference equations with exogenous inputs.Consider the difference equation where the state is x k ∈ M ⊆ R m×1 and the input is u k ∈ N ⊆ R n×1 .The lifting functions are now ψ : M × N → R and the Koopman operator U : H → H now satisfies where = u k if the input has state-dependent dynamics, or = 0 if the input has no dynamics [20, §6.5].Let the vector-valued lifting function ψ : M × N → R p×1 be partitioned as where the state-dependent lifting functions are ϑ : M → R p ϑ ×1 , the input-dependent lifting functions are υ : M × N → R p υ ×1 , and the lifting function dimensions satisfy p ϑ + p υ = p.With an exogenous input, substituting the state and input into (6) results in [20, §6.5.1] where U = A B .Expanding (8) yields the familiar linear state-space form, When designing lifting functions, the first m lifted states are often chosen to be the state of the original difference equation, . . .
This choice of lifting functions makes it easier to recover the original state from the lifted state.
Another common design decision is to leave the input unlifted when identifying a Koopman model for control.That is [12], However, recent work has also shown that bilinear input-dependent lifting functions are a better alternative for control affine systems [21].
An output equation, where , can also be considered to incorporate the Koopman operator into a true linear system with input, output, and state.While D = 0 in all cases, C can be chosen to recover the original states, or any other desired output.If the original state is not directly included in the lifted state, C can instead be determined using least squares [12, §3.2.1].The Koopman system is therefore where min ∼ denotes a minimal state space realization [22, §3.2.1].

Approximating the Koopman operator from data
To approximate the Koopman matrix from a dataset D = {x k , u k } q k=0 , consider the lifted snapshot matrices where ).The Koopman matrix that minimizes where (•) † denotes the Moore-Penrose pseudoinverse.
Often, Tikhonov regularization is included to improve the numerical conditioning of U by penalizing its squared Frobenius norm [23].The cost function then becomes where β is the regularization coefficient.

Closed-loop approximation of the Koopman operator
In this section the proposed closed-loop Koopman operator approximation method is outlined in detail.First, the Koopman representation of the closed-loop system is derived in terms of the known state-space matrices of the controller and the unknown Koopman matrix of the plant.Then, a method for extracting the unknown plant Koopman matrix from the identified closed-loop system is described.This method is inspired by the joint-input output subspace identification method of [19, §11.2.2].However, since the state-space matrices of the controller are required to recover the plant model, it is closer to an indirect system identification approach.

Formulating the closed-loop Koopman system
Consider the Koopman system modelling the plant, along with the linear system modelling the known controller, In the Koopman system, only A p and B p are determined from experimental data.The remaining state-space matrices are chosen to be such that C p recovers the original states of the nonlinear system being modelled, which are the first m states in which yields the series interconnection of the controller and the plant with a feedforward signal f k .This feedforward signal is entirely exogenous, and could be generated by a nonlinear inverse model of the plant.The new Koopman system's input includes the controller input and feedforward input, resulting in while its output is simply the plant output, ζ k = ζ p k .The new system's state includes the controller state and plant state, resulting in The cascaded plant and controller systems are depicted in Figure 2. The state space representation of the plant becomes The state space representation of the series-interconnected system is therefore or equivalently, Next, consider the negative feedback interconnection, where r k is an exogenous reference signal.The input of this new feedback-interconnected system is defined to be as depicted in Figure 3. Substituting (36) into (26) yields Substituting (39) into (35) results in a new output equation, where the feedback interconnection is well-posed if and only if . Substituting (41) into (39) yields Substituting (43) back into (34) results in a new state equation, Rearranging yields the state-space representation of the closed-loop system, with the state equation being and the output equation being Substituting in D p = 0 from (24) implies that Q = 1 and yields the simplified state-space repre-sentation, which is always a well-posed feedback interconnection, since

Extracting the plant Koopman system
Consider the closed-loop lifted dataset D = {ϑ k , υ k } q k=0 , where the matrix approximation of the Koopman operator is Since C p is known, A p can be recovered using U f 22 and U f 23 with Similarly, B p , U f 21 , U f 23 , and U f 24 are related via the expression As a result, B p can be recovered with where (•) † denotes the Moore-Penrose pseudoinverse.

Training setup
The advantages of the proposed closed-loop Koopman operator approximation method are demonstrated using an experimental dataset gathered from a motor drive with Harmonic Drive gearbox.The motor drive, pictured in Figure 4, exhibits nonlinear vibrations due to the dynamics of the gearbox.Since it is not practical or safe to operate such a motor drive in open-loop, it is operated using a simple linear position and velocity controller.This section describes how the dataset was collected and details the control system used to run the experiments.The sections that follow show that a linear model of the closed-loop system is unable to identify the nonlinear gearbox vibrations, while a Koopman model of the closed-loop system is able to model them.

Dataset
The complete motor drive dataset consists of 18 training episodes and 2 test episodes.The reference position, reference velocity, measured position, measured velocity, and measured torque are recorded at 1 kHz.The positions and velocities are recorded in rad and rad/s respectively at the output shaft of the gearbox, while the measured torque is recorded as a fraction of the drive's full-scale torque.
Figure 5 shows a portion of one of the test inputs used to characterize the closed-loop system.The drive's control software accepts position checkpoints, and generates smoothed trapezoidal velocity trajectories to reach them.Each episode contains 10 pseudorandom position checkpoints within one revolution of the gearbox output shaft in either direction.The resulting velocity profile resembles a smoothed pseudorandom binary sequence.For all episodes, the maximum allowed velocity and acceleration settings are selected in the drive's control software.

Controller
The drive's motor controller consists of an inner proportional-integral velocity loop paired with an outer proportional position loop.While the controller gains are known, its state-space matrices must be derived for later use in the closed-loop Koopman operator approximation procedure.
The controller can be represented by the multi-input single-output differential equation,  where θe (s) = θr (s) − θ(s) is the velocity error, θ e (s) = θ r (s) − θ(s) is the position error, θr (s) is the velocity reference, θ r (s) is the position reference, θ(s) is the motor velocity, θ(s) is the motor position, and τ (s) is the motor torque.The velocity controller gains are K θ p and K θ i , while the position controller gain is K θ p .Discretizing (55) using the forward Euler approximation, where and ∆t is the timestep, yields In the time domain, (58) is To find a state-space representation of (60), let the controller state be The state-space equations are then where T and y c k = τ k .Using the state-space matrices A c , B c , C c , and D c , the plant's Koopman representation can be extracted from that of the closed-loop system.

Training design choices
Two key design choices must be made to identify an accurate Koopman model of the motor drive system.First, a set of lifting functions must be selected.In this case, a simple set of sinusoidal lifting functions based on the dynamics of the Harmonic Drive gearbox are chosen.Second, a hyperparameter optimization scheme must be developed minimize the model's trajectory prediction errors for previously unseen datasets.

Choice of lifting functions
Harmonic Drive gearboxes are known to generate vibrations at a frequency once and twice the input frequency [25,26].Let the outputs of the motor drive system be the gearbox output angle θ(t) in rad and velocity θ(t) rad/s, and let the input be the desired torque τ (t) as a fraction of the full-scale torque.The nonlinear vibration can be modelled as an exogenous load torque disturbance, where α 1 and α 2 are vibration amplitudes, φ 1 and φ 2 are phase shifts, and r is the reduction ratio of the gearbox.In this case, the gearbox multiplies torque by a factor of 100, so the reduction ratio is r = 100.
These vibrations manifest themselves as tracking errors in the position and velocity of the drive.They also appear in the torque command due to the feedback action of the controller.Figure 6 shows the power spectral density of the drive's velocity tracking error during a constant 0.5 rev/s movement, which is the drive's maximum speed.Since the gearbox reduction ratio is known to be r = 100, the resulting gearbox input velocity is 50 rev/s.Velocity tracking errors therefore occur at integer multiples of 50 Hz, with decreasing power as frequency increases.By several orders of magnitude, the most significant tracking errors occur at 50 Hz.Inspired by (64), the plant's lifting functions are chosen to be Since the 50 Hz error in Figure 6 is by far the largest, only that sinusoidal term is included in the lifting functions.The vibration amplitude α is not required, as it becomes a part of the Koopman matrix during the regression process.Additionally, it was found experimentally that including higher-frequency harmonics has a minimal effect on the system's overall prediction error while complicating the hyperparameter optimization procedure.
For the sake of comparison, a set of linear lifting functions is also considered, that being For both choices of ϑ p k , the lifting functions of the closed-loop system are where the controller state x c k is computed from τ k , θ k , and θ r k using (63).Note that the controller has no feedforward signal.

Hyperparameter optimization
Extended DMD [27] with Tikhonov regularization is used to identify a Koopman model of the closed-loop system.The regularization coefficient β, along with the lifting function phase φ, are free parameters that must be chosen to minimize prediction errors on previously unseen datasets.Bayesian optimization [28] is used to automate the hyperparameter tuning procedure, but the cost function to use is an important design choice.In this section, two potential approaches are discussed: one using the closed-loop prediction error and one using the plant prediction error.
The closed-loop hyperparameter optimization approach consists of minimizing the prediction error of the closed-loop system using Bayesian optimization, and only extracting the plant system after the optimization procedure has completed.Using absolute error as the scoring metric results in the closed-loop cost function where the input to the closed-loop system is υ k = θ r k θr k T .Note that the model prediction errors are ∆x c k , ∆θ k , and ∆ θk , and are not to be confused with the controller tracking errors, θ e k and θe k .In some situations, minimizing the closed-loop cost function leads to an unnecessarily large regularizer, which degrades the performance of the extracted plant system.An improved hyperparameter optimization approach is to identify the closed-loop system, extract the plant, and evaluate its prediction error at each iteration of the Bayesian optimization algorithm.This approach requires a different cost function, that being where the plant input is υ p k = τ k .Given that the true quantity of interest is the prediction error of the plant in isolation of the controller, the cost function in (72) is preferred.In both (71) and (72), the absolute error of each state is computed, then averaged across the states to form the score.The mean squared error of each state could also be used to similar effect.In the results that follow, all hyperparameter optimization procedures minimize the prediction error of the plant rather than the closed-loop system.

Experimental results
In this section, the predicted trajectories of the identified Koopman system using the lifting functions (65) and (66) are compared with those of a linear system using the lifting functions (67) and (68).The identified systems models are compared both on their own and in the context of the complete closed-loop system.In each plot, the lifted states are retracted and re-lifted after every prediction iteration.As such, these errors follow the definition of local error from [29].

Closed-loop system
In this section, the predictions of the Koopman model and the linear model are compared in the context of the closed-loop system.Recall that the hyperparameters of both systems are identified using the open-loop prediction errors of the plant systems.The linear model of the closed-loop system is unable to generate the sustained vibrations found in the controller state and velocity signals, and instead settles around the mean value.In contrast, the Koopman model of the closed-loop system is able to predict these oscillations at the cost of increased transient error.to the input in Figure 5.The corresponding relative prediction errors are shown in Figure 8.The linear model accurately predicts the low-frequency dynamics of the closed-loop system, but it is not able to predict the high-frequency gearbox vibrations.The prediction instead settles near the mean value of the signal during constant-velocity segments of the trajectory.The Koopman model improves upon the linear model by predicting gearbox vibrations in addition to the low-frequency dynamics.This improvement is also reflected in the power spectral densities of the prediction errors shown in Figure 9.

Plant system
In this section, the predictions of the Koopman model and the linear model are compared after the corresponding plants have been extracted from their closed-loop systems.Figure 10 shows the predictions of the linear plant model and the Koopman plant model, where the motor torque is now the input.Due to the correlation between the plant's input and output introduced by the feedback controller, the torque input contains vibrations.As such, both the linear plant model and the Koopman plant model are able to predict high-frequency vibrations in the velocity.This is reflected in the similarities between the relative errors in Figure 11 and the power spectral densities of the errors in Figure 12.In fact, the Koopman plant model performs slightly worse than the linear model in an open-loop setting.
While both plant models predict velocity oscillations when subject to an oscillating torque input, only the Koopman system truly models the nonlinear behaviour of the plant.This fact is demonstrated in Figure 13, wherein both plant models are tested with a synthetic constant torque input.Only the Koopman model is able to predict the nonlinear gearbox vibrations when subject to an input that does not already contain them.

Conclusion
When identifying closed-loop systems, it is often not possible to neglect the effects of the control loop.This holds true for Koopman operator approximation as well as linear system identification.In this paper, a closed-loop Koopman operator identification method is presented, where the closedloop system is first identified as a whole, and then the Koopman system modelling the plant is extracted.The effectiveness of this method is demonstrated using a motor drive with nonlinear gearbox vibrations.While both a linear model and a Koopman model of the system yield similar prediction errors, only the Koopman system captures the nonlinear oscillations introduced by the gearbox dynamics.In future work, this plant Koopman model can be used to design an improved controller for the system.

Figure 2 :
Figure 2: Series interconnection of the controller and plant.

Figure 3 :
Figure 3: Feedback interconnection of the controller and plant.

Figure 4 :
Figure 4: Motor drive used to generate the training data, which consists of a motor with Harmonic Drive gearbox.The gearbox introduces nonlinear oscillations into the system, leading to tracking errors at specific frequencies.

Figure 5 :
Figure 5: Reference input to the closed-loop drive system from the first test episode.The velocity reference is a slightly smoothed pseudorandom binary sequence.It was generated by setting pseudorandom position checkpoints at the maximum allowable velocity and acceleration in the drive's control software.

Figure 6 :
Figure 6: Power spectral density of output velocity tracking error during a constant-velocity trajectory segment.The velocity at the gearbox input is 50 rev/s, leading to vibrations at integer multiples of 50 Hz.The most severe tracking errors occur at 50 Hz.A logarithmic scale is used to emphasize the high-frequency vibrations.

Figure 7 :
Figure 7: Predicted controller state, position trajectory, and velocity trajectory from the first test episode.The linear model of the closed-loop system is unable to generate the sustained vibrations found in the controller state and velocity signals, and instead settles around the mean value.In contrast, the Koopman model of the closed-loop system is able to predict these oscillations at the cost of increased transient error.

Figure 7
Figure7shows the trajectories predicted by the closed-loop linear and Koopman models subject

Figure 8 :
Figure 8: Controller state, position, and velocity prediction errors from the first test episode.The closedloop Koopman model's prediction results in lower steady-state velocity error compared to the closed-loop linear model.

Figure 9 :
Figure 9: Power spectral densities of the controller state, position, and velocity prediction errors from the first test episode.The Koopman model of the closed-loop system reduces controller state and velocity prediction errors at and below 50 Hz.A linear scale is used to emphasize the differences in the prediction errors of the models.

Figure 10 :
Figure10: Predicted position and velocity trajectories from the first test episode.Here, the torque signal is used as the input to the linear and Koopman plant models.Because of the feedback action of the controller, the torque signal has 50 Hz oscillations in it.As such, the predictions of both models oscillate at 50 Hz.

Figure 11 :
Figure 11: Position and velocity prediction errors from the first test episode.The prediction errors of the linear plant model and the Koopman plant model are similarly small, aside from the velocity oscillation of the Koopman model in the first segment of the velocity trajectory.

Figure 12 :
Figure 12: Power spectral densities of the position and velocity prediction errors from the first test episode.The prediction errors of the linear plant model and the Koopman plant model are indistinguishable aside from the reduced low-frequency velocity error of the linear model.A linear scale is used to emphasize the differences in the prediction errors of the models.

Figure 13 :
Figure 13: Predicted position and velocity trajectories for a synthetic constant torque input.The linear plant model does not oscillate at all, while the Koopman plant model predicts velocity oscillations due to the Harmonic Drive gearbox.
c) Recover plant system Figure1: Overview of the proposed closed-loop Koopman operator identification method.To simplify the presentation, no feedforward signal is used.(a) First, the controller reference, controller state, and plant state are recorded during a series of experiments.The controller state space matrices are known, and C p and D p are fixed.If the controller state is not directly available, it can be computed from its input and output.Only the Koopman matrix of the plant U p is unknown.(b) The plant state is lifted and augmented with the controller state and reference.The Koopman matrix of the closed-loop (CL) system, U f , is approximated.(c) The known controller state-space matrices are combined with the identified Koopman matrix to compute the Koopman matrix of the plant, U p .