Soft robots incorporate non-rigid materials into their morphology to facilitate compliant interactions with the external world. This compliance allows them to manipulate delicate objects, adapt to unstructured environments, and interact safely with coexisting humans [1, 2, 3]. Since their utility is derived from their compliance, control methods that preserve and/or exploit this property are desirable.
Accurate models facilitate better control performance. When an accurate model is available, predictive controllers can be built by using the model to calculate a feedforward term, then adding a feedback term to account for minor model uncertainty and disturbances. If an accurate model is unavailable, feedback must be relied upon more heavily. This poses several problems for soft robots. First, feedback requires sensing, but the morphology of soft robots precludes the use of most conventional sensors. Suitable alternatives are currently in development [4, 5, 6, 7], but are not yet readily available. Second, relying heavily on feedback to compensate for an inaccurate model has been illustrated to reduce the compliance of soft robotic systems . That is, excessive feedback negates the desirable compliance of a soft robot by replacing its natural dynamics with those of a slower, stiffer system. Therefore, accurate models are required to control soft robots in a manner that reduces dependence on feedback and simultaneously preserves compliance.
Models for soft robots can be separated into two categories: physics-based models and data-driven models. Physics-based models are constructed from observations of component material properties and first-principles, while data-driven models are constructed from observations of system behavior. Physics-based models have the ability to make predictions about a system’s behavior before the system is constructed. Thus, they are often used to inform the design and construction of soft robots intended for particular tasks [9, 10, 11, 12, 13]
. However, the infinite degrees of freedom and hysteretic behavior of soft materials make it difficult to construct accurate physics-based models of soft robots without making simplifying assumptions such as constant curvature[14, 15], quasi-static [16, 17, 18], or simplified geometry [19, 20, 21]. These models only describe behavior well in the subset of robot configurations where the assumptions hold; hence, they are limited in applicability.
Data-driven models are more broadly applicable since they do not make structural assumptions about the system. Data-driven models are constructed from observations of system behavior rather than from first-principles. Hence, given sufficient data, they are capable of describing system behavior well over the entire range of observations. By virtue of their soft bodies, it is often possible to safely command arbitrary control inputs without risk of damaging a soft robot or nearby human operators. Soft robots are therefore amenable to data collection over their entire operation range, making them particularly well suited for data-driven system identification techniques.
To capture the characteristic nonlinear behaviors of most soft robots, identification of a nonlinear model is necessary. Unfortunately, identifying a nonlinear model from data typically consists of solving a nonlinear, non-convex optimization problem, for which global convergence is not guaranteed . Furthermore, most nonlinear system identification methods require the manual initialization and tuning of training parameters, which have an obscure impact on the resulting model. A neural network, for example, may be able to capture the nonlinear behavior of a soft robot 
; however, its accuracy depends on the number of hidden layers, number of nodes per layer, activation function, and termination condition used during training, which must be selected through trial and error until acceptable results are achieved.
Linear model identification, on the other hand, does not suffer from the typical shortcomings of nonlinear identification since linear models can be identified via linear regression . However, linear models are poorly suited to capture the behavior of soft robots since their characteristic behavior is distinctly nonlinear .
In this work, we employ Koopman operator theory to generate nonlinear, data-driven models of soft robots via linear regression (see Fig. 1). This approach relies on the idea of lifting nonlinear dynamical systems to an infinite dimensional space, where those systems have a linear representation. In this space, it is possible to describe the dynamic behavior of a system by a linear operator rather than a nonlinear vector field . This linear operator, called the Koopman operator, is identified via linear regression , so it does not suffer from the convergence and tuning problems that are characteristic of neural networks and other nonlinear system identification methods.
Our primary contribution is demonstrating, on a real system, a data-driven method for constructing globally valid nonlinear models that does not require the manual tuning of multiple training parameters. To do so, we apply the Koopman based system identification method described in  and  to create a dynamic model of a soft robot arm and verify that it captures the system’s true dynamic behavior better than the models generated by several other state-of-the-art nonlinear system identification methods including a neural network, a nonlinear auto-regressive with exogenous inputs (NLARX), a nonlinear Hammerstein-Wiener, and a state space model. Koopman operator theory has previously been used to perform model-based control of nonlinear systems . However, to the author’s knowledge, this technique has never before been used to identify the dynamic model of a real system, much less a soft robot. We believe that this system identification method applied to soft robots will enable the rapid development of new and effective control strategies by making accurate nonlinear dynamic models easier to construct.
The rest of this paper is organized as follows: In Section II we formally introduce the Koopman operator and describe the system identification method. In Section III we describe the soft robot and experimental procedure used for collecting input-output data of the system. In Section IV we summarize the results of applying various nonlinear system identification techniques to the collected data and compare the performances of the models generated. In Section V concluding remarks and perspectives are provided.
Ii Koopman Operator Method for System Identification
The system identification method utilized in this work exploits the fact that any finite-dimensional nonlinear system has an equivalent infinite-dimensional linear representation in the space of real-valued functions of the system’s state and input. In this space of real-valued functions, the (linear) Koopman operator describes the flow of functions along trajectories of the system. The relationship between the finite and infinite dimensional representations of the system is bijective and well-defined . This enables us to approximate the Koopman operator via linear regression on observed data, then extract the equivalent nonlinear system representation. The remainder of this section summarizes the system identification method presented in  and  applied to a system with known input, which is later employed and validated on a real soft robotic system.
Ii-a Overview of Lifting Technique
Consider a dynamical system
where is the state of the system, is the input and is continuously differentiable in . Denote by the solution to (1) at time when beginning with the initial condition at time and a constant input applied for all time between and . For simplicity, we denote this map, which is referred to as the flow map, by instead of .
The system can be lifted to an infinite dimensional function space where and are compact subsets and is the space of square integrable real-valued functions with domain . Elements of are called observables. In , the flow of the system is characterized by the set of Koopman operators , for each , which describes the evolution of the observables along the trajectories of the system according to the following definition:
As desired, is a linear operator even if the system (1) is nonlinear, since for and
One can show that there is a one-to-one correspondence between the infinite-dimensional Koopman operator and finite-dimensional vector field. To understand this relationship, consider the derivative of an observable, , along trajectories of the system:
where the second relation follows since is held constant and where is the gradient with respect to . Since this equation holds for all observables, we introduce the infinitesimal generator of the Koopman operator [30, Equation 7.6.5] which is defined in terms of the vector field as
The infinitesimal generator thus describes the dynamics of the observables along trajectories of the system. Recalling that the Koopman operator describes the flow of observables, one can show that the Koopman operater is associated with its infinitesimal generator via the following relation:
By providing a linear representation of a system with a one-to-one correspondence to its nonlinear representation, Koopman operator theory enables linear system identification of nonlinear systems. This process proceeds in three steps which are summarized by Algorithm 1. In step one, measured states of the system are lifted to the space of observables. Step two consists of performing least-squares regression on the lifted data to obtain an approximation of the Koopman operator. In step three, an approximation of the nonlinear vector field is obtained via (7) and (8). The following three subsections describe each of these steps in more detail.
Ii-B Step 1: Lifting the Data
The first step of the Koopman-based system identification method consists of converting empirical data into a form that can be used to identify a linear model in the space of observables. Theoretically this process would consist of “lifting” state measurements into the infinite-dimensional space of observables . To be implementable, however, measurements can only be lifted into a finite-dimensional subspace. Define to be the subspace of spanned by linearly independent basis functions (e.g. monomials, sinusoids, exponentials). Any observable can be written as a linear combination of elements of its basis
Note that the vector of coefficients provide a vector representation for any . To evaluate at a given state and constant input , we introduce the lifting function defined as:
Then, can be expressed in vector form as
We refer to as an -dimensional “lifted” version of , since multiplying by the vector representation of an observable yields the value of the observable at .
Ii-C Step 2: Approximating the Koopman Operator
The second step of the Koopman-based system identification method is to identify the Koopman operator that best describes the flow of the lifted versions of measured data points. While the Koopman operator is theoretically infinite-dimensional, for practical purposes we identify a finite-dimensional approximation of it in which we denote by . Note, that can be represented by an matrix which operates on observables via matrix multiplication:
where are vector representations of observables in . Our goal is to find a that describes the action of the infinite dimensional Koopman operator as accurately as possible in the -norm sense on the finite dimensional subspace of all observables. Therefore, to perfectly mimic the action of acting on an observable in , the following should be true
Since this is a linear equation, it follows that for a given , solving (13) for yields the best approximation of on in the -norm sense:
where superscript denotes the least-squares pseudoinverse.
To approximate the Koopman operator from a set of experimental data, we take discrete state measurements with sampling period . Under the assumption that the control input is constant between samples, we separate the data into a set of so-called “snapshot pairs” of the form where
and denotes measurement noise. For our basis of , we choose the basis of monomials of and with total degree less than or equal to , which implies [27, Section III]. We then lift all of the snapshot pairs according to (10) and compile them into the following matrices:
is chosen so that it yields the least squares best fit to all of the observed data, which, following from (14), is given by
Ii-D Step 3: Obtaining the Vector Field
The final step of the Koopman-based system identification method is to identify the nonlinear vector field by making use of the one-to-one correspondence between the infinite and finite dimensional system representations. As earlier, our goal is to find an that describes the behavior of the vector field as accurately as possible in the -norm sense on the finite dimensional subspace .
The vector field is related to the Koopman operator through its infinitesimal generator according to Equation (7). With the approximation of the Koopman operator found in Section II-C, we can solve for the infinitesimal generator of the set of Koopman operators by inverting (8):
where denotes the principal matrix logarithm [31, Chapter 11]
. Note that the principal matrix logarithm is only defined for matrices whose eigenvalues all have non-negative real components, and thatmay have zero or negative eigenvalues when the number of data points is too small . Therefore this method might fail if the number of data points is insufficient. In this instance, more system measurements can be taken to resolve the issue.
We seek the vector field such that (20) holds as well as possible in the -norm sense for all observed data. Therefore we choose the least-square solution to (20) over the set of all observed data points which is given by
Iii System Identification of a Soft Robot
To demonstrate and evaluate the performance of the method outlined in Section II, we applied it to a continuously deformable soft robot arm and compared the resulting model to those generated by other nonlinear system identification techniques. In the following, we explain in detail the used robotic hardware, experimental setup, measurement procedure, and the data processing involved in the system identification process, as well as the process by which performance was evaluated and compared across models.
Iii-a Hardware Description
The soft robot used in this experiment consisted of three pneumatically driven McKibben actuators (also known as Pneumatic Artificial Muscles or PAMs) adhered together by latex rubber and connected to a common base mount on one end and to an end effector on the other (see Fig. 2). During the trials, the pressure inside each actuator was varied using a pneumatic pressure regulator (Enfield TR-010-g10-s), and the displacement and velocity of the end effector was measured at using a commercial motion capture system (Phase Space Impulse X2E).
For this robot, it is our primary interest to control the motion of the end effector. Hence, we chose a state representation of the system which is capable of describing the dynamics of the end effector as an ordinary differential equation, namely the position and velocity of the end effector with respect to a global coordinate frame, as shown in Fig.2
Iii-B Randomized Control Input
To generate a representative sampling of the system’s behavior over its entire operation range, a randomized input was applied. The control input into the system was a set of three signals into the pressure regulators corresponding to actuator pressures of
Before each trial, a table
of uniformly distributed random numbers between zero and ten was generated to be used as an input lookup table. Each control input was smoothly varied between elements in consecutive columns of the table over a transition period, with a time offset of between each of the three control signals
where is the current index into the lookup table at time .
Iii-C Data Collection
Data collection proceeded in 4 trials, each lasting between 8 and 40 minutes with an input transition period between 3 and 8 seconds. After collecting data, raw position and velocity measurements were put through a moving average filter with window size of s to reduce noise, then sampled uniformly with a period . Sampled velocity measurements were put through a second moving average filter with window size of s, due to the higher noise content of the velocity signal. The time-series data from each trial was partitioned into training and validation sets. Three 10 second validation sets were extracted from each trial and the remainder of the trial data was used for training.
Iii-D Model Comparison
We generated a state space model from the technique described in section II using a monomial basis of maximum degree . We then evaluated its accuracy by comparing model simulations to validation data (Fig. 3). Goodness of fit for the trajectory of a state was calculated using the normalized root-mean-square error (NRMSE), defined:
where is the simulated value of the state, is the total number of points, and are the measured minimum/maximum values of the state observed over all trials.
The performance of our model was benchmarked against a linear state space, nonlinear Hammerstein-Wiener, nonlinear auto-regressive with exogenous inputs (NLARX), and a neural network model. All models were trained and evaluated on the same training and evaluation data and generated using either the Matlab System Identification Toolbox or Neural Network Toolbox . The state space model was generated using the subspace method [24, Chapter 7] and specified to be 6 dimensional, i.e. the same dimension as the state defined in 22
. The neural network model was trained using the Levenberg-Marquardt backpropogation algorithm and sigmoid activation functions. It was trained several times using combinations of 10-30 hidden neurons and 1-10 delays. Only the results for the best of these models, corresponding to 10 hidden neurons and 10 delays, is displayed in Fig.4 and Table I.
The model generated by the Koopman system identification method has a total RMSE averaged across position states and velocity states of 5.98 mm and 3.66 mm/s, respectively. As shown in Table I, this corresponds to a total NRMSE averaged across all states of 2.1%. By this metric, it performs more than twice as well as the best competing linear and nonlinear models which have average NRMSEs of 4.6% and 4.5%, respectively (see Fig. 4
). The Koopman model also exhibits the smallest standard deviation of the NRMSE across states. This implies that the Koopman model more consistently captures the real behavior of all six states of the system, rather than just a subset of them. Fig.3 illustrates the ability of the Koopman-based model to predict the position of the end effector over a 30 second time horizon.
V Discussion and Conclusion
We have successfully applied a system identification technique based on Koopman operator theory to a soft robot and shown that the model generated outperforms those constructed by several other state-of-the-art nonlinear system identification methods. Perhaps unsurprisingly, the linear state space model was unable to capture the nonlinear dynamics of the robot as well as the Koopman model. As for the nonlinear models, there are several likely reasons why the performance of the Koopman model was superior. Since the Koopman model is a state space model, simulations can be initialized from the same initial condition as the real system. This is not the case for the other learned nonlinear models which do not have an internal state corresponding to the physical state of the robot. Rather, they act as black-box models only capable of mapping inputs to outputs.
Another advantage of the Koopman model is that its quality does not depend on an initial model estimate or tuning parameters. By iterating over the set of all initializations and tuning parameters, one may be able to generate better performing models than those shown; unfortunately, this multivariate trial-and-error process may not affect results in a predictable way. In contrast, the only tuning parameter involved in the Koopman method is the degree of the monomial basis, which has a magnitude that is directly proportional to model accuracy.
While the results here are promising, there are practical challenges to extending the Koopman approach to higher dimensional systems. As the dimension of the state space increases, so does the size of the monomial basis of the finite-dimensional subspace of observables. This greatly increases the size of the matrix equations that must be solved, leading to computational intractability for sufficiently high dimensional systems. However, if some information about the system is known beforehand, this issue could be counteracted by choosing a more suitable basis for the observables. For example, if the system exhibits oscillatory motion, a lower dimensional fourier basis may be more suitable than a monomial basis to represent the behavior. Such an extension of the method is left to future work.
Soft robots are notoriously difficult to model, but amenable to large-scale data collection and data-driven modeling methods. This paper demonstrates the utility of Koopman operator theory to make accurate nonlinear dynamical models easier to construct, enabling the rapid development of control strategies that exploit the unique characteristics of soft robots. Future work will aim to generalize this approach to higher dimensional models, non-polynomial models, and models that account for external loading and contact forces.
-  D. Rus and M. T. Tolley, “Design, fabrication and control of soft robots,” Nature, vol. 521, no. 7553, p. 467, 2015.
-  C. Majidi, “Soft robotics: a perspective—current trends and prospects for the future,” Soft Robotics, vol. 1, no. 1, pp. 5–11, 2014.
-  H. Lipson, “Challenges and opportunities for design, simulation, and fabrication of soft robots,” Soft Robotics, vol. 1, no. 1, pp. 21–27, 2014.
-  W. Felt, S. Lu, and C. D. Remy, “Modeling and design of “smart braid” inductance sensors for fiber-reinforced elastomeric enclosures,” IEEE Sensors Journal, vol. 18, no. 7, pp. 2827–2835, 2018.
-  W. Felt, M. J. Telleria, T. F. Allen, G. Hein, J. B. Pompa, K. Albert, and C. D. Remy, “An inductance-based sensing system for bellows-driven continuum joints in soft robots,” Autonomous Robots, pp. 1–14, 2017.
-  S. Yang and N. Lu, “Gauge factor and stretchability of silicon-on-polymer strain gauges,” Sensors, vol. 13, no. 7, pp. 8577–8594, 2013.
-  D.-H. Kim, N. Lu, R. Ma, Y.-S. Kim, R.-H. Kim, S. Wang, J. Wu, S. M. Won, H. Tao, A. Islam, et al., “Epidermal electronics,” science, vol. 333, no. 6044, pp. 838–843, 2011.
-  C. Della Santina, M. Bianchi, G. Grioli, F. Angelini, M. Catalano, M. Garabini, and A. Bicchi, “Controlling soft robots: balancing feedback and feedforward elements,” IEEE Robotics & Automation Magazine, vol. 24, no. 3, pp. 75–83, 2017.
-  J. Bishop-Moser and S. Kota, “Design and modeling of generalized fiber-reinforced pneumatic soft actuators,” IEEE Transactions on Robotics, vol. 31, no. 3, pp. 536–545, 2015.
-  D. Bruder, A. Sedal, R. Vasudevan, and C. D. Remy, “Force generation by parallel combinations of fiber-reinforced fluid-driven actuators,” IEEE Robotics and Automation Letters, vol. 3, no. 4, pp. 3999–4006, Oct 2018.
-  F. Renda, M. Giorelli, M. Calisti, M. Cianchetti, and C. Laschi, “Dynamic model of a multibending soft robot arm driven by cables,” IEEE Transactions on Robotics, vol. 30, no. 5, pp. 1109–1122, 2014.
-  W. Felt and C. D. Remy, “A closed-form kinematic model for fiber-reinforced elastomeric enclosures,” Journal of Mechanisms and Robotics, vol. 10, no. 1, p. 014501, 2018.
-  S. Neppalli, M. A. Csencsits, B. A. Jones, and I. D. Walker, “Closed-form inverse kinematics for continuum manipulators,” Advanced Robotics, vol. 23, no. 15, pp. 2077–2091, 2009.
-  R. J. Webster III and B. A. Jones, “Design and kinematic modeling of constant curvature continuum robots: A review,” The International Journal of Robotics Research, vol. 29, no. 13, pp. 1661–1683, 2010.
-  B. A. Jones and I. D. Walker, “Kinematics for multisection continuum robots,” IEEE Transactions on Robotics, vol. 22, no. 1, pp. 43–55, 2006.
-  T. George Thuruthel, Y. Ansari, E. Falotico, and C. Laschi, “Control strategies for soft robotic manipulators: A survey,” Soft robotics, vol. 5, no. 2, pp. 149–163, 2018.
-  I. A. Gravagne, C. D. Rahn, and I. D. Walker, “Large deflection dynamics and control for planar continuum robots,” IEEE/ASME transactions on mechatronics, vol. 8, no. 2, pp. 299–307, 2003.
-  D. Trivedi, A. Lotfi, and C. D. Rahn, “Geometrically exact models for soft robotic manipulators,” IEEE Transactions on Robotics, vol. 24, no. 4, pp. 773–780, 2008.
-  D. Bruder, A. Sedal, J. Bishop-Moser, S. Kota, and R. Vasudevan, “Model based control of fiber reinforced elastofluidic enclosures,” in Robotics and Automation (ICRA), 2017 IEEE International Conference on. IEEE, 2017, pp. 5539–5544.
-  A. Sedal, D. Bruder, J. Bishop-Moser, R. Vasudevan, and S. Kota, “A constitutive model for torsional loads on fluid-driven soft robots,” in ASME 2017 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference. American Society of Mechanical Engineers, 2017, pp. V05AT08A016–V05AT08A016.
-  J. Bishop-Moser, G. Krishnan, C. Kim, and S. Kota, “Design of soft robotic actuators using fluid-filled fiber-reinforced elastomeric enclosures in parallel combinations,” in Intelligent Robots and Systems (IROS), 2012 IEEE/RSJ International Conference on. IEEE, 2012, pp. 4264–4269.
-  S. Boyd and L. Vandenberghe, Convex optimization. Cambridge university press, 2004.
-  M. T. Gillespie, C. M. Best, E. C. Townsend, D. Wingate, and M. D. Killpack, “Learning nonlinear dynamic models of soft robots for model predictive control with neural networks,” in 2018 IEEE International Conference on Soft Robotics (RoboSoft). IEEE, 2018.
-  L. Ljung, System identification: theory for the user. Prentice-hall, 1987.
-  M. Budišić, R. Mohr, and I. Mezić, “Applied koopmanism,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 22, no. 4, p. 047510, 2012.
-  M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, “A data–driven approximation of the koopman operator: Extending dynamic mode decomposition,” Journal of Nonlinear Science, vol. 25, no. 6, pp. 1307–1346, 2015.
-  A. Mauroy and J. Goncalves, “Linear identification of nonlinear systems: A lifting technique based on the koopman operator,” arXiv preprint arXiv:1605.04457, 2016.
-  ——, “Koopman-based lifting techniques for nonlinear systems identification,” arXiv preprint arXiv:1709.02003, 2017.
-  I. Abraham, G. de la Torre, and T. Murphey, “Model-based control using koopman operators,” in Proceedings of Robotics: Science and Systems, Cambridge, Massachusetts, July 2017.
-  A. Lasota and M. C. Mackey, Chaos, fractals, and noise: stochastic aspects of dynamics. Springer Science & Business Media, 2013, vol. 97.
-  N. J. Higham, Functions of matrices: theory and computation. Siam, 2008, vol. 104.
-  MATLAB, version 7.10.0 (R2017a). Natick, Massachusetts: The MathWorks Inc., 2017.