Non-linear System Identification from Partial Observations via Iterative Smoothing and Learning
Authors: Kunal Menda
Jean de Becdelièvre
Jayesh K Gupta
Ilan Kroo
Mykel J. Kochenderfer
Zachary Manchester
Venue:  (2020)
Keywords: System Identification
Dynamical Systems
Partial Observations
Non-linear Programming
Expectation Maximization
Neural Networks

Abstract: System identification is the process of building a mathematical model of an unknown system from measurements of its inputs and outputs. It is a key step for model-based control, estimator design, and output prediction. This work presents an algorithm for non-linear offline system identification from partial observations, i.e. situations in which the system's full-state is not directly observable. The algorithm presented, called SISL, iteratively infers the system's full state through non-linear optimization and then updates the model parameters. We test our algorithm on a simulated system of coupled Lorenz attractors, showing our algorithm's ability to identify high-dimensional systems that prove intractable for particle-based approaches. We also use SISL to identify the dynamics of an aerobatic helicopter. By augmenting the state with unobserved fluid states, we learn a model that predicts the acceleration of the helicopter better than state-of-the-art approaches.

Under review as a conference paper at ICLR 2020
NON-LINEAR SYSTEM IDENTIFICATION FROM PARTIAL OBSERVATIONS VIA ITERATIVE SMOOTHING AND LEARNING
Anonymous authors Paper under double-bling blind review
ABSTRACT
System identification the process of building a mathematical model of an un- known system from measurements of its inputs and outputs. It is key step for model-based control, estimator design, and output prediction. This work presents algorithm for non-linear offline system identification from partial observations, i.. situations in which the system's 's full-state not directly observable. The algo- rithm presented, called SISL., iteratively infers the system's 's full state through non- linear optimization and then updates the model parameters. We test our algorithm on simulated system of coupled Lorenz attractors, showing our algorithm' 's abil- ity to identify high-dimensional systems that prove intractable for particle-base approaches. We also use SISL to identify the dynamics of aerobatic helicopter. By augmenting the state with unobserved fluid states, we learn model that pre- dicts the acceleration of the helicopter better than state-of-the- approaches.
INTRODUCTION
The performance of controllers and state-estimators for non-linear systems depends heavily on the quality of the model of system dynamics (Hou & Wang 2013). System-identification tackles the problem of learning or calibrating dynamics models from data (Ljung 1999), which is often a time- history of observations of the system and control inputs. In this work, we address the problem of learning dynamics models of partially observed, high-dimensional non-linear systems. That is, we consider situations in which the system's 's state cannot be inferred from single observation, but instead requires inference over time-series of observations.
The problem of identifying systems from partial observations arises in many roboties domains (Pun- jani & Abbeel, 2015) Cory & Tedrake 2008 Ordonez et al., 2017). Though we often have direct measurements of robot's pose and velocity, in many cases we cannot directly observe relevant quantities such as the temperature of actuators or the state the environment around the robot. Con- sider learning a dynamics model for an aerobatic helicopter. AAbbell et al. (2010) attempted to map only the helicopter's pose and velocity its acceleration and they found their model be inaccurate when predicting aggressive maneuvers. They posited that the substantial airflow generated by the helicopter affected the dynamics. Since it is often impossible to directly measure the state of the airflow around a vehicle, identification must be performed in a partially observed setting.
irstem-identiffcation is mature field with rich history (Ljung) 1999 2010). Various techniques can be classified by whether they apply to linear or non-linear systems, with partially or fully ob- served states. Additionally, techniques are applied in an online or batch-offline setting. This work presents an approach to offline identification of non-linear and partially observed systems. When a system is fully observed, i.e. its full state observed but corrupted by noise, a set of techniques called equation-error methods are typically employed (Aström & Eykhoff 1971). In such cases, we can consider observations as independent, and minimize the error between the observed state- derivatives and those predicted by the model given the control input and observed states. In par- tially observed settings, merely knowing the current input is insufficient to accurately predict the observation. Several black-box approaches exist to predict observations from time-series of in- puts. Autoregressive approaches directly map a time-history of past inputs to observations (Billings, 2013). Recurrent neural networks (Bailer-lones al. 1998, Zimmermann & Neuneier 2000) and

---

