Global, Unified Representation of Heterogenous Robot Dynamics Using Composition Operators: A Koopman Direct Encoding Method

The dynamic complexity of robots and mechatronic systems often pertains to the hybrid nature of dynamics, where governing equations consist of heterogenous equations that are switched depending on the state of the system. Legged robots and manipulator robots experience contact-noncontact discrete transitions, causing switching of governing equations. Analysis of these systems have been a challenge due to the lack of a global, unified model that is amenable to analysis of the global behaviors. Composition operator theory has the potential to provide a global, unified representation by converting them to linear dynamical systems in a lifted space. The current work presents a method for encoding nonlinear heterogenous dynamics into a high dimensional space of observables in the form of Koopman operator. First, a new formula is established for representing the Koopman operator in a Hilbert space by using inner products of observable functions and their composition with the governing state transition function. This formula, called Direct Encoding, allows for converting a class of heterogenous systems directly to a global, unified linear model. Unlike prevalent data-driven methods, where results can vary depending on numerical data, the proposed method is globally valid, not requiring numerical simulation of the original dynamics. A simple example validates the theoretical results, and the method is applied to a multi-cable suspension system.


I. INTRODUCTION
OBOTS perform complex dynamic tasks, interacting with the environment.As a legged robot interacts with a floor, its dynamics change discontinuously [1].As a manipulator robot manipulates an object, its fingers and arm are subject to non-holonomic geometric constraints, which may dynamically change [2].These robot dynamics are not unimodal; they are multi-faceted.Overall, their governing dynamics are represented as a combination of heterogenous equations of motion.When an independent set of generalized coordinates are used, the equations of motion are segmented into multiple regions, each representing specific dynamics subject to specific constraints in that region.Depending on contact conditions, constraints, and state locations, a different set of dynamic equations must be applied.This heterogeneous nature of governing equations poses a fundamental challenge for both dynamic modeling and control.The lack of a unified representation that is valid globally for all diverse state locations hinders theoretical analysis of global behaviors.
In system dynamics and control literature, Lyapunov functions have been used for guaranteeing the global stability of a particular class of hybrid systems.Methods have been developed for searching for or learning appropriate Lyapunov functions for given hybrid systems [3], [4].However, it remains a challenge to find a Lyapunov function for a general class of hybrid systems, especially robotics problems.Poincare's map has been applied to examine convergence and stability of periodic orbits, as observed in legged robots and others [5].This method requires computation of a series of points on a transversal section, which entails numerical simulations of the original trajectories.Furthermore, this method requires stability analysis of the resultant nonlinear system.Reachability has been studied extensively, and rigorous conditions for reachability have been established [6], [7].Although reachability is a fundamental property, more detailed analytic tools are required in robotics.Theoretical methodologies for hybrid systems are still limited and analysis is largely dependent on extensive simulations and numerical computations.A global, unified representation of a robot's governing equations could open the door to an alternative approach, which would supplement the existing theories and methodology.
The goal of the current work is to establish a methodology for obtaining a global, uniform representation of multi-faceted, heterogeneous robot dynamics that is amenable to analysis.We aim to construct a global, unified representation directly from the governing equations of heterogeneous robot dynamics.It is uniform in that a single model subsumes all the cases.No segmentation and switching are required.Such a unified representation facilitates the analysis of global behaviors, including stability and convergence.
Our approach utilizes a lifting of the dynamics using supernumerary state variables [8].Unlike the standard pointwise linearization, this lifting linearization has the potential to construct a global linear model.As illustrated in Fig. 1, governing equations in the space of independent state variables are nonlinear, heterogenous, and are segmented into local regions.As the system representation is lifted to a highdimensional space, it can be shown that the system behaves linearly and, more importantly, behaves uniformly with no segmentation nor switching.We exploit this linearity of representation for obtaining a unified representation.Although the governing equations are nonlinear and segmented in the original state space, they can be unified with a linear representation in a lifted space.All the heterogenous dynamics, including nonholonomic constraints and discontinuous state transitions among them, may be embedded into a single highdimensional linear equation.
The theoretical foundation of lifting dynamic modeling was established in the era of Littlewood [9] and Koopman [10] in the 1920's and 30's.Their results showed that a real, single valued function :  → , called an observable, which is compositional with a self-map :  → , denoted  ∘ , can be written as a linear transformation of the observable,  ∘  = .This linearity holds although the map F is a nonlinear function.Based on this, a nonlinear, autonomous system,  !"# = ( ! ) and its continuous time form, can be represented as a linear state equation in an infinite-dimensional space.This had a significant impact upon quantum mechanics in the 1930's and 40's.However, no practical application in the engineering field had been found until Igor Mezic pioneered an effective method for applying the Koopman operator to various fields of engineering problems three-quarters of century after the foundational work of Littlewood and Koopman [11].
Mezic and his peers have established the method based on spectral properties of the Koopman operator [12], [13].The Koopman operator theory has been integrated with Dynamic Mode Decomposition (DMD) and data-driven methods [14], [15].While the Koopman operator theory underpins the algorithm of DMD from the function analysis viewpoint [16], DMD provides the Koopman method with a practical datadriven technique.Based on the Koopman method, DMD has also been extended to a more general, effective algorithm.Furthermore, the Koopman-DMD approach has been extended to systems with control, or non-autonomous systems [17].It has opened the door to control system design, in particular, Model Predictive Control (MPC) [18].
Many successful applications of the Koopman-DMD approach have been reported.These include fluid mechanics [13], power systems [12], optimal control [19], and computer vision [20], [21].In robotics, too, the Koopman-DMD approach has made significant contributions in the last several years.Among others, Abraham and Murphy applied it to active learning for linearizing complex nonlinear dynamics [22].The method has also been applied to the modeling and control of soft robotic systems, where highly nonlinear, large deformations of soft materials have been modeled accurately [23], [24] and modeling of human-machine systems [25].
It must be noted, however, that a few major drawbacks and limitations exist to the DMD approach.First, it is difficult to guarantee global validity of a model, unless a "compact" set of data is available.Similar to other data-driven approaches, including machine learning, the results are dependent on data used for tuning the system.It is a challenge to collect all the data that guarantee the complete global behaviors of the system.Poor performance may result from insufficient data or from an inappropriate choice of observables.Thus, it is difficult to determine the root cause of performance issues, be it data collection or observable choice.
Another major limitation to the DMD approach is that the algorithm simply assumes that a linear dynamic relationship exists in a high-dimensional space [26].It does not address whether linearity exists or not.It is questionable particularly when the method is applied to heterogenous dynamical systems that possess discontinuities in state transition and switching of governing equations.
The method of this paper is an alternative to the DMD and other data-driven approaches based on DMD.A linear representation that is valid globally will be obtained directly from governing nonlinear dynamics, which are assumed available.In the following, the composition of an observable with a nonlinear state transition function  will be represented as a linear transformation of the observable function by using an integral kernel.A new formula called Direct Encoding will be developed to obtain the linear state transition matrix of the nonlinear system directly from inner products of observables and their composition with the nonlinear function .Conditions for applying the direct encoding formula to heterogenous nonlinear systems with discontinuity will be discussed.
Most prior works assume that the state transition map :  →  is smooth and continuous.This assumption is valid in many applications previously dealt with.However, the state transition map :  →  is not always smooth nor continuous when dealing with complex heterogenous dynamic equations.The goal of the current work is to establish a methodology for building a global and unified model of complex dynamical systems, which may include switching among diverse sets of dynamical equations and discontinuous transitions of state.
To this end, properties of the state transition map :  →  and a system of observable functions will be examined for finding conditions under which a linear state transition equation exists.Three major results will be obtained; a) a proper choice of observables matched with a given nonlinear state transition function  is required for obtaining a linear dynamic equation in a lifted space; b) if the conditions are met, a global linear model can be constructed through inner product computations over the state space; and c) although the conditions are not satisfied for a class of heterogenous systems, a finite-order model can be used to approximate the heterogenous dynamics under a specific mild condition.Because the method does not use simulation data for training a linear model, it does not possess the drawbacks that the DMD approach and data-driven methods suffer from.The method will be applied to a univariate hybrid dynamical system and a multi-cable heterogenous dynamical system.The implication of the proposed method and remaining open questions will be discussed at the end.

