Learning Discrepancy Models From Experimental Data

Lear ning Discr epancy Models Fr om Experimental Data Kadierdan Kaheman 1 , Eurika Kaiser 1 , Benjamin Strom 1 , J. Nathan Kutz 2 and Ste ven L. Brunton 1 Abstract — First principles modeling of physical systems has led to significant technological advances across all branches of science. For nonlinear systems, ho wever , small modeling errors can lead to significant de viations from the true, measured behavior . Even in mechanical systems, where the equations are assumed to be well-known, ther e are often model discr epan- cies corr esponding to nonlinear friction, wind resistance, etc. Discovering models for these discrepancies r emains an open challenge f or many complex systems. In this work, we use the sparse identification of nonlinear dynamics (SINDy) algorithm to discover a model f or the discrepancy between a simplified model and measurement data. In particular , we assume that the model mismatch can be sparsely represented in a library of candidate model terms. W e demonstrate the efficacy of our approach on several examples including experimental data from a double pendulum on a cart. W e further design and implement a feed-forward controller in simulations, showing improv ement with a discrepancy model. I . I N T R O D U C T I O N There are currently tw o dominant approaches in dynamical systems modeling: 1) deriving equations of motion from gov erning physical principles, and 2) data-dri ven system identification. These two approaches are typically used in isolation to dev elop an end-to-end model. Data-driv en mod- els, based on emerging techniques in machine learning, are promising for the description of complex systems [1]. Ho w- ev er , these models are often opaque, and it is difficult to cap- ture the correct dynamical system structure to satisfy basic constraints and conservation laws. In contrast, first-principles physics models capture constraints and conservation laws by design [2], b ut for complex systems, these models are either ov erly simplistic or exceedingly expensiv e to resolve all multiscale interactions. In this work, we propose a hybrid modeling approach, where we employ data-dri ven techniques to model the discrepancy between a simplified or insufficient physical model and observ ed measurements. W e emplo y the sparse identification of nonlinear dynamics (SINDy) framew ork [3] to disco ver parsimonious and interpretable discrepancy models. Thus, we lev erage prior knowledge of the simplified physics, while more accurately modeling details of the true system. As an illustrativ e example for this paper , we consider a double pendulum on a cart experiment. W e typically model this type of experiment as a simple mechanical system with a few degrees of freedom, described by either a Hamiltonian 1 Kadierdan Kaheman, Eurika Kaiser, Benjamin Strom, and Steven L. Brunton are with the Department of Mechanical Engineering, Univ ersity of W ashington, Seattle, W A, 98105, USA. 2 J. Nathan Kutz is with the Department of Applied Mathematics, Univ ersity of W ashington, Seattle, W A, 98105, USA. E-mails: {kadierk, eurika, strombw , kutz, sbrunton}@uw .edu or Lagrangian [2]. Ho wev er , this model neglects real-world effects such as nonlinear friction, bearing chatter , and wind resistance. These factors may all reduce control performance, for example, in model predictiv e control (MPC) [4], [5], [6], where prediction errors adversely affect rob ustness [7]. For this specific system, designing a feed-forward control law to swing up the pendulum from rest requires a highly accurate model. Even a small mismatch in the system model may result in a considerable deviation in the computed trajectory , since it is a chaotic system. It is also challenging to obtain a data-driv en model of this system with the correct structure, so it is beneficial to incorporate prior physical knowledge in the form of a simplified Hamiltonian or Lagrangian. A. Problem statement: Discrepancy modeling There are sev eral reasons why model discrepancies oc- cur [8], [9]. First, there may be measurement noise and exogenous disturbances. In this case, the Kalman filter may be thought of as a discrepancy model where the mismatch between a simplified model and observations is assumed to be a Gaussian process [10]. Next, the parameters of the system may be inaccurately modeled. Even worse, the structure of the model may not be correct, either because important terms are missing or erroneous terms are present. This is known as model inadequacy or model structure mis- match. Other challenges include incomplete measurements and latent v ariables, delays, and sensitive dependence on initial data in chaotic systems. In this work, we focus on parameter and structural un- certainties, although a broader frame work is the subject of ongoing work. W e consider dynamical systems of the form d dt x = f ( x , u ; µ ) , (1) where x ∈ R n is the state, u ∈ R r is the control input, µ ∈ R p are the parameters, and f are the dynamics. W e assume full-state measurements, although generally the state must be estimated from limited measurements. The discrepancy modeling problem seeks to model the difference between a quantity of interest φ ( t ) from a physical model φ m ( t ) and the observed v alue φ o ( t ) : δ φ ( t ) = φ o ( t ) − φ m ( t ) , (2) where δ φ is the discrepancy . Here, we consider the quantity of interest φ ( t ) to be the dynamics themselv es, or more precisely the rate of change of the state in time. Howe ver , this framework is more general and can incorporate sev eral other forms of model discrepancy . For example, φ ( t ) may be a conserved quantity , such as the Hamiltonian, from which the dynamics are deriv ed. B. Related work Accurate system modeling is a critical task in science and engineering, including for autonomous robotics [11], [12] and process control [13], [14], [15]. A model that is able to accurately predict the future state of a system is imperativ e for prediction and control. Ho we ver , the high-dimensional, nonlinear , and multi-scale nature of many systems renders modeling a challenging task. System identification has reached a high degree of ma- turity encompassing myriad techniques to identify linear and nonlinear systems [16], [17] from data, including state- space modeling via the eigensystem realization algorithm (ERA) [18] and other subspace identification methods, V olterra series [19], [20], linear and nonlinear autoregressiv e models [21] (e.g., ARX, ARMA, NARX, and N ARMAX), and neural network models [22], [23], to name only a fe w . Machine learning techniques such as manifold learning and non-parametric modeling have also been useful for identify- ing nonlinear systems [24], [25]. There is an increasing shift from black-box modeling to dev eloping models that are physically intuiti ve and inter- pretable, as well as models that are constrained models with known prior information. For instance, genetic programming can be used to infer go verning equations from data [26], [27]. The recent SINDy approach [3] identifies parsimonious and interpretable models by promoting sparsity , and it has been extended to incorporate the effect of control [6] and to take into account expert kno wledge, such as symmetries and conservation laws [28]. Ho we ver , it may be challenging to identify certain types of systems, such as those containing rational function nonlinearities [29]. Instead of modeling the system entirely from data, when partial information is av ail- able, such as an idealized Hamiltonian, data-dri ven modeling procedures may focus only on modeling the discrepancy . One way to compensate for model discrepanc y is with Bayesian approaches [8], where a model discrepancy func- tion is learned and the model output uncertainty is quantified from data [30]. Howe ver , this method requires the selec- tion of a prior form for the model discrepancy function, which is dif ficult due to the lack of physical kno wledge about the model structure [31]. Furthermore, this method may potentially introduce bias in the model parameters if calibration is performed without kno wledge of the model discrepancy [32], [33]. Alternativ ely , reinforcement learning can be used to learn the mismatch between a model and the actual dynamics, for example in robotics [34]. The learned dynamics are then used with the conceptual model to control the real system. Howe ver , neither Bayesian methods nor reinforcement learning can provide an interpretable represen- tation for the model mismatch, and thus conceal the physical meaning of the discrepancy model discov ered. C. Contributions of this work In the present work, we leverage SINDy to compensate for model parameter and structure mismatch giv en an im- perfect model of the system. This serves se veral purposes: (1) Prior partial knowledge on the system or a previously learned model may be a v ailable and can be incorporated to aid the modeling process and improve prediction accu- racy . (2) The SINDy algorithm suf fers from the curse of dimensionality due to the growth of the library size with increasing number of v ariables in the system, which makes it challenging to discover the full governing equations when a large library size is required. Howe ver , learning the model mismatch significantly reduces this burden, focusing SINDy only on modeling the mismatch. (3) Learning interpretable representations of the model mismatch may inform physical intuition and generalize beyond the training data. W e will demonstrate that SINDy is able to model discrepancies such as incorrect system parameters and model inadequacy (or structure) errors. The learned SINDy discrepancy model can then be used to enhance the imperfect model, providing an improv ed description of the system dynamics. I I . B AC K G R OU N D : S PA R S E I D E N T I FI C A T I O N O F N O N L I N E A R D Y N A M I CS Here, we briefly introduce the SINDy algorithm for iden- tifying dynamics from data [3]. W e consider a system of coupled nonlinear ordinary differential equations (ODE): d dt x ( t ) = f ( x ( t ); µ ) , (3) with time-dependent state x ( t ) = [ x 1 ( t ) , . . ., x n ( t )] T ∈ R n . The function f ( x ( t ); µ ) describes the underlying dynamics of the system; it is possible to extend this formulation to include control [6], although we omit it here for simplicity . f is often sparse in the space of candidate functions. If this assumption holds, sparse regression [3], [35] can be used to infer f from measurement data as explained in the following. Step 1: F orm a data matrix X ∈ R m × n and matrix of deriv atives ˙ X ∈ R m × n sampled at times t 1 , . . . , t m : X = [ x ( t 1 ) x ( t 2 ) · · · x ( t m )] T , (4a) ˙ X = [ ˙ x ( t 1 ) ˙ x ( t 2 ) · · · ˙ x ( t m )] T . (4b) Step 2: Construct a library of candidate functions Θ ( X ) Θ ( X ) = [ θ 1 ( X ) θ 2 ( X ) θ 3 ( X ) · · · θ p ( X )] , (5) where θ i ( x ) is a candidate function that may be present in the dynamics f . For e xample, trigonomet- ric functions θ i ( x ) = sin( x j ) or polynomial func- tions θ i ( x ) = x j x 2 k may be considered. The data matrix Θ is obtained by ev aluating these functions at all m entries of the data matrix X . Step 3: F orm the sparse regression problem such that ˙ X = Θ ( X ) Ξ , (6) where the matrix Ξ = [ ξ 1 , . . . , ξ n ] ∈ R p × n com- prises sparse vectors ξ i ∈ R p × 1 that select the active terms in the library Θ . These are found via sparse optimization [35] to minimize the following: ξ i = argmin ξ 0 i    Θ ( X ) ξ i − ˙ X    2 + λ k ξ i k 1 . (7) I I I . D A TA - D R I V E N D I S C OV E RY O F M O D E L D I S C R E PA N C Y In this section, we use SINDy to model discrepancies. A. Discrepancy modeling for systems without control Although it is possible to model any physical quantity , we consider modeling the dynamics themselves [3]. Consider noisy measurements from a true dynamical system f φ o ( t ) = f ( x ( t ); µ ) + , (8) where  ∈ R n is the measurement noise. The model output for this system is represented as φ m ( t ) = f m ( x ( t ); µ 1 ) . (9) Parameter error, µ 6 = µ 1 , model inadequacy , f 6 = f m , and measurement error will cause a model mismatch, given by δ φ ( t ) = φ o ( t ) − φ m ( t ) = f ( x ( t ); µ ) − f m ( x ( t ); µ 1 ) . (10) W e then construct a model for the discrepancy δ φ ( t ) as a function of the state and new parameters µ 2 , capturing the parameter and structure mismatch: δ φ ( t ) = g ( x ( t ); µ 2 ) . (11) Model error data is collected in δ Φ = [ δ φ ( t 1 ) δ φ ( t 2 ) · · · δ φ ( t m )] T , (12) where δ Φ ∈ R m × n . Giv en data δ Φ and X , we use SINDy to sparsely represent g ( x ( t ); µ 2 ) in a library of candidate functions Θ ( X ) . Con- structing this library typically requires some prior kno wledge of the system to select a suitable basis in which the discrep- ancy model will be sparse. The sparse regression problem is then formulated as δ Φ = Θ ( X ) Ξ . (13) Realistically , we will only have access to measurements of x ( t ) , from which we may derive d dt x ( t ) , and the quantity of interest becomes φ ( t ) ≈ d dt x ( t ) . The model output is φ m = f m ( x ; µ 1 ) , the discrepancy is δ φ = d dt x ( t ) − f m ( x ; µ 1 ) , and we hav e δ ˙ X := ˙ X − f m ( X ; µ 1 ) = Θ ( X ) Ξ . (14) W e solve for the sparsest coef ficient matrix Ξ that satisfies Eq. (14) using the SINDy approach. A schematic of the method is displayed in Fig. 1. B. Discrepancy modeling for systems with control W e now extend the formulation above to identify dynamics that are affected by a control input u ( t ) ∈ R q : φ o ( t ) = f ( x ( t ) , u ( t ); µ ) + , (15) The model output for this system is represented as φ m ( t ) = f m ( x ( t ) , u ( t ); µ 1 ) . (16) As with the uncontrolled case, due to the model inadequacy or parameter error , there will be a discrepancy δ φ ( t ) = φ o ( t ) − φ m ( t ) = f ( x ( t ) , u ( t ); µ ) − f m ( x ( t ) , u ( t ); µ 1 ) . (17) Our goal is to model δ φ ( t ) and we further assume it is a function of the system state x , control input u , and new parameters µ 2 : δ φ ( t ) = g ( x ( t ) , u ( t ); µ 2 ) . (18) T o achieve this goal, we form a sparse regression problem δ Φ = Θ ( X , U ) ξ , (19) where U ∈ R m × r is U = [ u ( t 1 ) u ( t 2 ) · · · u ( t m )] T . (20) As before, we will actually be observing φ ( t ) ≈ d dt x ( t ) with model φ m = f m ( x , u ; µ 1 ) . The discrepancy becomes δ φ = d dt x ( t ) − f m ( x , u ; µ 1 ) , resulting in the following regression problem: δ ˙ X := ˙ X − f m ( X , U ; µ 1 ) = Θ ( X , U ) Ξ . (21) Note that the library Θ ( X , U ) has the form Θ ( X , U ) = [ θ 1 ( X , U ) θ 2 ( X , U ) · · · θ v ( X , U )] , (22) where θ i ( X , U ) ∈ R m × 1 is a candidate function that may explain the discrepancy δ φ ( t ) . The functions θ i ( X , U ) can be any combination of X and U . For example, θ i ( X , U ) = sin( X ) cos( U ) , θ i ( X , U ) = X U 2 . By solving Eq. (19) we are able to model g ( x ( t ) , u ( t ); δ µ 2 ) . I V . E X A M P L E S In this section, we demonstrate applications of the pro- posed approach to discover model discrepancies. W e start with an illustrati ve example and then apply this method to experimental data from a double pendulum on a cart. A. V an der P ol oscillator T o begin, we will focus on an illustrati ve example, the V an der Pol oscillator , and show how the SINDy method can be used to compensate for both parameter errors and model inadequacy . The V an der Pol oscillator is given by ˙ x 1 = x 2 , (23a) ˙ x 2 = α  1 − x 2 1  x 2 − x 1 , (23b) with parameter α = 0 . 5 . 1) P arameter mismatc h: Suppose we do not know the true parameter α , but instead have an approximation α 1 = 0 . 1 . It is possible to compensate for this model discrepancy caused by parameter mismatch. W e first gather the measurement data of the actual system, in this case by integrating the true dynamics using a fourth order Runge Kutta scheme. W e integrate for 25 time units with time step ∆ t = 0 . 01 and initial condition is x 0 = (0 . 5 , 0) . Random Gaussian measurement noise with amplitude 0 . 01 is added to the data. W e ev aluate the model dynamics f ( x ( t ); α 1 ) , with the measured trajectory x ( t ) and the inaccurate parameter α 1 . The discrepancy between the true system output and the model is then given by δ ˙ x ( t ) := d dt x ( t ) − f ( x ( t ); α 1 ) . The augmented error matrix is formed as δ ˙ X =  δ ˙ x ( t 0 ) δ ˙ x ( t 1 ) δ ˙ x ( t 2 ) . . . δ ˙ x ( t m )  T , (24) Fig. 1: Illustration of SINDy to discover the system-model mismatch for the V an der Pol oscillator . The system is simulated to generate the measurement data. Measurement data x 1 and x 2 are provided to calculate the output estimated by our imperfect model. Sparse regression is used to infer the discrepancy model for the difference between the actual output and estimated output. The discrepancy model is then combined with the imperfect model to provide a better estimation of system dynamics. The model is then cross-validated on a new initial condition x 0 = ( − 0 . 2 , − 0 . 3) . As we can see, the model discovered by SINDy is able to compensate for the discrepancy between the actual system and flawed model. and the augmented state is X =  x ( t 0 ) x ( t 1 ) x ( t 2 ) . . . x ( t m )  T . (25) The next step is to ev aluate the library on the data. Here, we construct the library as Θ ( x 1 , x 2 ) = Time      y Functions − − − − − − − − − − − − − − − − − − − − − − − − − − →   | | | | | | | 1 x 1 x 2 x 1 x 2 x 2 1 x 2 x 1 x 2 2 x 2 1 x 2 2 | | | | | | |   , (26) with polynomial terms up to second order . Finally , we form the sparse regression problem δ ˙ X = Θ ( x 1 , x 2 ) Ξ , (27) and solve for the coefficient matrix Ξ . Results are shown in the top row of Fig. 1. SINDy successfully identifies the model discrepancy . 2) Model Inadequacy (or structur e mismatch): W e now assume that the model discrepancy is caused by model inadequacy , where the model f m is missing the linear term in the second equation of the V an der Pol system: ˙ x 1 = x 2 , (28a) ˙ x 2 = α  1 − x 2 1  x 2 . (28b) W e assume that the parameter α = 0 . 5 is correct. The data, x , ˙ x , δ ˙ x is collected and the library of functions is constructed as before. The sparse coef ficient matrix Ξ is determined via SINDy , and the model discrepancy is suc- cessfully modeled as shown in the bottom ro w of Fig. 1. T o summarize, the SINDy method successfully identifies the discrepanc y between the underlying go verning equation and an inaccurate model. The imperfect model is then combined with the SINDy discrepancy model, resulting in an accurate prediction of the true system dynamics. B. Experimental double pendulum on a cart W e now demonstrate the use of SINDy to identify the time-dependent Hamiltonian function for the double pendu- lum on a cart from experimental data. A non-dissipati ve pendulum is a quintessential example of a conservati ve system, which admits in variants such as the Hamiltonian function, from which the gov erning equations can be deri ved. The conservati v e Hamiltonian is given by H c ( q , p ) = T + V , where q and p are the generalized position and momentum of the system, and T and V represent the kinetic and potential energy of the system; the Hamiltonian H c is the total ener gy , which is constant along a trajectory . Since we measure q and compute the time deriv ati ve ˙ q , we will represent H c as a function of q and ˙ q in the following. Real-world mechanical systems are, ho wev er , generally not conservati ve, but instead exhibit friction and damping from joints and wind resistance. Thus, the total energy is not conserved, and instead decays o ver time without additional exogenous ener gy input. While the potential and kinetic Fig. 2: Double pendulum on a cart system with schematic (bottom left) and zoom of the two-link pendulum (bottom right). The pendulum cart is locked at its position to pre vent horizontal movement. energy terms can be easily formulated, dissipation terms can be more challenging to deri ve. In this work, we seek to iden- tify the time-dependent dissipativ e effects, by modeling the difference between the conservati ve model Hamiltonian and the measured energy of the system. In practice, the observed energy is obtained by ev aluating the idealized conservati v e Hamiltonian, consisting of the kinetic and potential terms T and V , on the measured trajectory . The model energy is giv en by ev aluating the idealized Hamiltonian on the initial condition. Thus, the difference gi ves the ener gy dissipation, which we will model with SINDy: δ H ( q ( t ) , ˙ q ( t )) := H m ( q ( t ) , ˙ q ( t )) − H m ( q (0) , ˙ q (0)) . (29) 1) Problem F ormulation: W e consider the double pendu- lum on a cart as sho wn in Fig. 2. The kinetic and potential energy of the double pendulum (assuming a locked cart position) are T = 1 2 ( m 1 ( ˙ x 2 1 + ˙ y 2 1 ) + m 2 ( ˙ x 2 2 + ˙ y 2 2 )) + 1 2 ( I 1 ˙ ϕ 2 1 + I 2 ˙ ϕ 2 2 ) , (30a) V = ( m 1 y 1 + m 2 y 2 ) g , (30b) where I 1 and I 2 are the inertia, m 1 and m 2 are the masses, and l 1 and l 2 are the lengths of each pendulum arm, respectiv ely . ϕ 1 and ϕ 2 are the angles and ˙ ϕ 1 and ˙ ϕ 2 are the angular velocity of the first and second pendulum arm. The relativ e lengths to the center of mass of the first and second pendulum arm are giv en by a 1 and a 2 . The position of each pendulum arm’ s center of mass ( x 1 , y 1 ) and ( x 2 , y 2 ) in (30) can be determined as x 1 = a 1 sin( ϕ 1 ) , (31a) x 2 = l 1 sin( ϕ 1 ) + a 2 sin( ϕ 2 ) , (31b) y 1 = a 1 cos( ϕ 1 ) , (31c) y 2 = l 1 cos( ϕ 1 ) + a 2 cos( ϕ 2 ) . (31d) W e assume that we can accurately determine the kinetic and potential energy of the system H m ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) = T + V , (32) which represents the insufficient model for the total ener gy . The discrepanc y model δ H ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) then comprises the dissipativ e energy terms. The frictional torque of the pendulum arm can be modeled as Γ 1 = k 1 ˙ ϕ 1 and Γ 2 = k 2 ( ˙ ϕ 1 − ˙ ϕ 2 ) , where k 1 and k 2 are damping coef ficients. Then, δ H ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) is giv en by δ H ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) = Z t 0 Γ 1 ˙ ϕ 1 + Γ 2 ( ˙ ϕ 1 − ˙ ϕ 2 ) dt = Z t 0 k 1 ˙ ϕ 2 1 + k 2 ( ˙ ϕ 1 − ˙ ϕ 2 ) 2 dt. (33) In many engineering applications, the direct measurement of the frictional term is difficult if not impossible. Thus, we would like to use our data-driv en approach to learn a model for the dissipated energy (33). Suppose that all the states ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 can be measured or estimated, then H m ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) can be immediately calculated. Also, we can use the energy at the initial measurement time as reference for the total energy of the system, H ( ϕ 1 ( t 0 ) , ϕ 2 ( t 0 ) , ˙ ϕ 1 ( t 0 ) , ˙ ϕ 2 ( t 0 )) = E 0 . This allows us to calculate the dissipated ener gy as δ H ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) = E 0 − H m ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) . (34) T o model this energy discrepancy , we define a library Θ ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) that contains polynomial and Fourier terms up to the third order . Then, the discrepancy can be represented by δ H ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) = Θ ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) ξ , (35) where ξ is a sparse vector that contains the coefficients of the activ e terms in the library . 2) Results: Here, we present results for modeling the dissipation δ H ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) . The parameters of our ex- perimental system are displayed in T able. I. All parameters are obtained using a parameter estimation technique [36]. The pendulum mass and length are constrained during the parameter estimation so that there won’t be a large deviation from the measured value. The experiment w as initialized at a random position, and the angles ϕ 1 and ϕ 2 − ϕ 1 were collected with a sampling rate of ∆ t = 0 . 001 for a duration of 20 s. W e used a US Digital HUBDISK-1 1” transmissiv e rotary encoder disk, which pro vides 5000 counts per revolution. Then the angular velocities ˙ ϕ 1 and ˙ ϕ 2 are T ABLE I: The parameters of the experimental double pendulum system. Pendulum Mass ( kg ) Center of Mass ( m ) Inertia ( kg · m 2 ) Length ( m ) Damping Coefficient Gravity Constant ( m/s 2 ) 1st Arm 0 . 2704 0 . 1910 0 . 003 0 . 2667 7 . 24 × 10 − 4 9 . 818 2nd Arm 0 . 2056 0 . 1621 0 . 0011 0 . 2667 1 . 65 × 10 − 4 0 10 20 -1.16 -1.14 -1.12 -1.1 0 10 20 0 0.02 0.04 0.06 (a) Left: Initial energy of the system E 0 and the sum of kinetic and potential energy H m ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) , which decreases due to friction. Right: The difference between E 0 and H m ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) representing the ener gy dissipated by friction. 0 5 10 15 20 0 0.02 0.04 0.06 (b) The magnitude of ener gy dissipated by the friction δH ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) and the SINDy identified dynamics. Note that the discrepancy is positiv e, since it is defined as the ideal conserved energy minus the measured dissipativ e ener gy . 0 5 10 15 20 -1 0 1 10 -4 (c) The SINDy estimation error on the training data. The percentage is calculated using δH ( ϕ 1 ,ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) − δH SINDy max( | δH ( ϕ 1 ,ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) | ) × 100% . Fig. 3: Results for the identified model discrepancy caused by friction, demonstrated on training data. 0 5 10 15 20 0 0.02 0.04 0.06 (a) The magnitude of energy dissipated by the friction δH ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) and the SINDy estimation e valuated on validation data. 0 5 10 15 20 0 2 4 6 (b) The SINDy estimation error on v alidation data. The percentage is calculated using δH ( ϕ 1 ,ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) − δH SINDy max( | δH ( ϕ 1 ,ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) | ) × 100% . Fig. 4: Results for the identified model discrepancy caused by friction, demonstrated on validation data. approximated by taking numerical deri v ati ves on the ra w data. T o mitigate noise, we first smooth the raw ϕ 1 and ϕ 2 data using the Savitzk y-Golay filter and afterwards compute the deriv ative by numerical differentiation. Results of the discrepancy model, identified with SINDy on the time series of δ H ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) , are shown in Fig. 3 for the training data and in Fig. 4 for validation data. Note that in Fig. 3a the signal H m does not decrease asymptoticly . W e observe that the table on which the pendulum in mounted oscillates with the pendulum, storing a small amount of potential energy . The error of the SINDy prediction and the measured δ H ( ϕ 1 , ϕ 2 , ˙ ϕ 1 , ˙ ϕ 2 ) is shown in Fig. 3b. From Fig. 3c we see that SINDy accurately describes the effect of friction, explaining the model discrepanc y . Crossvalidation is critical to av oid overfitting. Hence, we test the identified model on a new validation dataset unseen during the training stage. The performance of the SINDy- discov ered model on the validation data is shown in Fig. 4a. Although the errors are larger on the validation dataset, the SINDy model is quite accurate in modeling the missing friction term, demonstrating the ability of the discrepancy model to generalize to new test cases. C. Double pendulum on a cart with control In the previous section, SINDy is utilized to discover discrepancy models in the total energy caused by dissipation when giv en an insufficient conservati v e energy model. W e now consider the case of an actuated double pendulum on a cart, where a control input is added to the system. This is an important generalization, since many real-life system are affected by control, which requires additional modeling. T o discov er the model discrepancy in systems with control, we follow the steps outlined in Sec. III-B. 1) Problem formulation: In this example, we consider data from a numerical simulation of the double pendulum on a cart, shown in Fig. 2, with parameters provided in T ab. I. The equations of motion can be represented as ˙ x = f ( x ) + [0 , 0 , 0 , 0 , 0 , 1] T u, (36) where x = [ ϕ 1 , ϕ 2 , s, ˙ ϕ 1 , ˙ ϕ 2 , ˙ s ] T is the state vector and s represents the displacement of the pendulum cart. W e choose the acceleration of the pendulum cart as control input so that u = ¨ s . W e refer to [36] for a deriv ation of the equations of motion and a technique to estimate parameters from experimental data. W e seek to perform the swing-up control of the double pendulum. Howe v er , the employed model is fla wed, as it contains an additional term and an incorrect parameter: f m ( x ) = f ( x ) + [sin( ϕ 1 ) , 0 , 0 , 0 , 0 , 0] T + [0 , 0 , 0 , 0 , 0 , 0 . 95] T u. (37) Fig. 5: Swing-up of the double pendulum on a cart using feed-forward control. T op: The feed-forward control input is calculated using the imperfect model, which is in then applied to the actual system. The pendulum cannot reach the up-right position. Bottom: System’ s response under the control input calculated using the hybrid model, which is comprised of the imperfect model and the identified discrepancy model. Since the discrepancy is identified, the ne w model mimics the true dynamics of the double pendulum on a cart system. The discrepancy model is then δ ˙ x = ˙ x − f m ( x ) = [ − sin( ϕ 1 ) , 0 , 0 , 0 , 0 , 0 . 05 u ] T , (38) which we seek to model. First, we generate data for the identification using the imperfect model. Specifically , we design a feed-forward control input to swing up the double pendulum based on the imperfect model f m ( x ) . When this control signal is applied to the actual system, the observed behavior de viates from the pre-planned trajectory due to the model discrepancy . This difference is then used to form a sparse regression problem, as in Eq. (22), to identify the discrepancy δ ˙ x . Finally , a ne w swing-up trajectory is designed using the hybrid model consisting of the imperfect model, augmented with the discrepancy model. 2) Results: Simulation results for the swing-up of the double pendulum are sho wn in Fig. 5. The feed-forward trajectory is determined as the solution to an optimization problem [37]. The simulation time step is ∆ t = 0 . 001 , the prediction horizon is 2 , and the total swing-up time is chosen to be T = 2 . 7 . The weight matrices are Q = diag [10 , 10 , 20 , 1 , 1 , 0 . 1] and R = 1 for the state and input, respectiv ely . It is demonstrated in Fig. 5 that SINDy correctly identifies the model mismatch shown in Eq. (38) and that the new model mimics the actual dynamics of the system, which results in a successful swing-up of the pendulum. V . C O N C L U S I O N In this work, we present a data-driven frame work to model discrepancies between observ ations and simplified or incor - rect physical models. In particular , we lev erage the sparse identification of nonlinear dynamics (SINDy) [3] algorithm to discover model discrepancies caused by parameter mis- match or model inadequac y . W e demonstrate this approach on several systems, including the V an der Pol oscillator and experimental and numerical data from a double pendulum on an actuated cart. Our results suggest that a hybrid discrepancy modeling approach, in volving a physics-based model and a data- driv en correction, may hav e sev eral benefits. First, we in- corporate prior kno wledge and enforce conserv ation laws and constraints that are notoriously challenging in data- driv en approaches. In addition, focusing SINDy on the model mismatch is a much simpler task than trying to model all system dynamics at once, which would in volv e a lar ge library of candidate terms and a potentially ill-conditioned inv erse problem. This ill-conditioned problem often leads to the mis- estimation of parameters, such as the mass and length of the pendulum arms, which are accurately measured ahead of time. Instead, we are able to focus the data-driv en effort on modeling the few terms and parameters that cause the discrepancy . There are sev eral future directions suggested by this work. First, it will be important to extend this framework to include noise and exogenous disturbances, as in the Kalman filter , along with partial measurements. It will also be interesting to use our discrepancy models for the experimental swing- up control of the double pendulum on a cart, which is the subject of ongoing work. A C K N O W L E D G M E N T SLB would like to acknowledge funding support from the Army Research Office (AR O W911NF-17-1-0306, W911NF- 17-1-0422, and W911NF-19-01-0045). EK gratefully ac- knowledges support by the “W ashington Research Foun- dation Fund for Innov ation in Data-Intensiv e Discovery" and a Data Science Environments project award from the Gordon and Betty Moore Foundation (A ward #2013-10- 29) and the Alfred P . Sloan F oundation (A ward #3835) to the Univ ersity of W ashington eScience Institute, and funding through the Mistletoe F oundation. JNK and SLB acknowledge support from the Defense Advanced Research Projects Agency (D ARP A- P A-18-01-FP-125). W e would like to acknowledge v aluable discussions with Brian DeSilva, T on y Piasko wy , and Aditya Nair . R E F E R E N C E S [1] S. L. Brunton and J. N. Kutz, Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Contr ol . Cambridge Univ ersity Press, 2019. [2] J. E. Marsden and T . S. Ratiu, Introduction to mechanics and symme- try , 2nd ed. Springer-V erlag, 1999. [3] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems, ” Pr oceedings of the National Academy of Sciences , vol. 113, no. 15, pp. 3932–3937, 2016. [4] E. F . Camacho and C. B. Alba, Model pr edictive control . Springer Science & Business Media, 2013. [5] U. Eren, A. Prach, B. B. K oçer , S. V . Rakovi ´ c, E. Kayacan, and B. Açıkme ¸ se, “Model predictive control in aerospace systems: Current state and opportunities, ” J. Guid. Control & Dyn. , 2017. [6] E. Kaiser, J. N. Kutz, and S. L. Brunton, “Sparse identification of nonlinear dynamics for model predictive control in the low-data limit, ” Pr oc. Roy. Soc. A , vol. 474, no. 2219, 2018. [7] V . Botelho, J. O. Trierweiler , M. F arenzena, and R. Duraiski, “Method- ology for detecting model–plant mismatches affecting model pre- dictiv e control performance, ” Industrial & Engineering Chemistry Resear ch , vol. 54, no. 48, pp. 12 072–12 085, Nov 2015. [8] M. C. Kennedy and A. O Hagan, “Bayesian calibration of computer models, ” Journal of the Royal Statistical Society: Series B (Statistical Methodology) , vol. 63, no. 3, pp. 425–464, Aug 2001. [9] P . Arendt, D. W . Apley , and W . Chen, “Quantification of model un- certainty: Calibration, model discrepancy , and identifiability , ” Journal of Mechanical Design , vol. 134, p. 100908, 09 2012. [10] R. E. Kalman, “ A new approach to linear filtering and prediction problems, ” J. Fluids Eng. , vol. 82, no. 1, pp. 35–45, 1960. [11] B. Kim, D. Necsulescu *, and J. Sasiadek, “ Autonomous mobile robot model predictive control, ” International Journal of Control , vol. 77, no. 16, pp. 1438–1445, Nov 2004. [12] A. Rezaee, “Model predictive controller for mobile robot, ” T rans. En vir on. & Elec. Eng. , vol. 2, no. 2, p. 18, 2017. [13] O. Aarna, “Dynamic chemical plant model and its application, ” IF AC Pr oceedings V olumes , vol. 11, no. 1, pp. 279 – 286, 1978, 7th T riennial W orld Congress of the IF A C on A Link Between Science and Applications of Automatic Control, Helsinki, Finland, 12-16 June. [14] B. Pyakillya and S. Kladiev , “Identification of chemical reactor plant’s mathematical model, ” MATEC W eb of Conferences , vol. 37, p. 01045, 2015. [15] K. Arulmaran and J. Liu, “Handling model plant mismatch in state estimation using a multiple-model-based approach, ” Indust. & Eng. Chem. Resear ch , vol. 56, no. 18, pp. 5339–5351, 2017. [16] O. Nelles, Nonlinear system identification: fr om classical approac hes to neural networks and fuzzy models . Springer , 2013. [17] L. Ljung, “Perspectiv es on system identification, ” Annual Reviews in Contr ol , vol. 34, no. 1, pp. 1–12, 2010. [18] J. N. Juang and R. S. Pappa, “ An eigensystem realization algorithm for modal parameter identification and model reduction, ” J. of Guidance, Contr ol, and Dynamics , vol. 8, no. 5, pp. 620–627, 1985. [19] R. W . Brockett, “V olterra series and geometric control theory , ” Auto- matica , v ol. 12, no. 2, pp. 167–176, 1976. [20] B. R. Maner, F . J. Doyle, B. A. Ogunnaike, and R. K. Pearson, “ A nonlinear model predicti ve control scheme using second order Volterra models, ” in A CC , vol. 3, 1994, pp. 3253–3257. [21] H. Akaike, “Fitting autoregressi ve models for prediction, ” Ann Inst Stat Math , vol. 21, no. 1, pp. 243–247, 1969. [22] R. Lippmann, “ An introduction to computing with neural nets, ” IEEE Assp magazine , vol. 4, no. 2, pp. 4–22, 1987. [23] A. Draeger, S. Engell, and H. Ranke, “Model predictive control using neural networks, ” IEEE Contr ol Systems Magazine , vol. 15, no. 5, pp. 61–66, 1995. [24] J. C. Principe, L. W ang, and M. A. Motter , “Local dynamic modeling with self-org anizing maps and applications to nonlinear system iden- tification and control, ” Pr oceedings of the IEEE , vol. 86, no. 11, pp. 2240–2258, 1998. [25] J. Ko, D. J. Klein, D. Fox, and D. Haehnel, “Gaussian processes and reinforcement learning for identification and control of an autonomous blimp, ” in IEEE ICRA . IEEE, 2007, pp. 742–747. [26] J. Bongard and H. Lipson, “ Automated rev erse engineering of nonlin- ear dynamical systems, ” Proc. Natl. Acad. Sciences , vol. 104, no. 24, pp. 9943–9948, 2007. [27] M. Schmidt and H. Lipson, “Distilling free-form natural laws from experimental data, ” Science , vol. 324, no. 5923, pp. 81–85, 2009. [28] J.-C. Loiseau and S. L. Brunton, “Constrained sparse Galerkin regres- sion, ” Journal of Fluid Mechanics , vol. 838, pp. 42–67, 2018. [29] N. M. Mangan, S. L. Brunton, J. L. Proctor , and J. N. Kutz, “Inferring biological networks by sparse identification of nonlinear dynamics, ” IEEE T ransactions on Molecular , Biological, and Multi-Scale Com- munications , v ol. 2, no. 1, pp. 52–63, 2016. [30] J. A. McMahan and R. C. Smith, “Parameter estimation for predictiv e simulation of oscillatory systems with model discrepancy , ” IF A C- P apersOnLine , vol. 49, no. 18, pp. 428 – 433, 2016, 10th IF AC Symposium on Nonlinear Control Systems NOLCOS 2016. [31] Y . Ling, J. Mullins, and S. Mahade van, “Selection of model dis- crepancy priors in Bayesian calibration, ” J ournal of Computational Physics , v ol. 276, pp. 665 – 680, 2014. [32] J. Brynjarsdóttir and A. O Hagan, “Learning about physical parame- ters: the importance of model discrepancy , ” Inverse Problems , vol. 30, no. 11, p. 114007, Oct 2014. [33] J. A. Burns, E. M. Cliff, and K. Farlo w , “Parameter estimation and model discrepancy in control systems with delays, ” IF AC Proceedings , vol. 47, no. 3, pp. 11 679 – 11 684, 2014. [34] I. K oryakovskiy , M. Kudruss, H. V allery , R. Bab uska, and W . Caarls, “Model-plant mismatch compensation using reinforcement learning, ” IEEE Robot. & Aut. Lett. , vol. 3, no. 3, pp. 2471–2477, Jul 2018. [35] P . Zheng, T . Askham, S. L. Brunton, J. N. Kutz, and A. Y . Aravkin, “Sparse relaxed regularized regression: SR3, ” IEEE Access , vol. 7, no. 1, pp. 1404–1423, 2019. [36] K. Graichen, M. Treuer , and M. Zeitz, “Swing-up of the double pendu- lum on a cart by feedforward and feedback control with experimental validation, ” Automatica , vol. 43, no. 1, pp. 63 – 71, 2007. [37] H. Deng and T . Ohtsuka, “ A highly parallelizable Newton-type method for nonlinear model predicti ve control, ” IF AC-P apersOnLine , vol. 51, no. 20, pp. 349 – 355, 2018, 6th IF AC Conference on Nonlinear Model Predictiv e Control NMPC 2018.