Under review as a conference paper at ICLR 2020
cebsspac--idetificclle methods (Van Overschee & De Moor 1994) also be used to learn black- box dynamical systems from this data.
However, in many cases prior knowledge can be used to specify structured, parameterized models of the system (Gupta et al. 2019). Such models can be trained with less data and used with a wider ar- ray of control and state-estimatie esttmmtttitiion techniques than non-linear black-box models (Gupta al. 2019
Lutter et al. 2019b a). Techniques used to identify partially observed structured models are often based on Expeetation-Maxxmizaatooo (EM) (Dempster al., 1977 Schön al., 201 Kantas al.) 2015 Ghahramani & Roweis 1999). An alternating procedure is performed in which a smoothing step uses the current system dynamics estimate to infer the distribution over state tratttttaaectoiiss, and is followed by learning step that uses this distribution to update the system dynamics estimate.
In the non-linear or Gaussian case, it is typically not possible analytically characterize the distribution over trajectories, and thus methods based on Sequential Monte-Carle such as Parti- cle Smoothing (PS) (Schön al., [2011] Kantas al. 2015). or Extended Kalman Smoothing (EKS) (Ghahramani & Roweis| 1999] are employed in the E-step. Though considered state-of-the- art for this problem, both methods become intractable high-dimensional state spaces. PS suffers from the curse of dimensionality, requiring an intractably large number of particles if the state space is high- high-dimensional (Snyder et al., 2008 Kantas et al. 2013), and an M-step that can be quadratic in complexity with respect to the number of particles (Schön et al. 2011). EKS-based methods are fast during the E-step, but the M-step requires approximations to integrate out state uncertainty, such fitting lineariiiitis with Radial Basis Function approximators, and scales poorly with the dimension of the state-space (Ghahramani & Roweis 1999).
In this work, we present system-identiiication algorithm that is suited for high-dimensional, non- linear, and partially observed systems. By assuming that the systems are close to deterministic, as often the case in robotics, we approximate the distribution over unobserved states using only their maximum-like likellhood (ML) point-estimate. Our algorithm, called SISL -System-identification via Iterative Smoothing and Learning) performs the following two steps until convergence:
In the smoothing or E-step, we use non-linear programming to tractably find the ML point- estimate of the unobserved states.
In the learning or M-step, we use the estimate of unobserved states to improve the estimate system parameters.
The idea to use an ML point-estimate in lieu of the distribution over unobserved states in the EM procedure's E-step not new, and, in general, does not guarantee monotonic convergence to a lo- cal optimum (Celeux & Govaert 1992). However, such an approximation is equivalent to regular EM if the ML point-estimate the only instance of unobserved variables with non-negligible prob- ability (Celeux & Govaert) 1992 Neal & Hinton 1998). We apply this idea to the problem of system-identifiiation for nearly deterministic systems, in which MI. point-estimates can serve as surrogates for the true distribution over unobserved state-trajectories.
The primary contribution of this work an algorithm for identifying non-linear, partially observed we specify the assumptions underpinning our algorithm and discuss the computational methodology for using it. In Section 3 empirically demonstrate that it is able to identify the parameters of nigh-dimensiioal system of coupled Lorenz attractors, a problem that proves intractable for particle- methods. We also demonstrate algorithm the problem of identifying the dynamics of aerobatic helicopter, and compare against various approaches including the faate-of-the-ad approach (Punjani & Abbeell 2015).
systems that is able to scale high-dimensional problems. In Section 2
METHODOLOGY
2.1 FORMAL PROBLEM STATEMENT
In this work, assume that given batch of trajectories containing observations V1.7 RT of a dynamical system as evolves over a time horizon input sequence -1 1. We assume that this dynamical system has state R that evolves and
T, possibly forced by some known
2

---