II. LINEAR REPRESENTATION OF NONLINEAR DYNAMICS
This section will present the basic formulation of the problem, including a minimum background of functional analysis and composition operators [27].Although the Koopman Operator theory has been established in a general Banach space [28], Direct Encoding entails inner products.As such, the following formulation assumes that all functions are involved in a Hilbert space.This results in a simple, straightforward transformation of nonlinear dynamics into lifted linear state equations using basic function analysis.This formulation based on inner products leads to the Direct Encoding formula and lays grounds for addressing the applicability of the method to heterogenous dynamical systems in the subsequent section.

A. Representation of a Compositional Operator Using a Kernel
The fundamental linearity given by the Koopman operator theory can be manifested with use of an integral kernel.Consider a discrete-time, autonomous dynamical system: where  !∈  ⊂ ℛ $ is an independent state vector of dimension n, ( ! ) is a real-valued vector function :  → , which is a self-map called a state transition function,  ! is output, and ( ! ) is an output function.In Koopman operator theory, ( ! ) is called an "observable".Combining (1) and (2) yields In this section we aim to examine the properties of [( !)] from a function analysis viewpoint.In doing so we treat [( !)] as a real-valued function on  ⊂ ℛ $ , and deal with a scalar function.In the following, we drop the subscript t to emphasize transformation of the entire function over  ⊂ ℛ $ , and de-emphasize time evolution.The connection between function transformation and time evolution will be established in the succeeding sub-section.
Consider a Hilbert space ℋ with an inner product 〈−, −〉: ℋ × ℋ → , expressed as where  # ,  % ∈ ℋ, and  % = represents the complex conjugate of  % .The function in (3) is a composition of function g with a selfmap F, which is represented using a composition operator as In the case of heterogenous dynamical systems, the state transition function F may be segmented into multiple regions, called "state locations", and the governing equations must be changed depending on the state location.Despite such complexities, there exists a linear relationship between observable function g and its composition with F -that is,  ∘  -as shown next using a kernel.

Proposition 1
Let ℋ be a Hilbert space on  ⊂ ℛ $ , and let { # ,  % ,  ' , ⋯ } be an orthonormal set of basis functions spanning ℋ.Let F be a state transition function :  →  and  be a function in the Hilbert space.

𝑔 ∈ ℋ
Then, there exists a linear relationship between function  and its composition with F, where  is an integral variable, and :  ×  →  is a kernel given by

Proof
From ( 6),  can be expanded in the orthonormal basis By using this expansion, the compositional function  ∘  can be written as Distributing the compositional operator yields From the definition of inner product (4) Note that the kernel (, ) encodes the state transition function () with the basis functions { # ,  % ,  ' , ⋯ }.The kernel does not depend on the function  ∈ ℋ, and it can be applied to arbitrary functions in ℋ: Concatenating these in an infinite dimensional vector of functions  + ∘  yields These observables  # ,  % , ⋯collectively represent the state of the system in a lifted space of infinite dimension.Denoting the concatenated observables by we aim to obtain a linear state transition equation.Under certain conditions, such a linear state transition equation exists, as discussed next.

B. Linear State Transition
While the kernel representation given by eq.( 7) manifests the linear transformation of observables to their compositions with the nonlinear state transition function , it is not in the form of linear state equation, which we need in engineering applications.The following Proposition provides the conditions for the existence of a linear state equation.

Proposition 2
Let { # ,  % ,  ' , ⋯ } be an orthonormal set of basis functions in a Hilbert space ℋ on  ⊂ ℛ $ , and let F be a state transition function :  → .If the compositional function  + ∘  is in the Hilbert space, then there exists a linear state transition matrix A of infinite dimension that transfers a lifted state in the form of , where { # ,  % ,  ' , ⋯ } forms an orthonormal basis in ℋ, to the lifted state associated to the state transition function F, where , and the state transition matrix A is given by
A special case of Proposition 2 is to use { # ,  % ,  ' , ⋯ } for the second set of orthonormal basis functions { # ,  % ,  ' , ⋯ }.In this case, the matrix A is simplified.
In Proposition 2, the observables { # ,  % ,  ' , ⋯ } must be an orthonormal set of basis functions.While this condition is restrictive, it can be relaxed.Let { # ,  % ,  ' , ⋯ } be an independent and complete set of observables spanning the Hilbert space.Using the orthonormal basis functions, each observable can be expanded to  + = ∑〈 + ,  ( 〉 ( , or collectively expressed as where Because both { # ,  % ,  ' , ⋯ } and ( # ,  % ,  ' , ⋯ ) are independent and complete basis functions, the matrix C is nonsingular.Thus, the following corollary allows us to use merely independent basis functions.

Corollary 2
In Proposition 2, if { # ,  % ,  ' , ⋯ } are an independent and complete set of basis functions, the state transition matrix A in (17) is given by where  ̅ has been given by (22) and  by (24).

III. MAIN RESULTS
This section presents a new formulation of the Koopman operator through the use of inner products.The linear representation is globally valid and does not depend on data.
Its applicability to heterogenous dynamical systems will be discussed.

A. The Direct-Encoding Method Using Inner-Products
Incorporating Proposition 2 and Corollaries 1 and 2 and constructing inner products, now we obtain the following main proposition.

Proposition 3
If observables { # ,  % ,  ' , ⋯ } form an independent and complete set of basis functions, spanning a Hilbert space ℋ, and their compositions with a state transition function :  →  are involved in the Hilbert space: then the following state transition equation exists with a lifted state and a state transition matrix given, respectively, by where the matrices  and  consist of inner products given by
Post-multiplying a row vector function [ # ,  % ,  ' , ⋯ ] to (27) and integrate over X yields, where the matrix  is factored out.Each integral in the vector outer product is the inner product of the corresponding two functions.Therefore, Since { # ,  % ,  ' , ⋯ } are independent, matrix R is nonsingular.Therefore, (28) is obtained.Q.E.D.
Note that computation of  =  .#does not require { # ,  % , ⋯ } nor  and  .#involved in (25).Once we know the existence of ( 27), matrix A can be computed directly by taking inner products between  + and its composition with .The inner products 〈 + ,  -〉 and 〈 + ∘ ,  -〉 exist and are welldefined because both  + and  -∘  belong to the Hilbert space.
Note that this formula for obtaining a linear system matrix  is an alternative to the prevailing Koopman DMD method, which assumes the existence of a linear system matrix  in data and computes the matrix by Least Squares Estimate.The  matrix based on the Direct Encoding is globally valid, while the DMD and its variants are data dependent.

B. Weak Existence Conditions for Finite-Dimensional Approximation
In Proposition 2, it is assumed that the compositional function  + ∘  is in a Hilbert space for all i, including ∞.As will be detailed in the succeeding section, this condition is difficult to satisfy when applying the compositional operator method to heterogenous systems, specifically to hybrid systems with discontinuous state transition.This section explores the possibility to apply the compositional operator method to hybrid systems by replacing the strong condition,  + ∘  ∈ ℋ, ∀  ∈ ℕ, by its weaker condition,  + ∘  ∈ ℋ, ∀  < ∞.Namely, the compositional function  + ∘  must be in a Hilbert space only for an arbitrary, finite i.This implies truncation of the infinite-dimensional lifted space, resulting in approximation of the linear state transition.From Proposition 2 and Corollary 1, we obtain:

Proposition 4
Let { # ,  % ,  ' , ⋯ } be an orthonormal set of basis functions in a Hilbert space ℋ on  ⊂ ℛ $ , and let F be a state transition function .If the first m compositional functions  + ∘  are in the Hilbert space, then there exists a linear state transition matrix  /) that transfers an infinite-dimensional lifted state in the form of , to its finite-dimensional state consisting of composite functions.
where  /) is a m by ∞ matrix given by TMECH-08-2022-14068.R1 Proposition 4 shows that, with a finite number of compositional functions being in the Hilbert space, the one-step time-evolution of the finite number of lifted state variables can be predicted exactly.Since (36) predicts only the finite number of state variables, the prediction of consecutive state transitions may not be exact.Nonetheless, truncated state variables and system matrix may yield reasonable accuracy for a large m.

IV. AN ILLUSTRATIVE EXAMPLE
Consider a univariate nonlinear dynamical system on  ∈  = [0,1].Note that this piecewise linear system has a discontinuity at  =  * .It is heterogeneous and a type of hybrid system.
From (8) the kernel (, ) that encodes this state transition equation () is obtained by using the following exponential trigonometric functions, which form orthonormal basis functions spanning a Hilbert space.
The kernel should be able to map an arbitrary observable in the Hilbert space to its composition with ().For this kernel to be transformed to the time-evolution, state transition matrix A, the conditions of Proposition 2 must be satisfied.In particular, Using the exponential trigonometric function, we obtain: This compositional function can be expanded in the exponential trigonometric basis for an arbitrary finite k.Despite the discontinuity due to the jump in (), function  ( ∘  satisfies the following four conditions required for the existence of Fourier series [29]: 1)  ( ∘  is absolutely integrable.
2) There are a finite number of discontinuities within a finite interval.3) All the discontinuities are finite.4)  ( ∘  is of bounded variation; there are no more than a finite number of maxima and minima in any finite interval.With the existence of a Fourier series,  ( ∘  is in a Hilbert space for any finite k.This allows us to construct the statetransition matrix based on Proposition 4.
Suppose that the following real-valued observables are selected as an independent set of basis functions, which are complete, Conversion of  ( to  -is given by a nonsingular matrix C and its inverse  .#, too, can be computed easily.See Appendix A. These matrices can be truncated to  ×  matrices, and the time-evolution state transition matrix  / =  /  ̅ /  / .# is obtained. This  / matrix should be able to predict the time-evolution of the lifted state  != [ # ,  % , ⋯ ,  / ] , mapped to  !"# .Fig. 3 shows the comparison the state transition based on the linear model,  !"# =  !, to the true transition through the original nonlinear function  !"# = ( ! ).Due to the discontinuity in (), the true time evolution of the system shows a series of pronounced jumps, as shown by the solid lines in the figures.Nonetheless, the linear model can track this response.For small numbers of observables, truncated at m = 17 and 33, the prediction errors are significant.However, as the number of the observables increases, m = 1025, the linear model can almost perfectly predict the complex response.Fig. 3-(e) shows mean squared error vs. the number of observables.It is confirmed as the number of observables increases the error approaches zero in the mean.With finite m, this is still an approximation.However, the accuracy is high enough for practical applications.Note that the existence conditions hold for any large m that is finite.This linear model is global; it is valid for the entire domain.Unlike data-driven methods, where results depend on how data are collected, this model is derived directly from the governing model.This model is also uniform.Despite the discontinuous switching of the two linear models involved in the system, the linear model represents the heterogeneous system behavior with a single system of linear dynamic equations uniformly in the lifted space.
This linear model facilitates dynamic analysis of the hybrid system.Fig. 4 shows the plot of poles obtained from the  %34 matrix.that all the poles are on or within the unit circle.Clearly, the system is marginally stable.A group of poles are on the unit circle, indicating that the system has undamped oscillation modes.This agrees with the time trajectories in Fig. 3, where the system converges to a periodic orbit.

V. MULTI-CABLE DYNAMIC MANIPULATION
Now the method is applied to a robotics problem.Contactnoncontact switching of dynamics is a challenging problem but is a fundamentally important issue in broad fields of robotics.In multi-cable manipulation, for example, an object suspended with cables experiences discrete switching of governing equations as each cable alters between slack and taut conditions [30].As shown in Fig. 5, a mass m is connected to a pair of cables.When the point mass is released from a certain height, it drops like a projectile object, but it bounces back when the cable(s) become taut, like a ball bouncing on a floor.Each cable can bear only a unidirectional load; it goes slack when a compressive load acts on the cable.The constraints due to this unidirectional nature of cable suspension are therefore nonholonomic.Depending on which cable is taut or slack, different dynamic equations are switched.
Let  / and  / be coordinates of the point mass.Equations of motion are segmented to 4 regions.
where  5 and  6 are tensions of cables A and B, respectively,  5 and  6 are unit vectors pointing in the direction of the two cables, and  # , ⋯ ,  2 are regions of the point mass locations corresponding to D1: both cables are slack, D2: cable A is taut and cable B is slack, D3: cable B is taut and cable A is slack, and D4: both cables are taut.Both cables have a high extensile stiffness but zero stiffness for a compressive load.The cables are assumed to be massless and to have small damping.The independent state variables are: The dynamic range of the system is finite, given the length of each cable.Gaussian Radial Basis Functions (RBFs) are used as observables: where in each RBF  ( ,  ( is the center point.By construction, RBFs are guaranteed to be independent [31].The number of RBFs, that is, the order of the lifted space, must be finite for implementing the observables.In total, 1,500 RBFs are used for covering the 4-dimensional state space within the bounded dynamic range.They are placed at optimal locations by clustering sample points using the k means ++ algorithm.More details are described in Appendix B. The independence of the RBFs can be confirmed by examining whether the R matrix in ( 29) is non-singular.The Q matrix is determined by computing the inner products between individual RBFs and their composition with the state transition function (): From these R and Q matrices, the A matrix is obtained:  =  .# .Fig. 6 shows the comparison between the point mass trajectory of the original nonlinear model -ground truth -and the one predicted with the linear, lifted model -prediction.In this comparison, the point mass is released from an initial position where both cables are slack.The point mass bounces as one of the cables becomes taut and immediately goes slack.
With some damping at the cable, the converges to an equilibrium position after a series of bounces, as shown with a black solid line in Fig. 6-(a).This nonlinear dynamic response is predicted almost perfectly with the linear model in the lifted space shown in red. Figure 6-(b)-(e) show the trajectories of the individual state variables.The predicted trajectories shown in red agree with the ground truth black lines even after many bounces.Comparing to the prior work on multi-cable manipulation [30], where a data-driven method based on DMD was used, the proposed method shows an order-of-magnitude more accurate prediction, although the system structure is different.This linear model incorporates the four dynamic equations into a single unified representation.Again, the results are obtained directly from the governing equations, and the model is valid globally.

VI. DISCUSSION: APPLICABILITY OF THE KOOPMAN OPERATOR
TO HETEROGENEOUS DYNAMICAL SYSTEMS The current work addresses key conditions in applying the composition operator method to those heterogenous dynamics problems.To find such conditions, the current work took a simple, straightforward approach to converting function transformation equations to time-evolution dynamic equations based on composition operator theory.The conversion hinges on the conditions given by ( 16) in Proposition 2 and by (26) in Proposition 3.These conditions are also the key to examining the applicability of the Koopman Operator to a class of challenging robotics problems.If one uses the Koopman Operator for merely linearizing a nonlinear dynamical system, where the state transition function () is smooth and continuous, these conditions, ( 16) and (26), are satisfied for most observable functions,  ( ∈ ℋ.However, this is not automatically met for heterogeneous dynamical systems, including hybrid and switched systems.
In the proof of Proposition 2, the compositional function  -∘  was expanded to  -∘  = ∑ 〈  -∘ ,  ( 〉 ) (*#  ( .Through this expansion, observables  # ,  %,⋯ were pulled out.This is the key mechanism for creating a new extended state space with state variables  = [ # ,  % , ⋯ ] , , leading to the formation of matrix A. In Proposition 2, it was assumed that { # ,  %,⋯ } forms an orthonormal basis, but this condition was relaxed in Proposition 3, where  # ,  %,⋯ are an independent and complete set of basis functions. These existence conditions are dependent on the selection of observables as well as on the properties of state transition function ().The first example discussed above is a univariate heterogenous dynamical system with discontinuities in state transition.For this system, properties of the compositional function,  + ∘ , can be examined analytically.Although discontinuous transitions are involved, a Fourier series exists for an arbitrary j that is finite.This allows us to expand  + ∘  to ∑ 〈  + ∘ ,  ( 〉 ) (*#  ( and obtain the time-evolution dynamic equation.This does not necessarily mean lim +→)  + ∘  ∈ ℋ because the compositional function may have an infinite number of discontinuities in a finite interval, which is against the conditions for the existence of Fourier series expansion.Nonetheless, for any finite order  < ∞,  + ∘  ∈ ℋ.This property supports the high accuracy prediction that is achieved with a large number of observables, as confirmed in Fig. 3-E.
For practical implementation, the infinite-dimensional observable vector  must be truncated.In doing so, the condition: is a useful property that underpins high-accuracy, highdimensional truncation.Relaxing the conditions for the existence of a time-evolution A matrix to the above condition, ∀ < ∞, allows us to apply the composition operator method to broader systems, including hybrid and switched systems with discontinuities.In fact, the prediction based on the  / matrix under the conditions (35) showed an extremely high accuracy over a long period of time, as demonstrated in Fig. 3.
In general, it is difficult to analytically examine  + ∘  ∈ ℋ for more general, multivariate systems with complex heterogeneous dynamics.However, there are a few methods for numerically examining the conditions [27], [31].
• Using these metrics, we can examine the conditions.
The second example, the multi-cable suspension problem, is a multivariate, heterogenous dynamical system described with 4 sets of dynamic equations.This system does not have any discontinuity in state transition; both position and velocity of the mass are continuous, as long as the stiffness of the cables is finite.The 4 dynamic equations have been combined and transformed to a single linear dynamic equation, which has been shown extremely accurate.The linear model can predict the mass trajectory with high accuracy even after many bounces.The high accuracy modeling was achieved through direct encoding of the governing equations, rather than curve fitting to numerical data.
As the extensile cable stiffness increases to infinite stiffness, the velocity of the mass tends to change discontinuously when bouncing.To deal with discontinuous state transitions, the method based on Proposition 4 is required, as in the case of the univariate example.This entails truncation of the lifted space dimension, which leads to approximation.However, truncation is necessary for practical implementation.Therefore, this does not degrade the usefulness of the composition operator approach and its applicability to heterogenous dynamical system with discontinuous state transition.
Overall, the proposed composite operator method for global, unified representation of heterogenous robot dynamics is applicable to various robotics problems, including hybrid and switched systems with discontinuities.
A few caveats must be stated, however.First, a large number of observables are required for dealing with discontinuous state transitions.This can be alleviated by using observable functions that are locally tunable.RBFs, for example, can allocated densely to a region near the discontinuity [31].See Appendix B. Second, the Direct Encoding method is difficult to apply to a high-order nonlinear system due to the multivariate inner product integration.As the number of independent state variables, n, exceeds 6 or 7, it is difficult to perform the computation of inner products.This necessitates an effective computational algorithm.

VII. CONCLUSION
The current work demonstrates that complex robot dynamics given by heterogenous governing equations can be represented with a single, unified linear model in a lifted space.Although the original dynamics are highly nonlinear, segmented, and heterogeneous, the lifted linear model is valid globally.Although the system experiences a multitude of jumps and transitions to diverse dynamic behaviors, the prediction accuracy remains high.In robotics, heterogeneous dynamics are common, due to contact-noncontact transition, non-holonomic constraints, and other discrete transitions among diverse dynamics.The composition operator approach presented in the current work can be a promising methodology in dealing with those heterogeneous dynamics.
In the current work, only autonomous systems having no exogenous input are considered.Extending the original Koopman Operator theory to general control systems is a critically important challenge that has been addressed by several authors and is still an active research area in the field [17], [19], [32], [33].
The theory established for autonomous systems, including the current work, can be extended without any significant change to some class of non-autonomous systems, if they meet the following conditions.
• Control  !comes into a system as a linear term:  !"# = ( ! ) +  ! .In this case the effect of the input term  !can be separated, and the matrix B can be determined by solving a least squares estimate problem [8], [17], [33] ; or • If the control is a feedback regulator,  != ℎ( ! ), then the problem reduces to an autonomous system by embedding the control law into the system [17], [19].These are simple cases to which the existing methods for autonomous systems can be applied.For other general cases, the Koopman methods, including the current work, must be modified substantially.An elegant solution to this problem is to represent a control system as a bilinear system, where inputs are multiplied with functions of state [34].The resultant system is not linear, but the Koopman operator can be applied more rigorously.More research effort is needed to extend the current work to those general control cases.
Although the general non-autonomous system's problem has not yet been solved, it has been reported that the Koopmanbased lifting linearization is effective for controlling a class of nonlinear systems.The multi-cable system described in Fig. 4 has been controlled with MPC [30].The highly nonlinear, heterogenous dynamics are represented with a single linear state equation, which allows for efficient real-time MPC computation.Unlike nonlinear MPC, the linear MPC is a convex-optimization computation with no local minima problem.Reference [35] reported that the computation time reduced approximately 100 times compared to nonlinear MPC.
Finding an effective set of observables remains a challenging problem faced by the Koopman operator research community [36].In particular, the proposed method entails effective observables that meet the critical condition  + ∘  ∈ ℋ, which depends on both the properties of the observables and the heterogenous dynamics of the system.Effective algorithms must be developed for searching for observables that can approximate the system with fewer variables and that can keep the compositional function in a Hilbert space despite the heterogenous nature of the dynamics.
Despite the limitations to the current method, the proposed method is a promising approach to dealing with challenging robotics problems delineated as heterogeneous dynamics.The proposed method has the potential to provide the robotics community with a powerful tool: global, unified representation of otherwise complex heterogenous robot dynamics, which is also amenable to analysis and synthesis.The RBFs must be tuned properly so that a finite number of RBFs can approximate the system effectively.The most critical is the locations of center points  ( .A clustering technique was applied to find effective locations of the center points. The overall prediction accuracy of the linear model, ( + 1) = (), was evaluated by the predicted trajectories against the ones computed from the nonlinear model,  !"# = ( ! ).In doing so, several ground truth trajectories were created.Using these as sample points, we can find where RBFs should be placed densely or sparsely.We applied the k means clustering algorithm with an improved initial setting based on the k mean ++ algorithm.Fig. 7 shows the optimized locations of the center points.High-density center points can be found in the vicinity of the locations where one of the cables switches between taut and slack conditions.We specified the number of clusters to be 1,500 for placing 1,500 RBFs.

Figure 2 -
shows this state transition function with parameters, , , , and  * .This system consists of two distinct linear dynamics in two segmented regions or state locations.Note that this piecewise linear system has a discontinuity at  =  * .It is heterogeneous and a type of hybrid system.From (8) the kernel (, ) that encodes this state transition equation () is obtained by using the following exponential trigonometric functions, which form orthonormal basis functions spanning a Hilbert space.

Fig. 2 -
(b)~(e) confirm this property.All of the observables,  # = ,  % =  % ,  ' =  ' ,  2 = cos 2, are perfectly transformed to their compositional functions with the same kernel.Namely, the right-hand side of (7) computed with the kernel agrees with the true compositional function of the left-hand side.

Fig. 2
Fig. 2 Example of univariate heterogenous dynamical system (a) The state transition function is heterogenous, consisting of two linear systems combined at x*, where a discontinuous jump occurs.(b)~(e): These plots show that the Koopman operator can transform an arbitrary observable in the Hilbert space to its composition with F(x).The linear transformation shown with black circles can reproduce the correct nonlinear function g[F(x)] shown in red.Four different observables are tested; (b)  != , (c)  " =  " , (d)  # =  # , and (e)  $ = 2.

JumpFig. 3
Fig. 3 Evaluation of the linear lifted model in terms of timeevolution prediction accuracy Black lines in (a)~(d) show the true time trajectory of the hybrid dynamical system, compared to the ones predicted with the linear lifted model shown in blue.(a) With 1,025 observables, almost perfect prediction was achieved, while the prediction accuracy declines with 257 observables (b), 33 observables (c), and 17 observables (d).(e) summarizes the prediction accuracy in terms of root mean square error for different numbers of observables.The error converges to zero in the mean.

Figure 4
Figure 4 Plot of poles of  "%& on complex plane

Fig. 5
Fig. 5 Diverse states of a point mass suspended with two cables

Fig. 6
Fig. 6 Accuracy of time evolution prediction for the point mass suspended with two cables.(a) Trajectories in the xy-plane.The black solid line shows the trajectory generated with the original nonlinear, heterogenous dynamic equations, while the red dash-and-dot line shows the trajectory predicted with the linear lifted model.The two trajectories overlap, indicating high accuracy of prediction.(b) ~(e) show the time profiles of position x and its velocity and position y and its velocity, respectively.