Under review as a conference paper at ICLR 2020
generates observations according to the following equations,
1++1 = f(1,, ue, t) + Wi, We ~ (-) V/ ( ue, vt, Ut ~ p-(-)
where wi is referred to the process noise and assumed to be additive for notational simplicity, but this is not required assumption. Without loss of generality, we can drop the dependence on Ut, absorbing into the dependence t.
the observation noise. wy and i2 both
We further assume that we are provided a class of parameterized models fo(2, 1) and go () for 0 θ that approximate the dynamical system's evolution and observation processes. The goal of our algorithm is to find the parameters 0 that maximize the likelihood of the observations. That is, we seek find:
OM.. arg max p(ye: 0)
(2)
arg max
p(yer Ti.r, 0)p(x1r 0)drar
Assuming, the system is Markovian in we can factorize the distributions
p (0) p(x)) II Pul - fo( ()),
(3)
p FITT,
II Po - go(xx,
In order to tractably solve the maximization problem in Equation particle-based EM tech- niques typically approximate the integral as an expectation over a particle set (Schön et al. 2011 Kantas 2015). However, since particle-based methods can struggle with high-dimensional spaces (Snyder et al. 2008 Kantas et al. 2015), this work seeks point estimates for the likelihood- maximizing state sequence, which can be found using non-linear programming.
2.2 SURROGATE OBJECTIVE FUNCTION
Instead of solving the full maximum likelihood problem shown in Equation (2). solve surrogate problem:
SisL = arg max max p(yi:r #1rr, 0)p(r1: l0)
(4)
The maximized surrogate objective is equivalent to the original objective the distribution over F1.T given V1.T and 0 is well-approximated by its ML point-estimate, i.e. the distribution is Dirac delta function. We expect this assumption to hold if the trajectories are sufficiently long, the dynamics are close to deterministic, and all relevant dynamic modes are persistently excited. Making these assumptions, we now present an algorithm for solving Equation
2.3 SYSTEM-IDENTIIIATIION VIA ITERATIVE SMOOTHING AND LEARNING
By taking the logarithm of the Akelihood-objectivg Equation and using the factorization from Equation 3
we get the SISL objective J(1:1 O) as follows:
J(r.r 0) ((gg(()))
log p. go 1))
log Pm (#e+1 (2.,
(5)
Jointly maximizing this objective over Vi.r and yields ASISL.
Though this objective can be optimized as is by a non-linear optimizer, it is not necessarily efficient to do since 0 and 2ı.r are highly coupled, leading inefficient and potentially unstable updates. For this reason, we take an iterative approach akin to EM, which performs Block-corddinate ascent on the [
3

---

Under review as a conference paper at ICLR 2020
At iteration k, we first perform smoothing by holding constant, and finding the ML point- esttimt for 21.T as follows:
(k)
(k-1) arg max JZI.T + P.ll 21T
(6)
#1.r
Here Ps scales a soft trust-region regularizer similar to damping terms found in Levenberg- Marquardt methods (Levenberg) 1944] Marquardi I963). The smoothing problem can be solved in (T) by taking advantage of sparsity in the smoothing objective's Hessian matrix.
In the learning step, we hold F1.7 constant and find:
g(k) arg max J [x (<]17,00) 0) + pello 1(k-1) log
Again, Pe scales a soft trust-region regularizer, and specifying log allows us to regularize toward prior. The above optimization problem can be solved using any non-linear optimizer, such as a first or second-order gradient descent scheme. The SISL algorithm iterates between the smoothing and learning steps until convergence.
EXPERIMENTS
The objective of our experiments is to demonstrate that the SISL algorithm is capable of identi- fying high-dimensional non-linear systems partially observed settings. We do so simulation by identifying the parameters of system of partially observed coupled Lorenz attractors, well as by identifying the dynamies of real aerobatic helicopter. In the second experiment, we build on previous analysis of the dataset (Abbeel et al. 2010; Punjani & Abbeel (015) by attempting to characterize the interaction of the helicopter with the fluid around it, without having any direct observation of the fluid state.
3.1 IDENTIFICATION OF COUPLED LORENZ SYSTEMS
In this didactic experiment, we will show that:
SISL learns unbiased parameter estimates of systems that are close to deterministic, and, 2. SISL scales to dige-niiensional problems in which particle-smoothing is intractable.
To justify these claims, we system that is sufficiently non-linear and partially observable to make particle-based smoothing methods intractable. We choose system of coupled Lorenz attrac- tors for this purpose, owing their ability to exhibit chaotic behavior and their in non-linear atmospheric and fluid flow models (Bergé et al. 1984). Arbitrary increase in state dimensionality can be achieved by coupling multiple individual attractors. The state of a system with K coupled Lorenz attractors is x R follows:
{..., L1,k, 22,4, .... The dynamics of the system are as
δ1 = 0x(22, x1,k)
22 11,b(Pk T3,k) 22,4
(8)
ia, I1,4 BLTTRR
1 + Hx
where H is an R3K x3K matrix.
We nominally set the parameters (0), Pk. Bk) to the values (10, 28, 8/3). and randomly sample the entries of H from normal distribution to generate chaotic and coupled behavior between attractors, while avoiding self-coupling These parameters are estimated during identification. In order to make the system partially observed, the observation y E I3 22) is found from r as follows:
y = Cx + v, V 2~ N((,000
where C E (3K-2) is known matrix with full row-rank, and is the observation noise sampled from Gaussian with diagonal covariance I.I.. The entries of C are also randomly sampled from standard normal distribution. the following experiments, we simulate the system for T 128 timesteps at sample rate of At 0.04s, and integrate the system using 4hh Runge- Runge-Kutta method. Initial conditions for each trajectory are sampled such that 21,6 A(-6, ((-6.2.5 2.5²) 22,k (- 2.5²), 23,k (224, 2.52).
4

---

Under review as a conference paper at ICLR 2020
2 Attractors
6 Attractors
10 particies
Frequency
40 particles
100
150 panticlss
640 particles
Figure 1: Histograms of weighted RMSE of PS with varied number of particles.
3.1.1 UNBIASED ESTIMATION IN DETERMINISTIC SETTINGS
To test the conditions under which SISL learns
unbiased parameter estimates, we simulate single Lorenz system with H = O, and known
dr
||
P
0.001
0.01
10.011(0.012)
28.000(0.001)
2 2.667(0.000)
C' 2x2x3 We introduce and vary the process
0.010
0.01
10.017(0.012)
28.000(000011
2 (0.00)
noise w ~ (0, T), and vary the observa-
0.100
0.01
10.064(0.036)
27.990(0.0133
2.676 (0.004)
0.001
0.05
10.006(0.016)
27.908 (0.002)
2 2.666(00001)
tion noise coefficient ou, and then attempt to estimate the parameters (a, p. 8). Using initial guesses within 10% of the system's 's true param- eter values, we run SISL on single sampled trajectory. For each choice of ow and d., we repeat this process for 10 random seeds.
0.001
0.10
0.908().022)
(7.995(0.004)
2.665(0.001)
Table 1: Mean parameter estimates and standard errors for single Lorenz system simulated with various aw and o...
Table 1 shows the mean and standard errors of parameter estimates for various ow and v.. We highlight in red the mean estimates that are not within two standard errors their true value. We see that a and are estimated without bias for all scenarios. However, the estimate of /3 appears to become biased as the process noise is increased, but not as the observation noise increased. This supports the assumption that the objective used in SISL is sound when systems evolve close to deterministically, but can be biased if it is not.
3.1.2 INTRACTABILITY OF PARTICLE SMOOTHING
State-of-the-ar methods for parameter identification of partially observed systems rely on particle- smoothing (PS) to estimate distribution over states T1:7 given trajectory of observations y1.T. We will experimentally demonstrate that the performance of PS does not scale to high-dimensional systems. To do so, compare the number of particles required for PS reliably characterize the distribution over Z1.T for a system of two and six coupled Lorenz systems. The systems are simulated with observation noise a, 0.01 but without process noise.
We implement a PS as specified by |Schön et al. (2011). Let it,n, Wt,n be the estimated system-state and partiele weight corresponding the nith particle at the ith timestep. In order test whether PS can reliably characterize the posterior distribution over T1.7, we measure the weighted root mean square error (RMSE), & as follows:
T
Win 12. æell2
EN
(10)
NT
Σ Wt,n
[=1
n=
If PS is reliable using N particles, we should see that EN is tightly distributed across random seeds.
Figure shows histograms of EN for random seeds appear to have values of εNv greater than that of the true posterior state distribution only the case of N 10. However, in the case of the higher-dimensional six-attractor system, not even 160 particles are sufficient to reliably characterize the true posterior, shown by the fact that the majority of random seeds have large values for ÊN-
[10, 40, 160, 640). For the two-attractor system, many
These results experimentally demonstrate that, for a sufficiently non-linear and partially observ- able system, the number of particles required to reliably characterize the posterior distribution over hidden states grows intractably with the dimension of the system. Since particle-based EM meth- ods for system-ide iaatem-idnttification are typically super-linear in complexity with respect to the number of particles (Schön et al. 2011), these methods are ill-suited high-dimensional problems.
5

---

Under review as a conference paper at ICLR 2020
3.1.3 CONVERGENCE OF SISL ON SIGH-DIMENSIOAL PROBLEMS
To demonstrate that SISL is capable of identifying high-dimensional systems, we show that we can estimate the dynamics of an 18 dimensional system of six coupled Lorenz attractors. Moreover, as the number trajectories provided to SISL increases, it converges to more accurate estimates. To test this claim, sample 2, 4, and 8 trajectories from system with parameters 8 and 0, 0.01. We randomly initialize each element the parameters being optimized (Ø ; PI:K; fl.K H) to within 10% of the their value in Ouve- We then SISL each batch, tracking the error in the estimated dynamics as training proceeds. We measure this error, which we call E(0), as follows:
(00 Ezzp((z)) [l)))())) Jo (3r)
(11)
In the learning step, we do not regularize to a prior and set pe = 0.
In Figure [
the results of this experi- ment for four random seeds for each batch size. We can that, the number of trajectories used in training increases, the error in the es- timated dynamics tends toward zero. Further- more, we see that SISL convergences monoton- ically local optimum in all cases. This ex- periment supports our claim that SISL is able to identify the parameters of a digh-dimensional, non-lineat and partially observable system that intractable for particle-based methods.
2 Trajectories
4 Trajectories
109
8 Trajectories
10-1
500
1000
SISL Epoch
Figure 2: Error in estimated dynamics as SISL trains on a varied number of trajectories from system of six coupled Lorenz attractors.
The experiments conducted thus far have demonstrated that SISL can learn unbiased pa- rameter estimates of -eatiiiittriiiiittic sys- tems, and can scale to high-dimonn nigh-dimensional prob-
lems for which particle-based methods intractable. In the next experiment, use SISL to characterize the effect of unobserved states on the dynamics of an aerobatic helicopter.
3.2 CHARACTERIZING AEROBATIC HELICOPTER DYNAMICS
Characterizing the dynamics of helicopter undergoing aggressive aerobatic maneuvers widely considered to be a challenging pustem-identification problem (Abbeel et al.) 2010 Punjani & Abbeel 2015). The primary challenge that the forces on the helicopter depend on the induced state of the fluid around it. The state of the fluid cannot be directly observed and its dynamics model unknown. Merely knowing the state of the helicopter and the control commands at a given time does not contain enough information to accurately predict the forces that act on it.
In order to address this issue, Punjani & Abbeel (2015) use an approach based on Takens theorem, which suggests that system's state can be reconstructed with a finite number of lagged-observations of it (Takens] 198T). Instead of attempting to estimate the unobserved fluid state, they directly learn mapping from a 0.3s long history of observed state measurements and control commands to the forces acting on the helicopter.
This approach is sensible, and is equivalent considering the past 0.5 of observations as the system's 'ss state. However, it can require a very large number of lagged observations to represent complex phenomena. Having such a high dimensional state can make the control design and state- estimation more complicated. To avoid large input dimensions. a trade-off between the duration of the history and sample frequency is necessary. This trade-off will either hurt the resolution of low-frequency content or will alias high-frequencies. We attempt to instead explicitly model the unobserved states affecting the system.
3.2.1 OBJECTIVE AND DATASET
The objective of this learning problem is predict y/, the helicopter's acceleration time t, from an input vector u, containing the current measured state of the helicopter (its velocity and rotation rates) and the control commands.

---

Under review as a conference paper at ICLR 2020
We use data collected by the Stanford Autonomous Helicopter Project (Abbeel et al. 2010). Tra- jectories are split into 10s long chunks and then randomly distributed into train, test, and validation sets according to the protocol established by Abbeel et al. (2010); Punjani & Abbeel (2015) and
summarized in Appendix A. The train, test and validation sets respectively contain 466, 100 and 101 trajectories of 500 time-steps each.
simple success metric on a given trajectory is the root mean squared prediction error, RMSE =
(measured)
(measured) is the measured force from the dataset, (pred) where V/2
Vt 2
(pred) is the force predicted by the model, and T is the number of time-steps in each trajectory. Vi
3.2.2 PREVIOUS WORK AND BASELINES
We first consider naive baseline that does not attempt to account for the time- varying nature of the fluid-state. We train a neural-network to map only the current helicopter state and control commands to the accelerations: Y/- NNo., (u), where NNo,, is neural-network with parameters 0z, We refer to this model as the naive model.
We also compare the work of Punjani & Abbeel (2015). They prediet y using time-history Ut- Hat of H lagged observations of the helicopter's 's measured state and control commands. This input is passed through Rel.U-activated neural network with a single hidden-layer combined with what they call a Quadratic Lag Model. As a baseline, we reproduce their performance with a single deep neural network 8/+ NNo, (u- with parameters 0n. We call this neural network model the H25 model. Both of these models can be trained via stochastic gradient descent to minimize the Mean- Sauun--qared-E-ie Errorrr--rrrrreeee (MSE) of their predictions for y. The optimization methodology for these models is described in Appendix A.2
As third baseline, we compare with cubbsacc--denttffiiee 1994]. We let Vt = N. (ux) be the prediction errors of the trained naive model. We use the MATLAB command n iidi to fit a linear dynamical system of the following form:
methods (Van Overschee & De Moor
F++1 = A,t, + Bu4;
Gt = C. + D.ut
(12)
Here, x Rdd the unobserved state with arbitrary dimension The learned parameters are 8, [A,, B,, C2 Dz]. We use state dimension of 10 and call this model the SID model. The 44sid algorithm scales super-linearly with the amount of data supplied, and thus we train on 10 randomly sampled subsets of 100 trajectories each, and report the distribution in prediction performance.
Particle- EM methods are not presented as baselines because they are intractable on problems with large state-spaces, shown in Section 3.1.2
3.2.3 NON-LINEAR UNOBSERVED STATE MODRL
Similar to the parameterization used for subpace-identificttnn we fit the prediction errors of the naive model using the following dynamical system:
I+41 = A + Bx UL = CCLEI + DNL + NN (Ie, us)
(13)
where NN is a neural network, and 6x. = [ANL, BAL. CNL. DrL. Fxu.) are the learned parameters. While learning, we assume that both process and observation noise are distributed with diagonal Gaussian covariance matrices I and I respectively. The values of a and a, are treated as hyperparmeters of SISL. Here as well, we use state dimension of 10 and call this model the NL model. The optimization methodology for this model is described in Appendix A.2
It should be noted that the system we learn need not actually correspond an interpretable model of the fluid-state, but only of time-varying hidden-states that are useful for predicting the accelerations of the helicopter. Expert knowledge of helicopter aerodynamics could be used to further inform gray-box model trained with SISL.
3.2.4 EVALUATION METHODOLOGY
The test RMSE of the naive and H25 models can be evaluated directly on the test trajectories using next-step prediction. However, the SID) and NI. models require estimate of the unobserved state
7

---

Under review as a conference paper at ICLR 2020
NL (ours)
Inverted, ve tia swesps
forward atsenays Aig
freenlll
(m/66)
demas2
free lel
steg and gr
onferiation
_wwope, ofeeracion_wwess
dodging demess
demoss
onientation, sweeps
(m//at)
Minh
Ecoecs
-5
oedging demest
freesalk, apgressre
dodging, demasi
dedging demau2
H25
lourst
(m//*)
demos]
SD
0.0
1.5
Emor (me
time
3
(a)
(b)
Figure 3: (a): Test performance of optimized models on various trajectories. Error bars SID represent the standard deviation of performance of the 10 trained models. (b): Predicted acceleration along axis z, y and in the body frame for one of the harder test set trajectories. Larger versions of these plots, including rotational accelerations, can be found in Appendix A.3
before making a prediction. The natural analog of next-step prediction is extended Kalman filtering (EKF), during which states are recursively predicted and corrected given observations. At given time-step, a prediction Dt is made using the current estimate of xy, and is used the computation of RRSE.E The state-estimate is then corrected with the measured j/.
3.2.5 RESULTS
Figure 3a shows the RMSE of the compared models on trajectories in the test-set. We see that the NL model is able consistently predict the accelerations on the helicopter with better accuracy than any of the other models. The naive model performs on average 2.9 times worse than the H25 model, and its results can be found in Appendix AA.3 The SID model notably outperforms the state-of-the- H25 model, suggesting that a large linear dynamical system can be used to approximate a non-linear and partially observable system (Korda & Mezić 2018). However, introducing non-linearity as in the NL model noticeably improves performance.
Figure 3b depicts the errors in prediction over sample trajectory the test-set. Here, we also see that the NL model is able to attenuate the vime-varying error present in predictions made by the H25, suggesting that has accurately characterized the dynamics unobserved, time-varying states.
This experiment validates the effectiveness of SISL to identify non-linear dynamical model of unobserved states that affect the forces acting an aerobatic helicopter.
4 CONCLUSIONS
This paper presented algorithm for system identification of non-linear systems given partial state observations. The algorithm optimizes system parameters given a time history of observations by iteratively finding the most likely state-history, and then using it to optimize the system parameters. The approach is particularly well suited for high-dimensional and nearly deterministic problems.
In simulated experiments on partially observed system of coupled Lorenz attractors, we showed that our algorithm can perform identification on a problem that particle-based EM methods are fundamentally ill-suited for. We also validated that our algorithm is an effective replacement for identification methods based on EM if the system is close to deterministic, but can yield biased parameter estimates if it is not. We then used our algorithm to model the time-varying hidden- states that affect the dynamics of an aerobatic helicopter. Our approach outperforms state-of-tt ffate-of-the-n the-art methods because it is able to fit large non-linear models unobserved states.
8

---

Under review as a conference paper at ICLR 2020
We aim to apply our algorithm to system identification problems in number of domains. There has recently been interest in characterizing the dynamics of aircraft with high aspect ratios, for which the difficultto-obseeve bending modes substantially impact dynamics. Additionally, the inability to measure friction forces in dynamic interactions involving contact typically stands in the way of system identification, and thus requires algorithms that are capable of identification under partial observation.
REFERENCES
Pieter Abbeel, Adam Coates, and Andrew Y Ng. Autonomous helicopter aerobatics through ap- prenticeship learning. The International Journal of Robotics Research, 29(13):1608-1639 2010.
Karl Johan Aström and Peter Eykhoff. System identification survey. Automatica, 7(2):123-162, 1971.
Coryn AL Bailer-Jones, David JC MacKay, and Philip J Withers. A recurrent neural network for modelling dynamical systems. Network: Computation in Neural Systems, 9(4):311-44,, 1998.
Pierre Bergé, Yves Pomeau, and Christian Vidal. Order within Chaos: Towards a Deterministic Approach to Turbulence. Wiley & Sons, 1984.
Stephen A Billings. Nonlinear System Identification: NARMAX methods in the time, frequency and spatio-temporal domains. John Wiley Sons, 2013.
Gilles Celeux and Gérard Govaert. A classification EM algorithm for clustering and two stochastic versions. Computational statistics & Data analysis, 14(3):11-332,, 1992.
Rick Cory and Russ Tedrake. Experiments in fixed-wing UAV perching. In AIAA Guidance, Navi- gation and Control Conference and Exhibit, pp. 7256, 2008.
Arthur P Dempster, Nan M Laird, and Donald B Rubin. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 39(1): 1-22, 1977.
Zoubin Ghahramani and Sam T Roweis. Learning nonlinear dynamical systems using an EM algo- rithm. In Advances Neural Information Processing Systems (NeurIPS), pp. 431-437, 1999.
Jayesh K Gupta, Kunal Menda, Zachary Manchester, and Mykel Kochenderfer. general frame- work for structured learning of mechanical systems. arXiv preprint arXiv:1902.08700 2019.
Chong-Shens Hou and Zhuo Wang. From model-based control data-driven control: Survey, classification and perspective. Information Sciences, 235:3-35, 2013.
Eric Jones, Travis Oliphant, Pearu Peterson, et al. SciPy: Open source scientific tools for Python, 2001-. htRLhttp:/ //tt::://ww scipy. org/
Nikolas Kantas, Arnaud Doucet, Sumeetpal S Singh, Jan Maciejowski, and Nicolas Chopin. On particle methods for parameter estimation in state-space models. Statistical Science, 30(3):328- 351, 2015.
Diederik Kingma and Jimmy Ba. Adam: method for stochastic optimization. In International Conference on Learning Representations, 2015.
Milan Korda and Igor Mezić. Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control. Automatica, 13:149-160, July 2018. doi: 10. 10.1016/j autoomataa 2018.03.046.
Kenneth Levenberg. A method for the solution of certain non-linear problems in least squares. Quarterly of Applied Mathematics, 162):164 July 1944. doi: 10.1090/qam/10666.
Lennart Ljung. System identification. Wiley Encyclopedia of Electrical and Electronics Engineer- ing, pp. 1-19, 1999.
Lennart Ijuug. Perspectives system identification. Annual Reviews in Control, 34:1-12, 2010.

---

Under review as a conference paper at ICLR 2020
Michael Lutter, Kim Listmann, and Jan Peters. Deep Lagrangian networks for end-to-end learning of energy-based control for under-actuated systems. In International Conference on Intelligent Robots and Systems (IROS), 2019a.
Michael Lutter, Christian Ritter, and Jan Peters. Deep Lagrangian networks: Using physics model prior for deep learning. In International Conference on Learning Reprexentations, 2019b. URL ftttps://openrewee
net/ BorumiiddBBk1ijjccaa
Donald W. Marquardt. An algorithm for least least-squares estimation of nonlinear parameters. Journal of the Society for Industrial and Applied Mathematics, 11(2):31-441, June 1963. doi: 10.1137/ 01 01100
Radford Neal and Geoffrey Hinton. view of the EM algorithm that justifies incremental, sparse, and other variants. In Learning in Graphical Models, pp. 355-368. Springer, 1998.
Camilo Ordonez. Nikhil Gupta, Brandon Reese, Neal Seegmiller, Alonzo Kelly, and Emmanuel G Collins Jr. Learning of skid-steered kinematic and dynamic models for motion planning. Robotics and Autonomous Systems, 95:207-2211 2017.
Ali Punjani and Pieter Abbeel. Deep learning helicopter dynamics models. In IEEE International Conference Robotics and Automation (ICRA), pp.
2015.
Thomas B Schön, Adrian Wills, and Brett Ninness. System identification of nonlinear state-space models. Automatica, 47(1))39-49, 2011.
Chris Snyder, Thomas Bengtsson, Peter Bickel, and Jeff Anderson. Obstacles to high-dimensional particle filtering. Monthly Weather Review, 136(12):4629 04494,4600000 2008.
Floris Takens. Detecting strange attractors in turbulence. pp. 366-381. Springer, 1981.
Dynamical Systems and Turbulence,
Peter Van Overschee and Bart De Moor. N4SID: Subspace algorithms for the identification of combined deterministic-stochasii systems. Automatica, 30(1):75-93, 1994.
HG Zimmermann and R Neuneier. Modeling dynamical systems by recurrent neural networks. WIT Transactions on Information and Communication Technologies, 25, 2000.
10

---

Under review as a conference paper at ICLR 2020
A APPENDIX
A.1 DATASET COLLECTION AND PREPROCESSING
In this work we use the dataset gathered by Abbeel et al. 2010 and available at http:://hhhl stanford. edu/ A gas-poweres helicopter was flown by professional pilot to collect large dataset of 6290s of flight. There are four controls: the longitudinal and lateral cyclic pitch, the tail rotor pitch and the collective pitch. The state is measured thanks an accelerometer, gyroscope. a magnetometer and vision cameras. (Abbel et al. (2010) provide the raw data, as well as states es- timates in the Earth reference frame obtained with extended Kalman smoothing. Following Punjani & Abbeel (2015) s protocol, we use the fused sensor data and downsample it from 100Hz to 50Hz.
From the Earth frame accelerations provided in the dataset, we compute body frame accelerations (minus gyroscopic terms) which are the prediction targets for our training. Using the notations from Punjani Abbeel 2015), we can write the helicopter dynamics in the following form:
C120
2wq
8 = F(,
(14)
C
v (n(,,)) (s,
where E R13 is the helicopter state consisting of its position v, quaternion-attitue q, linear velocity e, angular velocity w. and 8 E I to be the control command. C12 is the rotation-matrix from the body to Earth reference frame, and and fw are the linear and angular accelerations caused by acrodynamic forces, and are what we aim to predict.
This notation connects with the one used in paper in the following way:
We define the concatenation of all inputs the model, including the relevant state variables e and w and control commands 8.
We define as the output predicted, which would correspond to concatenation of fo and fw.
We define x as the vector of unobserved flow states to be estimated and not present in their model.
A.2 TRAINING HELICOPTER MODELS
Neural-networks in the naive and H25 models have eight hidden layers of size 32 each, and tanh non-lincarities. We optimize these models using an Adam optimizer (Kingma & Ba 2015) with harmonic learning rate decay, and mini-bateh size of 512.
The neural network in the NL model has two hidden layers of size 32 each, and tanh non-linearity. We train the NL model with SISL, using P: Pe 0.5, Tw o. 1.0, and use Adam optimizer to optimize Equation 77) in the learning step. The learning rate for dynamics parameters in 6N. is 5.0 10-1 and observation parameters in 0sa. is 1.0 x 10-3 For its relative robustness, we optimize Equation (6) using a non-linear least squares optimizer with a Trust-Region Reflective algorithm (Jones al. | (2001-) in the smoothing step. This step can be solved very efficiently by providing the solver with the block diagonal sparsity pattern of the Jacobian matrix.
To evaluate the test metric, running an EKF is required. The output of EKF depends on several user- esoriprovided parameters:
Tp: value of the initial state
E covariance of error on initial state
Q: covariance of process noise
R: covariance of observation noise
In this work, we assume that Q, R and En are all set to the identity matrix. mo is assumed to on all dimensions.
11

---

Under review as a conference paper at ICLR 2020
A well-tuned EKF with inaccurate initial state value converges to accurate estimations in only few time steps of transient behavior. Since the H25 model needs 25 past inputs to predict its first output prediction, drop the first 25 predictions from the EKF when computing RMSE, thereby omitting some of the transient regime.
A.3 FIGURES
Naive
turn_demos1
H25
NL (ours)
vertical_sweeps
SID
stop_and_gg
dodging_demosi
orientation_sweepp
forward_sideways_flig
tictocs
chaos
bodging_demose
free_fall
freestyle_gentle
circles
flips_loops
dodging_demos2
Freestyle_aggressivee
turn_demos3
turn_demos2
inverted vertical_sweeps
ercentation_sseps__wwee
dodging_demne
3
5
RMS Error [ms
2]
Figure 4: Test performance of optimized models on various trajectories
12

---

Under review as a conference paper at ICLR 2020
Figure 5: Predicted acceleration along axis and 2 in the body frame. For subspace identification and models trained by SISL, this plot requires running extended Kalman filter. These figures can be reproduced for any other trajectory with the included code.
SID
NL (ours)
time (s)
H25
datasct
25
(120/ 2)
(va//" 2)
Dly
(m//²) 2
dx
13

---

Under teview as a conference paper at ICLR 2020
Figure 6: Predicted circular accelerations around axis y and in the body frame. For subspace identification and models trained by SISL, this plot requires running an extended Kalman filter. These figures can reproduced for any other trajectory with the included code.
SID
NL (ours)
time (s)
H25
2
2
dataset
0
10
10
10
10
50
50
(mafs 5)
(m// (m/s²)) 2)
(m/ )
Way
14


Reviews:
[{"review_id": "rygPOM5kqH", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper proposes a novel way to address the well studied system identification problem. More specifically, this paper focuses on non linear systems with partial observations. To this aim, the authors propose the \"System-indetification via Iterative Smoothing and Learning\" (SISL) algorithm. The gist of the approach is to formulate an ML problem to estimate the system model parameters from a parametric family of models to estimate both the dynamical system evolution (state updates) and the observation process. The exact optimization for this ML problem is not tractable. This work proposes to solve a surrogate point-estimate  ML problem. This optimization is solve using an alternate approach (between smoothing and learning). \n\nThe paper also presents experimental results both on synthetic  and real data aiming to demonstrate that SISL provides unbiased estimates in deterministic settings. The synthetic setting is also used to motivate the need for methods able to handle large system dimensionality for non-linear and partially observable systems. \n \nThe real data setting compares favorably the SISL approach with well established baselines. \n\nI believe that the paper can be improved in the following ways: \n1) There are intermediate steps in the derivation of the maths that could be better explained. For example, in equation 2, the “Markovian” aspect of the system is key to derive the formula (x_{t+1} only depends on x_{t}). Equation 4 is not really derived from Equation 3. The factorization shown in Equation 2 is relevant for Equation 4 although the 4) is for the surrogate objective function. \n\nAlso, I did not understand why the naive approach is not plotted on Figure 3 (a and b) (but it comes back in the appendix).  Is it just an oversight?", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "BJehPT0uKB", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "\n\n#######################\nRebuttal Response:\nThanks for these clarifications and updating the paper. I adapted my score to weak accept. To increase my score to accept, I would like to see evaluations of control experiments rather than just comparing the MSE. \n\n#######################\nReview: \n\nSummary: \nThe paper introduces an iterative learning scheme to perform system identification on partially observable systems. The proposed approach iteratively optimizes the unobserved state and the model parameters. The proposed approach is demonstrated on a 'high'-dimensional Lorenz oscillator and existing quadcopter data. Furthermore, the authors show that existing approaches don't scale to 'high' dimensions. \n\nConclusion;\nThe paper introduces a simple iterative scheme to learn a model with a hidden state. The approach is sound and reasonable and the paper is well written but I am uncertain whether this paper fits the scope of ICLR. Maybe R:SS would be a better fit? Furthermore, the experiments can be improved, which is especially important as this paper does not propose a fundamental new approach. It would be really interesting to see this model controlling an actual quadcopter. My current vote is borderline reject but can be improved by further highlighting why quadcopter dynamics are only partially obsevable and further experiments.\n\nExperiments:\nMy main concern with the paper is the experiments of more complex systems. The initial experiments demonstrating that the approach works and previous methods don't scale to high dimensions are very good and highlight important points. However, the more complex experiments might be not indicative of the performance. Especially for the quadcopter experiments, the proposed method only learns a non-linear noise model while the state transition is still linear and compares it against the purely linear case. Both methods seem qualitatively comparable and the quantitative improvement is only marginal. Therefore, I am doubtful whether the quadcopter data really requires to incorporate the partially observable model and I would prefer to see an evaluation of the proposed iterative scheme using a linear model. Within your evaluation you should compare the different optimization approaches (i.e., naive vs iterative) using the same model representation as your proposed algorithm only proposes a new optimization scheme but is not fixed to a representation. Furthermore, it would be interesting to see a non-linear state model and a non-linear observation model. Does the iterative approach also learn this model? Learning the robot dynamics model with contacts would be the best experiment to validate your approach, as this is only partially observable using the joint measurements and is non-linear. Could you setup an experiment with a simulated two-link robot, generate data containing samples with & without contact and learn the model just using q, qd, qdd, and tau. This experiment would clearly highlight that your iterative scheme can identify a partially observable model.  \n\nMinor Comments:\n\"We let ̃yt=yt−NNθn(ut)be the prediction errors of the trained naive model.\" -> it is a bit unclear whether the linear model approximates just the noise or the complete model? From the following text I assume that the complete model is modelled as linear model. Then the introduction of the naive model seems unnecessary as the naive model is not used. Could you please clarify this point?", "extra_scores": {"experience_assessment": "I have published one or two papers in this area.", "review_assessment:_checking_correctness_of_experiments": "I carefully checked the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I carefully checked the derivations and theory."}, "extra_text": {}}, {"review_id": "B1eBS0jm5B", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This is a strait forward paper that aims to address the lack of scalability of prior methods on system  identification.  The idea is simple yet the assumptions required are significant and restrict the applicability of the proposed algorithm to systems that are closed to deterministic.  Below are some pros and cons for this paper:\n\nCons:\n(+) The proposed method is simple and it could be  quickly adopted by robotics researchers and practitioners\n(+) The method seems to outperform previous algorithms\n\nPros:\n(-) The idea of working with deterministic systems is restrictive.  While the process noise may appear in some robotics applications to be small  I see no way of how one could reduce the effect of the observation noise.    Also in the context of systems that navigate in terrestrial environments  the process noise is definitely more present that in the case of flying vehicles. \n\n(-) It would definitely strengthen the paper if the authors could elaborate more on the concepts of controllability, observability and how these concepts relate to their method etc. How does the existing methods handles these properties? How is the performance of the method affected  when the true system is not controllable but it is stabilizable?\n\n\n(-) There seems to be no assumption on the shape of the trajectories. What kind of trajectories you need to perform system identification? Can you mathematically represent these assumptions? In the adaptive control literature there exists conditions related to \"sufficient rich signals\" see the book by Ioannou on Robust Adaptive Control. These are signals that can excite dynamics so that to ensure proper identification.  \n\n\n\n", "extra_scores": {"experience_assessment": "I have published one or two papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper is about nonlinear system identification in an EM-style learning framework. The idea is to use nonlinear programming for the E step (finding a MAP estimate) and then refine the model parameters. In flavor, this approach is similar to the work by Roweis and Ghahramani. 

However, this paper does not offer any new insights whatsoever and the (very short) methods section arrives at proposing to compute the maximum a posteriori estimate (eq. 5). While the motivation for this given in the paper is a bit hard to understand it is of course a very well-known and useful estimator. Besides the maximum likelihood estimator this is one of the most commonly used point estimators, see any textbook on statistical signal processing. There has been quite a bit of work in the signal processing community over the last 10 years, and a good overview can be found here:
https://web.stanford.edu/~boyd/papers/pdf/rt_cvx_sig_proc.pdf
This should give evidence that this is indeed a standard way of solving the problem and it does work really well. Given that we have so fast and good optimizers these days it is common to solve Kalman filtering/smoothing problems via this optimization problem.
The paper does not contain any analysis at all. The experiments do of course show that the method works (when there is low noise). Again, we know very well that the MAP estimate is a decent estimator for unimodal problems. The MAP estimator can also be made to work well for noisy situations.

As for the comments that the sequential Monte Carlo methods do not work in higher dimensions that is indeed true. However, there are now algorithms that work in much higher dimensions than those considered by the authors of this paper, e.g.
https://ieeexplore.ieee.org/document/8752074
which also contains an up-to-date survey on the topic. Furthermore, when it comes to particle smoothing there are also much more efficient smoothers than 10 years ago. The area of particle smoothing has also evolved rapidly over the past years. 

Summary:
The paper makes use of the well-known MAP estimator for learning nonlinear dynamical systems (states and parameters). This is by now a standard technique in signal processing. There are several throw-away comments on SMC that are not valid and that are not grounded in the intense research of that field over the past decade.

Decision: Reject