# Dynamical models for metabolomics data integration

Polina Lakrisenko  Daniel Weindl<EMAIL_ADDRESS>Institute of Computational Biology, Helmholtz Zentrum München, Neuherberg
85764, Germany  Center for Mathematics, Technische Universität München,
Garching 85748, Germany

###### Abstract

As metabolomics datasets are becoming larger and more complex, there is an
increasing need for model-based data integration and analysis to optimally
leverage these data. Dynamical models of metabolism allow for the integration
of heterogeneous data and the analysis of dynamical phenotypes. Here, we
review recent efforts in using dynamical metabolic models for data
integration, focusing on approaches that are not restricted to steady-state
measurements or that require flux distributions as inputs. Furthermore, we
discuss recent advances and current challenges. We conclude that much progress
has been made in various areas, such as the development of scalable simulation
tools, and that, although challenges remain, dynamical modeling is a powerful
tool for metabolomics data analysis that is not yet living up to its full
potential.

###### keywords:

metabolic modeling , kinetic modeling , mechanistic modeling , data
integration , dynamical model , metabolomics

††journal: Current Opinion in Systems Biology

## Introduction

Metabolism is a key determinant of cellular behavior, and metabolomics
approaches are being applied in a wide range of domains [1]. Nowadays,
metabolomics experiments yield rich datasets, which poses new challenges for
their analysis. Current metabolomics assays provide information on hundreds to
thousands of metabolites and large numbers of samples can be measured within
reasonable times [2]. Since holistic interpretation of the resulting datasets
using “mental models” becomes impossible, computational models are required
and increasingly used for data integration and interpretation. However,
observations from individual experiments are still commonly analyzed
independently and in a qualitative manner, instead of being integrated in an
overarching formal quantitative model.

A variety of metabolic modeling frameworks have emerged, which allow for the
integration of varying types of data, each with their specific advantages and
disadvantages. A broad overview of these existing metabolic modeling
approaches is provided in recent reviews [3, 4, 5]. In this article, we focus
on metabolomics data integration using dynamical models. These models, usually
specified in the form of ordinary differential equations (ODEs), are
particularly appealing because: 1) they allow for both the integration of
various types of metabolomics data and other types of data; 2) they allow for
the analysis of inherently dynamical phenotypes; and 3) they can provide
quantitative dynamical information on latent quantities, such as metabolic
fluxes or compartment-resolved metabolite levels. On the downside, dynamical
models come with a larger number of a priori unknown parameters and are
computationally more challenging, and therefore, are less scalable than other
modeling approaches. When only considering steady-state measurements, these
issues can be circumvented to a large extent, as exemplified by the recently
developed K-FIT algorithm [6]. However, one is often interested in transient
behavior, or a case where a steady-state cannot be attained. Therefore,
although such methods have provided great insights in many applications, they
are not discussed further here. We also exclude methods that require a known
flux distribution or a biochemical objective function as input. For an
overview of these methods and their applications, see, for example, [5, 7, 8].

In the following, we give a brief overview of different types of metabolomics
data, review recent examples of metabolomics data integration and analysis
using dynamical models, and discuss major challenges and recent advances.

## Metabolomics data

Metabolomics measurements are usually performed using mass spectrometry,
nuclear magnetic resonance, or enzymatic assays, with mass spectrometry being
the most widespread. Metabolomics datasets can be very rich and diverse,
depending on the analytical platforms and experimental protocols. There are
the classical label-free metabolomics approaches, which assess metabolite
levels, and there are stable-isotope-assisted approaches, which determine
isotopic enrichment after application of stable-isotope-labeled tracers [1].
Recent mass spectrometry methods can quantify isotopic enrichment of more than
100 metabolites in a single run [9]. Various types of samples are analyzed,
such as complete cell populations and extracellular media, and in some cases
compartment-specific metabolite pools are accessible [10]. Single-cell
metabolomics assays have also emerged [11, 12]. Measurements can be time-
resolved or for single timepoints, quantitative or semi-quantitative. Due to
this heterogeneity, various metabolic modeling approaches have been developed
that are, more or less, tailored or restricted to specific types of data [3].
However, in many studies, a combination of different types of data is acquired
and integration of all these data is possible with dynamical models [13, 14].

## Applications

Model-based data analysis aims at analyzing experimental observations in the
light of the current understanding of the observed system, as encoded in the
model. Experimental observations are either used as model inputs or to infer
model parameters. The behavior of the parameterized model is subsequently
analyzed to obtain biological insights. Here, we focus mainly on the case of
inferring model parameters from experimental data. Typical steps comprise 1)
data acquisition, 2) model construction, 3) parameter inference, 4)
uncertainty analysis, and 5) analyzing model fit and model-derived predictions
(Fig. 1).

Figure 1: Typical steps for data integration using dynamical models. A model
is developed based on current knowledge (brackets indicate optional
components). Experimental observations of various types are used as model
inputs or to estimate model parameters. If parameters are estimated, the model
fit and parameter intervals are assessed. Competing hypotheses or structural
uncertainty can be encoded in different candidate models to select the most
plausible one given the data. The parameterized model is used in various
downstream analyses to derive biological insights. Examples include
sensitivity analysis or metabolic control analysis to determine which
parameters have the highest impact on reaction fluxes or metabolite
concentrations. In silico experiments can be performed to derive testable
hypotheses and to design validation experiments. Considering uncertainties is
crucial in all analyses to derive meaningful conclusions.

Dynamical metabolic models have been applied in various contexts. Recurrently
pursued objectives are 1) understanding dynamical processes [15, 16], for
example, through comparison of competing models [17]; 2) inferring control
mechanisms and rate-limiting steps [18, 13]; and 3) leveraging those to push
some system of interest in specific directions, for example, for strain
optimization or drug target identification) [19, 20, 14].

Berndt et al. [16] developed a comparably large dynamical model of liver
metabolism to, among others, study the response to various perturbations. The
model describes the regulation of enzyme activities by allosteric effectors,
hormone-dependent reversible phosphorylation, and variable protein abundances.
A subset of parameters was estimated from measurements of intracellular and
extracellular metabolite levels under 25 different experimental conditions.
Carter et al. [17] performed model selection among four candidate models,
calibrated on time-resolved measurements of metabolite levels, to infer the
most likely inhibitory mechanism in a drug-drug interaction. Feldman-Salit et
al. [18] trained a dynamical model of sulfur assimilation in Arabidopsis
thaliana on steady-state metabolite measurements, which allowed them to infer
dynamical control patterns. Millard et al. [13] analyzed time-resolved
measurements of biomass, metabolite levels, and isotopic enrichment using a
coarse-grained dynamical model that links glucose uptake to acetate metabolism
and growth. Their analysis improved the understanding of control and
regulation of acetate overflow and suggested the existence of a yet unknown
regulatory program. Álvaro Marín-Hernández et al. [19] applied a dynamical
model of glycolysis, the pentose phosphate pathway, and glycogen metabolism of
cancer and non-cancer cells, with parameters inferred from various enzyme
assays, to derive suggestions for new therapeutic targets. Ou et al. [20]
applied a dynamical model, which integrated experimental proteomics,
metabolomics, and fermentation kinetics data, to identify key regulators and
guide metabolic engineering of n-butanol biosynthesis. Their comparison of the
results from the dynamic model to those derived from a static model indicated
higher accuracy of the former. Ramos et al. [14] integrated both extracellular
and intracellular time-resolved metabolite measurements from different
experimental settings using a dynamical model comprising 33 state variables
and describing central carbon metabolism and cell growth. They identified
different physiological states and performed in silico experiments to
demonstrate how use of their model can improve cell line engineering and
medium design. Moon et al. [21] used a dynamical model to predict the dynamics
of mitochondrial NADPH concentration and NADPH/NADP+ ratio in response to
oxidative stress. Parameters of the model were estimated from time-resolved
NADPH measurements. Lövfors et al. [15] developed a model describing hormonal
regulation of lipolysis. The model, comprising 15 state variables, was
calibrated to time-resolved in vivo and in vitro measurements of metabolite
and phosphoprotein abundance after different stimulations. Horvath et al. [22]
calibrated an exceptionally comprehensive model of E. coli cell-free protein
synthesis, comprising 148 metabolites, based time-resolved on metabolite and
protein product measurements. An ensemble model was subsequently used to
identify the pathways with high influence on product yield. Yilmaz et al. [23]
developed a coarse-grained model of protein synthesis in Chinese hamster ovary
cells based on elementary flux modes. The model was calibrated to time-
resolved measurements of cell density and levels of extracellular metabolites
and protein product.

These application examples demonstrate the potential of dynamical models to
answer diverse research questions, by integrating and jointly analyzing
various types of data. However, many current applications of dynamical models
that integrate time-resolved measurements only cover small parts of the
metabolic network or use very coarse-grained representations. This may be due
to various challenges discussed in the following sections.

## Constructing kinetic models

Depending on the available data and research question, the model scope has to
be defined and a model needs to be constructed. Model construction includes
choosing and specifying 1) the network topology, and 2) kinetic rate laws.
Retrieving or inferring model parameter values is covered in the following
section.

Genome-scale metabolic reconstructions are a valuable resource for deriving
the model topology. Such knowledgebases that describe the metabolic
capabilities of an organism have been created for many species through
automated reconstruction or extensive manual curation through community
efforts [24]. For example, the recently published Human1 [25] unifies two
lineages of human genome-scale metabolic reconstructions and provides a good
example of open and transparent curation. Although such reconstructions are
continuously improved, structural uncertainty remains, for example due to
enzymes of unknown function, enzyme promiscuity, or unclear subcellular
localization of enzymes. Automated model reduction algorithms for generating
targeted kinetic models from genome-scale metabolic reconstructions have been
devised that can significantly speed-up model development [26, 27]. However,
care has to be taken, that the (over-)reduced model complexity does not skew
analysis results Hameri et al. [28].

If stable-isotope-labeling data are to be analyzed, the model not only needs
to account for the pools and reactions of unlabeled metabolites, but also for
those of the additionally occurring isotopic isomers (isotopologues). Deriving
such isotopologue reaction networks requires knowledge of the mapping of
substrate atoms to product atoms, for all reactions involving potentially-
labeled compounds. A number of algorithms for automatically deriving such atom
mappings have been developed (for example, [29]) and recent genome-scale
metabolic reconstructions also include atom mappings [30]. However, correct
atom mappings are crucial for modeling, and manual curation is still required,
which is tedious for bigger models. The resulting isotopologue reaction
networks can, depending on the network and tracer, be vastly larger than the
usual metabolic networks due to the combinatorial complexity of stable isotope
incorporation.

Having established a reaction network, mathematical expressions for describing
reaction rates need to be chosen. Various approaches for specifying rate laws
exist (reviewed in [8]). The choice of rate laws is a trade-off between
“mechanisticness” on the one hand and model complexity and number of
parameters on the other hand. Using a thermodynamics-based formalism [31, 32]
helps to create physically feasible models. Furthermore, such a
parameterization can reduce the number of unknown parameters and simplify
parameter estimation. The required thermodynamic parameters can be derived
from experiments or approximated, by group contribution methods or quantum-
mechanical simulations [33].

## Determining model parameters

Dynamical models usually come with a comparably high number of parameters
whose values are not known a priori. Some parameter values can be retrieved
from databases, while others need to be estimated. Often, a combination of
both approaches is used.

Commonly used databases of experimentally determined kinetic parameters
include BRENDA [34] and SABIO-RK [35]. Additionally, the BioModels Parameters
database [36] provides structured access to parameter values that are used in
models contained in the BioModels repository. If parameter values for a
specific organism or experimental conditions are not available, values from a
closely related species or settings may still be useful [16]. Datanator
provides a simple interface to find such related parameter values based on
various similarity measures [37]. However, using parameters determined by in
vitro assays or other makeshift parameter values comes with the caveat that
they may not optimally reflect the in vivo situation [5]. As an alternative to
taking these the parameter values from databases as true values, they can
instead be used as prior information during parameter estimation [8, 18].

Given sufficiently informative data, model parameters can be inferred, for
example, via optimization- [16, 17, 18, 13, 19, 20, 14, 21, 15, 23] or
sampling-based [22] approaches. However, this is computationally costly as,
for most cases, it involves thousands to millions of numerical ODE
simulations. Nevertheless, it has been demonstrated that parameter estimation
is computationally tractable for large dynamical models. For example, through
leveraging scalable algorithms, optimization-based maximum likelihood
estimation was shown to be possible for a dynamical model of signaling
pathways comprising over 1000 states and 4000 unknown parameters [38]. In some
cases, the structure of the optimization problem can be exploited to further
improve convergence and reduce computational costs. This has been demonstrated
for optimization-based parameter inference from relative measurements [39],
which are quite common in metabolomics datasets.

Scalable and highly optimized algorithms for model simulation and sensitivity
analysis have been made available through easy-to-use toolboxes [40, 41] and
there exists a wide variety of algorithms and tools for parameter inference
[42, 43]. Many of these tools are able to exploit increasingly available high-
performance computing resources, which is key for moving towards larger
dynamical models [44, 39, 38]. The use of community standards for specifying
models and parameter estimation problems [45, 46, 47] gives easy access to a
majority of these tools.

Sometimes model parameters are estimated independently and only later included
in the full model [19]. Although this can be computationally cheaper, it might
not yield the best fit to the data.

## Dealing with uncertainty

Independently of the algorithm employed for parameter estimation, parameter
estimates will be subject to uncertainty. This parameter uncertainty
propagates to prediction uncertainty and can result in false conclusions [48].
Uncertainty analysis is, therefore, crucial.

Different methods to determine parameter confidence intervals exist. Commonly
applied methods include those based on the Fisher information matrix (FIM),
the profile likelihood approach, or sampling-based procedures [49]. The choice
of method is a trade-off between accuracy and computational complexity, with
FIM-based methods being the cheapest computationally, but least accurate, and
sampling being the most expensive but also most informative. Profile
likelihood is often the method of choice, since it provides accurate results
and is more scalable than sampling-based methods. Furthermore, as opposed to
other methods, profile likelihood can be applied in spite of, and to detect,
non-identifiable parameters [49]. Implementations of these algorithms are
available through several easy-to-use software tools [42].

Methods that are similar to those used for parameter uncertainty analysis are
applied for prediction uncertainty analysis [48, 49]. Sampling-based
approaches account for the prediction uncertainty by design, as a sample from
the prediction posterior is acquired while generating a sample from the
parameter posterior. Where sampling or prediction profile likelihood are
computationally too expensive, a cheaper ensemble-based approach can be
applied instead [48]. During parameter estimation, multiple parameter vectors
will be obtained that result in similarly good model fits. An ensemble of
models can, for example, be built from these parameter vectors and the spread
of the resulting simulations can be used as an estimate for prediction
uncertainty [48].

Recent applications of dynamical models of metabolism include examples of
identifiability analysis [50, 17], sampling-based assessment of parameter and
prediction uncertainties [13], as well as other approaches uncertainty
analyses [15].

It is important to note that large parameter uncertainties do not have to
manifest in large prediction uncertainties, but may still allow for deriving
robust predictions [18, 38, 44]. For example, Feldman-Salit et al. [18]
generated an ensemble of models that reproduce the available data equally well
and used it to derive predictions that were consistent across the ensemble.

## Conclusions

Increasingly detailed metabolic reconstructions for more and more organisms
provide the basis for constructing genome-scale models, but also for deriving
targeted sub-models. Efficient implementations of scalable algorithms enable
simulation and sensitivity analysis for dynamical models with as many as a few
thousand state variables, even on personal computers. Adopting community
standard for specifying models or parameter estimation problems not only
facilitates reproducibility and reusability of models and data [47], but also
grants easy access to a wide ecosystem of tools for model simulation,
parameter inference and further analyses.

Computationally challenging parameter estimation is often considered the major
bottleneck for data-based dynamical (metabolic) modeling [5]. It is true that
for very large models this still is intractable and more scalable algorithms
are required. Results from other fields of application suggest that a better
leveraging of both existing tools and computational resources can facilitate
work with models and datasets that are an order of magnitude larger than those
currently used in the field of dynamical modeling of metabolism.

Aside from any algorithms, successful model-based data analysis depends on
matching models and data. To obtain informative data, experiments should
ideally be designed with a specific modeling goal in mind, and the available
data should be taken into account when building models. Nevertheless,
parameter uncertainty will always remain. Robust predictions may still be
possible despite large parameter uncertainties. Therefore, dynamical modeling
should not be ruled out prematurely due to supposedly too sparse data.
However, awareness and assessment of parameter and prediction uncertainties is
crucial, although often neglected.

In summary, we feel that dynamical modeling is a valuable, yet underused, tool
for integrative metabolomics data analysis that could help to derive a more
comprehensive, quantitative, and dynamical understanding of metabolism in many
applications.

## Acknowledgements

This work was supported by the German Federal Ministry of Education and
Research (BMBF) within the e:Med funding scheme (junior research alliance
PeriNAA, grant no. 01ZX1916A).

## References

  * Jang et al. [2018] C. Jang, L. Chen, J. D. Rabinowitz,  Metabolomics and isotope tracing,  Cell 173 (2018) 822–837. doi:10.1016/j.cell.2018.03.055.
  * Liu et al. [2019] X. Liu, L. Zhou, X. Shi, G. Xu,  New advances in analytical methods for mass spectrometry-based large-scale metabolomics study,  TrAC Trends in Analytical Chemistry 121 (2019) 115665. doi:10.1016/j.trac.2019.115665.
  * Volkova et al. [2020] S. Volkova, M. R. A. Matos, M. Mattanovich, I. Marín de Mas,  Metabolic modelling as a framework for metabolomics data integration and analysis,  Metabolites 10 (2020). doi:10.3390/metabo10080303,
** Broad review of metabolic modelling approaches for data integration.

  * Antoniewicz [2021] M. R. Antoniewicz,  A guide to metabolic flux analysis in metabolic engineering: Methods, tools and applications,  Metabolic Engineering 63 (2021) 2–12. doi:10.1016/j.ymben.2020.11.002.
  * Foster et al. [2021] C. J. Foster, L. Wang, H. V. Dinh, P. F. Suthers, C. D. Maranas,  Building kinetic models for metabolic engineering,  Current Opinion in Biotechnology 67 (2021) 35–41. doi:10.1016/j.copbio.2020.11.010.
  * Gopalakrishnan et al. [2020] S. Gopalakrishnan, S. Dash, C. Maranas,  K-fit: An accelerated kinetic parameterization algorithm using steady-state fluxomic data,  Metabolic Engineering 61 (2020) 197–205. doi:10.1016/j.ymben.2020.03.001.
  * Strutz et al. [2019] J. Strutz, J. Martin, J. Greene, L. Broadbelt, K. Tyo,  Metabolic kinetic modeling provides insight into complex biological questions, but hurdles remain,  Current Opinion in Biotechnology 59 (2019) 24–30. doi:10.1016/j.copbio.2019.02.005.
  * Saa and Nielsen [2017] P. A. Saa, L. K. Nielsen,  Formulation, construction and analysis of kinetic models of metabolism: A review of modelling frameworks,  Biotechnology Advances 35 (2017) 981–1003. doi:10.1016/j.biotechadv.2017.09.005.
  * Shi et al. [2020] X. Shi, B. Xi, P. Jasbi, C. Turner, Y. Jin, H. Gu,  Comprehensive isotopic targeted mass spectrometry: Reliable metabolic flux analysis with broad coverage.,  Analytical chemistry 92 (2020) 11728–11738. doi:10.1021/acs.analchem.0c01767.
  * Nonnenmacher et al. [2017] Y. Nonnenmacher, R. Palorini, A. F. d’Herouël, L. Krämer, M. Neumann-Schaal, F. Chiaradonna, A. Skupin, A. Wegner, K. Hiller,  Analysis of mitochondrial metabolism in situ: Combining stable isotope labeling with selective permeabilization,  Metabolic Engineering 43 (2017) 147–155. doi:10.1016/j.ymben.2016.12.005.
  * Ali et al. [2019] A. Ali, Y. Abouleila, Y. Shimizu, E. Hiyama, S. Emara, A. Mashaghi, T. Hankemeier,  Single-cell metabolomics by mass spectrometry: Advances, challenges, and future applications,  TrAC Trends in Analytical Chemistry 120 (2019) 115436. doi:10.1016/j.trac.2019.02.033.
  * Thiele et al. [2019] C. Thiele, K. Wunderling, P. Leyendecker,  Multiplexed and single cell tracing of lipid metabolism,  Nature Methods 16 (2019) 1123–1130. doi:10.1038/s41592-019-0593-6.
  * Millard et al. [2021] P. Millard, B. Enjalbert, S. Uttenweiler-Joseph, J.-C. Portais, F. Létisse,  Control and regulation of acetate overflow in Escherichia coli,  eLife 10 (2021) e63661. doi:10.7554/eLife.63661.
  * Ramos et al. [2020] J. R. C. Ramos, A. G. Rath, Y. Genzel, V. Sandig, U. Reichl,  A dynamic model linking cell growth to intracellular metabolism and extracellular by-product accumulation,  Biotechnology and Bioengineering 117 (2020) 1533–1553. doi:10.1002/bit.27288.
  * Lövfors et al. [2020] W. Lövfors, J. Ekström, C. Jönsson, P. Strålfors, G. Cedersund, E. Nyman,  A multi-level model analysis of lipolysis and fatty acid release from adipocytes in vitro and from adipose tissue in vivo,  bioRxiv (2020). doi:10.1101/2020.12.18.423229,
* Model of lipolysis integrating time-resolved metabolite measurements.
  * Berndt et al. [2018] N. Berndt, S. Bulik, I. Wallach, T. Wünsch, M. König, M. Stockmann, D. Meierhofer, H.-G. Holzhütter,  HEPATOKIN1 is a biochemistry-based model of liver metabolism for applications in medicine and pharmacology,  Nature Communications 9 (2018) 2386. doi:10.1038/s41467-018-04720-9.
  * Carter et al. [2020] S. J. Carter, A. S. Ferecskó, L. King, K. Ménochet, T. Parton, M. J. Chappell,  A mechanistic modelling approach for the determination of the mechanisms of inhibition by cyclosporine on the uptake and metabolism of atorvastatin in rat hepatocytes using a high throughput uptake method,  Xenobiotica 50 (2020) 415–426. doi:10.1080/00498254.2019.1652781.
  * Feldman-Salit et al. [2019] A. Feldman-Salit, N. Veith, M. Wirtz, R. Hell, U. Kummer,  Distribution of control in the sulfur assimilation in Arabidopsis thaliana depends on environmental conditions,  New Phytologist 222 (2019) 1392–1404. doi:10.1111/nph.15704,
* Ensemble-based analysis of control patterns.
  * Álvaro Marín-Hernández et al. [2020] Álvaro Marín-Hernández, J. C. Gallardo-Pérez, M. A. Reyes-García, M. Sosa-Garrocho, M. Macías-Silva, S. Rodríguez-Enríquez, R. Moreno-Sánchez, E. Saavedra,  Kinetic modeling of glucose central metabolism in hepatocytes and hepatoma cells,  Biochimica et Biophysica Acta (BBA) - General Subjects 1864 (2020) 129687\. doi:10.1016/j.bbagen.2020.129687.
  * Ou et al. [2020] J. Ou, T. Bao, P. Ernst, Y. Si, S. D. Prabhu, H. Wu, J. J. Zhang, L. Zhou, S.-T. Yang, X. M. Liu,  Intracellular metabolism analysis of clostridium cellulovorans via modeling integrating proteomics, metabolomics and fermentation,  Process Biochemistry 89 (2020) 9–19. doi:10.1016/j.procbio.2019.10.032.
  * Moon et al. [2020] S. J. Moon, W. Dong, G. N. Stephanopoulos, H. D. Sikes,  Oxidative pentose phosphate pathway and glucose anaplerosis support maintenance of mitochondrial nadph pool under mitochondrial oxidative stress,  Bioengineering & Translational Medicine 5 (2020) e10184. doi:10.1002/btm2.10184.
  * Horvath et al. [2020] N. Horvath, M. Vilkhovoy, J. A. Wayman, K. Calhoun, J. Swartz, J. D. Varner,  Toward a genome scale sequence specific dynamic model of cell-free protein synthesis in Escherichia coli,  Metabolic Engineering Communications 10 (2020) e00113. doi:10.1016/j.mec.2019.e00113,
** Data integration using a large model of E. coli cell-free protein
synthesis.

  * Yilmaz et al. [2020] D. Yilmaz, S. J. Parulekar, A. Cinar,  A dynamic EFM-based model for antibody producing cell lines and model based evaluation of fed-batch processes,  Biochemical Engineering Journal 156 (2020) 107494. doi:10.1016/j.bej.2020.107494.
  * Fang et al. [2020] X. Fang, C. J. Lloyd, B. O. Palsson,  Reconstructing organisms in silico: genome-scale models and their emerging applications,  Nature Reviews Microbiology 18 (2020) 731–743. doi:10.1038/s41579-020-00440-4.
  * Robinson et al. [2020] J. L. Robinson, P. Kocabaş, H. Wang, P.-E. Cholley, D. Cook, A. Nilsson, M. Anton, R. Ferreira, I. Domenzain, V. Billa, A. Limeta, A. Hedin, J. Gustafsson, E. J. Kerkhoven, L. T. Svensson, B. O. Palsson, A. Mardinoglu, L. Hansson, M. Uhlén, J. Nielsen,  An atlas of human metabolism,  Science Signaling 13 (2020). doi:10.1126/scisignal.aaz1482.
  * van Rosmalen et al. [2021] R. van Rosmalen, R. Smith, V. Martins dos Santos, C. Fleck, M. Suarez-Diez,  Model reduction of genome-scale metabolic models as a basis for targeted kinetic models,  Metabolic Engineering 64 (2021) 74–84. doi:10.1016/j.ymben.2021.01.008,
* Automated reduction of genome-scale models.
  * Masid et al. [2020] M. Masid, M. Ataman, V. Hatzimanikatis,  Analysis of human metabolism by reducing the complexity of the genome-scale models using redHUMAN,  Nature Communications 11 (2020) 2821. doi:10.1038/s41467-020-16549-2,
* Workflow for automated generation of reduced models of human metabolism.
  * Hameri et al. [2021] T. Hameri, G. Fengos, V. Hatzimanikatis,  The effects of model complexity and size on metabolic flux distribution and control: case study in escherichia coli,  BMC Bioinformatics 22 (2021) 134. doi:10.1186/s12859-021-04066-y.
  * Jaworski et al. [2019] W. Jaworski, S. Szymkuć, B. Mikulak-Klucznik, K. Piecuch, T. Klucznik, M. Kaźmierowski, J. Rydzewski, A. Gambin, B. A. Grzybowski,  Automatic mapping of atoms across both simple and complex chemical reactions,  Nature Communications 10 (2019) 1434. doi:10.1038/s41467-019-09440-2.
  * Brunk et al. [2018] E. Brunk, S. Sahoo, D. C. Zielinski, A. Altunkaya, A. Dräger, N. Mih, F. Gatto, A. Nilsson, G. A. Preciat Gonzalez, M. K. Aurich, A. Prlić, A. Sastry, A. D. Danielsdottir, A. Heinken, A. Noronha, P. W. Rose, S. K. Burley, R. M. T. Fleming, J. Nielsen, I. Thiele, B. O. Palsson,  Recon3D enables a three-dimensional view of gene variation in human metabolism,  Nature Biotechnology 36 (2018) 272–281. doi:10.1038/nbt.4072.
  * Gawthrop et al. [2021] P. J. Gawthrop, M. Pan, E. J. Crampin,  Modular dynamic biomolecular modelling: The unification of stoichiometry, thermodynamics, kinetics and data,  bioRxiv (2021). doi:10.1101/2021.03.24.436792.
  * Mason and Covert [2019] J. C. Mason, M. W. Covert,  An energetic reformulation of kinetic rate laws enables scalable parameter estimation for biochemical networks,  Journal of Theoretical Biology 461 (2019) 145–156. doi:10.1016/j.jtbi.2018.10.041.
  * Joshi et al. [2021] R. P. Joshi, A. McNaughton, D. G. Thomas, C. S. Henry, S. R. Canon, L. A. McCue, N. Kumar,  Quantum mechanical methods predict accurate thermodynamics of biochemical reactions,  ACS Omega 6 (2021) 9948–9959. doi:10.1021/acsomega.1c00997,
* Prediction of thermodynamics parameters from first-principles.
  * Chang et al. [2020] A. Chang, L. Jeske, S. Ulbrich, J. Hofmann, J. Koblitz, I. Schomburg, M. Neumann-Schaal, D. Jahn, D. Schomburg,  BRENDA, the ELIXIR core data resource in 2021: new developments and updates,  Nucleic Acids Research 49 (2020) D498–D508. doi:10.1093/nar/gkaa1025.
  * Wittig et al. [2017] U. Wittig, M. Rey, A. Weidemann, R. Kania, W. Müller,  SABIO-RK: an updated resource for manually curated biochemical reaction kinetics,  Nucleic Acids Research 46 (2017) D656–D660. doi:10.1093/nar/gkx1065.
  * Glont et al. [2020] M. Glont, C. Arankalle, K. Tiwari, T. V. N. Nguyen, H. Hermjakob, R. S. Malik-Sheriff,  BioModels Parameters: a treasure trove of parameter values from published systems biology models,  Bioinformatics 36 (2020) 4649–4654. doi:10.1093/bioinformatics/btaa560.
  * Roth et al. [2020] Y. D. Roth, Z. Lian, S. Pochiraju, B. Shaikh, J. R. Karr,  Datanator: an integrated database of molecular data for quantitatively modeling cellular behavior,  Nucleic Acids Research 49 (2020) D516–D522. doi:10.1093/nar/gkaa1008,
* Web-based tool to retrieve various types of molecular data.
  * Fröhlich et al. [2018] F. Fröhlich, T. Kessler, D. Weindl, A. Shadrin, L. Schmiester, H. Hache, A. Muradyan, M. Schütte, J.-H. Lim, M. Heinig, F. J. Theis, H. Lehrach, C. Wierling, B. Lange, J. Hasenauer,  Efficient parameter estimation enables the prediction of drug response using a mechanistic pan-cancer pathway model,  Cell Systems 7 (2018) 567–579.e6. doi:10.1016/j.cels.2018.10.013,
* Software for scalable simulation and sensitivity analysis of ODE models.
  * Schmiester et al. [2019] L. Schmiester, Y. Schälte, F. Fröhlich, J. Hasenauer, D. Weindl,  Efficient parameterization of large-scale dynamic models based on relative measurements,  Bioinformatics 36 (2019) 594–602. doi:10.1093/bioinformatics/btz581.
  * Fröhlich et al. [2021] F. Fröhlich, D. Weindl, Y. Schälte, D. Pathirana, L. Paszkowski, G. T. Lines, P. Stapor, J. Hasenauer,  AMICI: high-performance sensitivity analysis for large ordinary differential equation models,  Bioinformatics (2021). doi:10.1093/bioinformatics/btab227, btab227.
  * Somogyi et al. [2015] E. T. Somogyi, J.-M. Bouteiller, J. A. Glazier, M. König, J. K. Medley, M. H. Swat, H. M. Sauro,  libRoadRunner: a high performance SBML simulation and analysis library,  Bioinformatics 31 (2015) 3315–3321. doi:10.1093/bioinformatics/btv363.
  * Mitra and Hlavacek [2019] E. D. Mitra, W. S. Hlavacek,  Parameter estimation and uncertainty quantification for systems biology models,  Current Opinion in Systems Biology 18 (2019) 9–18. doi:10.1016/j.coisb.2019.10.006.
  * Villaverde et al. [2018] A. F. Villaverde, F. Fröhlich, D. Weindl, J. Hasenauer, J. R. Banga,  Benchmarking optimization methods for parameter estimation in large kinetic models,  Bioinformatics 35 (2018) 830–838. doi:10.1093/bioinformatics/bty736.
  * Christodoulou et al. [2018] D. Christodoulou, H. Link, T. Fuhrer, K. Kochanowski, L. Gerosa, U. Sauer,  Reserve Flux Capacity in the Pentose Phosphate Pathway Enables Escherichia coli’s Rapid Response to Oxidative Stress,  Cell Syst 6 (2018) 569–578. doi:10.1016/j.cels.2018.04.009.
  * Schreiber et al. [2020] F. Schreiber, B. Sommer, T. Czauderna, M. Golebiewski, T. E. Gorochowski, M. Hucka, S. M. Keating, M. König, C. Myers, D. Nickerson, D. Waltemath,  Specifications of standards in systems and synthetic biology: status and developments in 2020,  Journal of Integrative Bioinformatics 17 (2020). doi:10.1515/jib-2020-0022.
  * Schmiester et al. [2021] L. Schmiester, Y. Schälte, F. T. Bergmann, T. Camba, E. Dudkin, J. Egert, F. Fröhlich, L. Fuhrmann, A. L. Hauber, S. Kemmer, P. Lakrisenko, C. Loos, S. Merkt, W. Müller, D. Pathirana, E. Raimúndez, L. Refisch, M. Rosenblatt, P. L. Stapor, P. Städter, D. Wang, F.-G. Wieland, J. R. Banga, J. Timmer, A. F. Villaverde, S. Sahle, C. Kreutz, J. Hasenauer, D. Weindl,  PEtab—interoperable specification of parameter estimation problems in systems biology,  PLOS Computational Biology 17 (2021) 1–10. doi:10.1371/journal.pcbi.1008646.
  * Porubsky et al. [2020] V. L. Porubsky, A. P. Goldberg, A. K. Rampadarath, D. P. Nickerson, J. R. Karr, H. M. Sauro,  Best practices for making reproducible biochemical models,  Cell Systems 11 (2020) 109–120. doi:10.1016/j.cels.2020.06.012.
  * Villaverde et al. [2019] A. F. Villaverde, E. Raimúndez-Álvarez, J. Hasenauer, J. R. Banga,  A comparison of methods for quantifying prediction uncertainty in systems biology,  IFAC-PapersOnLine (2019). doi:10.1016/j.ifacol.2019.12.234.
  * Wieland et al. [2021] F.-G. Wieland, A. L. Hauber, M. Rosenblatt, C. Tönsing, J. Timmer,  On structural and practical identifiability,  Current Opinion in Systems Biology 25 (2021) 60–69. doi:10.1016/j.coisb.2021.03.005,
* Methods for identifiability and uncertainty analysis.
  * Roy and Finley [2019] M. Roy, S. D. Finley,  Metabolic reprogramming dynamics in tumor spheroids: Insights from a multicellular, multiscale model,  PLOS Computational Biology 15 (2019) 1–36. doi:10.1371/journal.pcbi.1007053.


# Scalable Identification of Partially Observed Systems with Certainty-
Equivalent EM

Kunal Menda  Jean de Becdelièvre  Jayesh K. Gupta  Ilan Kroo  Mykel J.
Kochenderfer  Zachary Manchester

###### Abstract

System identification is a key step for model-based control, estimator design,
and output prediction. This work considers the offline identification of
partially observed nonlinear systems. We empirically show that the certainty-
equivalent approximation to expectation-maximization can be a reliable and
scalable approach for high-dimensional deterministic systems, which are common
in robotics. We formulate certainty-equivalent expectation-maximization as
block coordinate-ascent, and provide an efficient implementation. The
algorithm is tested on a simulated system of coupled Lorenz attractors,
demonstrating its ability to identify high-dimensional systems that can be
intractable for particle-based approaches. Our approach is also used to
identify the dynamics of an aerobatic helicopter. By augmenting the state with
unobserved fluid states, a model is learned that predicts the acceleration of
the helicopter better than state-of-the-art approaches. The codebase for this
work is available at https://github.com/sisl/CEEM.

System Identification, Expectation-Maximization, Robotics

[table]capposition=top capbtabboxtable[][0.85]

##  1 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 addresses 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 systems (shown in
Figure 1) that are high-dimensional and non-linear. We consider situations in
which the system’s state cannot be inferred from a single observation, but
instead requires inference over time-series of observations.

$x_{t}$$x_{t+1}$$u_{t}$$y_{t}$$y_{t+1}$$\cdots$$\cdots$$\cdots$$\hat{x}_{t}$$\hat{x}_{t+1}$$f_{\theta}(\hat{x},u,t)$$g_{\theta}(\hat{x},u,t)$$g_{\theta}(\hat{x},u,t+1)$$\cdots$$\cdots$SystemIdentification

Figure 1: A graphical model representing a partially observed dynamical
system. Gray-box identification algorithms attempt to search a model class of
dynamics and observation models for the model that maximizes the likelihood of
the observations.

The problem of identifying systems from partial observations arises in
robotics (Punjani & Abbeel, 2015; Cory & Tedrake, 2008; Ordonez et al., 2017)
as well as domains such as chemistry (Gustavsson, 1975) and biology (Sun et
al., 2008). In many robotic settings, we have direct measurements of a robot’s
pose and velocity, but in many cases we cannot directly observe relevant
quantities such as the temperature of actuators, the state the environment
around the robot, or the intentions of other agents. For example, Abbeel et
al. (2010) attempted to map the pose and velocity of an aerobatic helicopter
to its acceleration. They found their model to be inaccurate when predicting
aggressive maneuvers because of the substantial airflow generated by the
helicopter that affected the dynamics. Since it is often impossible to
directly measure the state of the airflow around a vehicle, identification
must be with only partial observability.

System identification is a mature field with a rich history (Ljung, 1999,
2010). A variety of techniques have been proposed to learn predictive models
from time series data. Autoregressive approaches directly map a time-history
of past inputs to observations, without explicitly reasoning about unobserved
states (Billings, 2013), and are the state-of-the-art approach to the
aforementioned problem of modeling the aerobatic helicopter (Punjani & Abbeel,
2015). In contrast, state-space models (SSM) assume an unobserved state
$x_{t}$ that evolves over time and emits observations $y_{t}$ that we measure,
as shown in Figure 1. Recurrent Neural Networks (RNNs) (Bailer-Jones et al.,
1998; Zimmermann & Neuneier, 2000) are a form of black-box non-linear SSM that
can be fit to observation and input time-series, and Subspace Identification
(SID) methods (Van Overschee & De Moor, 1994) can be used to fit linear SSMs.

However, in many cases prior knowledge can be used to specify structured,
parametric models of the system (Gupta et al., 2019, 2020) in state-space
form, commonly refered to as gray-box models. Such models can be trained with
less data and used with a wider array of control and state-estimation
techniques than black-box models (Gupta et al., 2019, 2020; Lutter et al.,
2019b, a).

To identify partially observed gray-box models, the unobserved state-
trajectory is often considered as a missing data, and techniques based on
Expectation-Maximization (EM) are used (Dempster et al., 1977; Schön et al.,
2011; Kantas et al., 2015; Ghahramani & Roweis, 1999). The smoothing step
(E-step) deals with state inference—the current system dynamics estimate is
used to infer the distribution of unobserved state-trajectories conditioned on
the observations $p(x_{1:T}\mid y_{1:T})$. This distribution is sometimes
called joint smoothing distribution in the literature, and it is used to
estimate the expected log-likelihood of the observations. In the learning step
(M-step), the system’s dynamics estimate is updated such that it maximizes the
expected log-likelihood. The smoothing step can typically be approached with
particle approximations of $p(x_{1:T}\mid y_{1:T})$, but naive implementations
of particle smoothers can be computationally intensive and become rapidly
intractable in high dimensions. Across various fields and disciplines,
numerous methods have been developed to alleviate the computational burden,
several of which are discussed in Section 2.3.

This work is motivated by robotics applications, in which the following
assumptions are often valid:

  * •

Systems evolve nearly deterministically, which implies that the process noise
is small and unimodal, and,

  * •

The distribution of states conditioned on the observations $p(x_{1:T}\mid
y_{1:T})$ is unimodal.

We study the benefits of making the certainty-equivalent approximation in the
E-step of the EM procedure, and we refer to this approach as CE-EM.
Specifically, we use non-linear programming to tractably find the maximum-
likelihood (ML) point-estimate of the unobserved states, and use it in lieu of
a fully characterized approximate distribution. The contributions of this
paper are to describe an efficient implementation of CE-EM, and to test the
approximation against state-of-the-art approaches on a variety of system-
identification problems. We demonstrate on a system of Lorenz attractors that:

  * •

CE-EM can be faster and more reliable than approaches using particle
approximations,

  * •

CE-EM scales to high-dimensional problems, and,

  * •

CE-EM learns unbiased parameter estimates on deterministic systems with
unimodal $p(x_{1:T}\mid y_{1:T})$.

We also demonstrate the algorithm on the problem of identifying the dynamics
of an aerobatic helicopter. We show that a non-linear SSM can be trained with
CE-EM that outperforms various approaches including the most recent work done
on this dataset (Punjani & Abbeel, 2015). A codebase implementing CE-EM and
other supplementary material can be found at our website:
https://sites.google.com/stanford.edu/ceem/.

##  2 Background

This section states the nonlinear system identification problem with partial
observations and discusses approaches that use EM to solve it.

###  2.1 Formal Problem Statement

In this work, we assume that we are given a batch of trajectories containing
observations $y_{1:T}\in\mathbb{R}^{m\times T}$ of a dynamical system as it
evolves over a time horizon $T$, possibly forced by some known input sequence
$u_{1:T}$. We assume that this dynamical system has a state
$x\in\mathbb{R}^{n}$ that evolves and generates observations according to the
following equations,

$\displaystyle x_{t+1}$ $\displaystyle=f(x_{t},u_{t},t)+w_{t},$ $\displaystyle
w_{t}\sim p_{w}(\cdot)$ (1) $\displaystyle y_{t}$
$\displaystyle=g(x_{t},u_{t},t)+v_{t},$ $\displaystyle v_{t}\sim p_{v}(\cdot)$

where $w_{t}$ is referred to as the process noise and $v_{t}$ as the
observation noise. Both $w_{t}$ and $v_{t}$ are assumed to be additive for
notational simplicity, but this is not a required assumption.111In order to
substitute $w_{t}$ or $v_{t}$ in the arguments of $p_{v}$ and $p_{w}$, the
only requirement is to be able to express them as functions of the other
terms. Hence, the formulation is amenable to any method of injecting noise
into the system, so long as we can estimate the probability of the noise given
the other terms in their respective equations. Without loss of generality, we
can drop the dependence on $u_{t}$, absorbing it into the dependence on $t$.

We further assume that we are provided a class of parameterized models
$f_{\theta}(x,t)$ and $g_{\theta}(x,t)$ for $\theta\in\Theta$ that approximate
the dynamical system’s evolution and observation processes. The goal of our
algorithm is to find the parameters $\theta$ that maximize the likelihood of
the observations. That is, we seek to find:

$\displaystyle\theta_{\text{ML}}$
$\displaystyle=\underset{\theta}{\operatorname*{arg\,max}}~{}p(y_{1:T}\mid\theta)$
(2) $\displaystyle=\underset{\theta}{\operatorname*{arg\,max}}~{}\int
p(y_{1:T},x_{1:T}\mid\theta)dx_{1:T}$
$\displaystyle=\underset{\theta}{\operatorname*{arg\,max}}~{}\int
p(y_{1:T}\mid x_{1:T},\theta)p(x_{1:T}\mid\theta)dx_{1:T}$

Using the graphical model shown in Figure 1, we can rewrite this integral as:

$\displaystyle\theta_{\text{ML}}=\underset{\theta}{\operatorname*{arg\,max}}~{}\int$
$\displaystyle\prod_{t=1}^{T}p(y_{t}\mid x_{t},\theta)\cdot$ (3)
$\displaystyle\prod_{t=1}^{T-1}p(x_{t+1}\mid
x_{t},\theta)p(x_{1}\mid\theta)dx_{1:T}$

Finally, using the notation chosen to describe a dynamical system in Equation
1, we obtain:

$\displaystyle\theta_{\text{ML}}$
$\displaystyle=\underset{\theta}{\operatorname*{arg\,max}}\,\int~{}\left(\prod_{t=1}^{T}p_{v}(y_{t}-g_{\theta}(x_{t},t))\right)\cdot$
(4) $\displaystyle\hskip
40.0pt\left(p(x_{1})\prod_{t=1}^{T-1}p_{w}(x_{t+1}-f_{\theta}(x_{t},t))\right)dx_{1:T}$

The Expectation-Maximization (EM) algorithm has been used in the literature to
solve problems of this form.

###  2.2 Expectation-Maximization for SSM Identification

The EM algorithm is a two-step procedure that copes with the missing state
information by forming an approximation $Q(\theta,\theta_{k})$ of the joint
state and observation likelihood $p_{\theta}(x_{1:T},y_{1:T})$ at the $k$th
iteration of the algorithm:

$\displaystyle Q(\theta,\theta_{k})=\int\log{}$ $\displaystyle
p(x_{1:T},y_{1:T}\mid\theta)\cdot$ (5) $\displaystyle p(x_{1:T}\mid
y_{1:T},\theta_{k})dx_{1:T}$

Classical proofs show that maximizing $Q(\theta,\theta^{k})$ with respect to
$\theta$ results in increasing $p_{\theta}(x_{1:T},y_{1:T})$ (Dempster et al.,
1977). The EM algorithm iteratively performs two steps:

  1. 1.

Compute $Q(\theta,\theta_{k})$

  2. 2.

Update
$\theta_{k+1}=\underset{\theta}{\operatorname*{arg\,max}}\;Q(\theta,\theta_{k})$

The evaluation of $Q(\theta,\theta_{k})$ in the E-step requires the estimation
the $p(x_{1:T}\mid y_{1:T},\theta_{k})$, as well as the integration of the
likelihood over this distribution. In the linear Gaussian case, EM can be
performed exactly using Kalman smoothing (Rauch et al., 1965; Ghahramani &
Hinton, 1996). In the more general case, there is no analytic solution, and
previous work on approximating the E-step is summarized in the next section.

###  2.3 Related Work

There are general approaches that are based on particle representations of the
joint states-observations likelihood which use Sequential Monte-Carlo (SMC)
methods such as Particle Smoothing (PS) (Schön et al., 2011; Kantas et al.,
2015). With enough particles, PS can handle any system and any joint
distribution of states and observations. However, PS suffers from the curse of
dimensionality, requiring an intractably large number of particles if the
state space is high-dimensional (Snyder et al., 2008; Kantas et al., 2015). In
the simplest form, both the E-step and the M-step can be quadratic in
complexity with respect to the number of particles (Schön et al., 2011). An
important body of work has attempted to alleviate some of this burden by using
Nested SMC (Naesseth et al., 2015) and forward filtering-backward simulation
(FFBSi), (Lindsten & Schön, 2013), and conditional particle filtering
(Lindsten, 2013) in the E-step. PS can also require variance reduction
techniques, such as full-adaptation, to perform reliably, and likelihood
estimates in the E-step can be noisy (Svensson et al., 2018). Stochastic
Approximation EM can be used to stabilize learning in the presence of this
noise (Delyon et al., 1999; Svensson & Lindsten, 2018).

Another group of methods is based on linearizing the dynamics around the
current state estimates to obtain a time-dependant linear Gaussian dynamical
system, to which we can apply Kalman smoothing techniques. Ghahramani & Roweis
(1999) use Extended Kalman Smoothing (EKS) for a fast E-step and show that,
when using radial basis functions to describe the system dynamics, no
assumption is required in the M-step. The drawbacks of this approach are that
the number of radial basis functions required to accurately represent a
function grows exponentially with the input dimension, that the user cannot
specify a gray-box SSM, and that EKS performance can vary dramatically with
the hyperparameters. Goodwin & Aguero (2005) proposed a method called MAP-EM,
which linearizes around the state-trajectory that maximizes joint state and
observation likelihood, approximates process and observation noise as
Gaussian, and performs a local version of EM using this linearization. In
comparison, the method we consider assumes that the maximum likelihood state-
trajectory estimate concentrates all of the probability mass. Goodwin & Aguero
(2005) call this simplification Certainty Equivalent EM (CE-EM) and report
worse results than their approach on simple examples.

This paper shows that CE-EM can actually be a fast, simple, and reliable
approach. The method we use here is tailored to systems that are not dominated
by process noise, in which the state-trajectory distribution is unimodal, and
whose high dimensional state-space make other methods intractable. Such
scenarios are regularly encountered in robotics, but perhaps less so in other
applications domains of system identification such as chemistry, biology, or
finance.

##  3 Methodology

This section introduces the CE-EM approximation and presents an efficient
algorithm to identify high-dimensional robotic systems.

###  3.1 Certainty-Equivalent Expectation Maximization

Under the certainty-equivalence condition, the distribution of states
conditioned on observations is assumed to be a Dirac delta function. That is,
we assume:

$p(x_{1:T}|y_{1:T},\theta)=\delta_{x_{1:T}^{\text{ML}}(\theta)}(x_{1:T})$ (6)

where

$\displaystyle x_{1:T}^{\text{ML}}(\theta)$
$\displaystyle=\operatorname*{arg\,max}_{x_{1:T}}p(x_{1:T}\mid
y_{1:T},\theta)$ (7)
$\displaystyle=\operatorname*{arg\,max}_{x_{1:T}}p(x_{1:T},y_{1:T}\mid\theta)/p(y_{1:T}\mid\theta)$
$\displaystyle=\operatorname*{arg\,max}_{x_{1:T}}p(x_{1:T},y_{1:T}\mid\theta)$

In other words, this assumption implies that there is only one state
trajectory that satisfies the system dynamics, and is coherent with the data.

The E-step can be rewritten as:

$\displaystyle Q(\theta,\theta_{k})$
$\displaystyle=\int\log~{}p(x_{1:T},y_{1:T}\mid\theta)\cdot$ (8)
$\displaystyle\hskip
50.00008pt\delta_{x_{1:T}^{\text{ML}}(\theta_{k})}(x_{1:T})dx_{1:T}$
$\displaystyle=\log p(x^{\text{ML}}_{1:T}(\theta_{k}),y_{1:T}\mid\theta)$

By observing the maximization of the expression in Equation 7 as the E-step,
and maximizing the expression in Equation 8 as the M-step, we see that CE-EM
is simply block coordinate-ascent on the single joint-objective:

$\begin{split}J(x_{1:T},\theta)=&\log p(x_{1:T},y_{1:T}\mid\theta)\\\ =&\log
p(x_{1})+\sum_{t=1}^{T}\log p_{v}(y_{t}-g_{\theta}(x_{t},t))\\\
&+\sum_{t=1}^{T-1}\log p_{w}(x_{t+1}-f_{\theta}(x_{t},t))\end{split}$ (9)

Jointly maximizing this objective over $x_{1:T}$ and $\theta$ yields
$\theta_{\text{ML}}$.

Though this objective can be optimized as is by a non-linear optimizer, it is
not necessarily efficient to do so since $\theta$ and $x_{1:T}$ are highly
coupled, leading to inefficient and potentially unstable updates. For this
reason, we still opt for the block coordinate-ascent approach similar to EM,
where $\theta$ and $x_{1:T}$ are each sequentially held constant while the
other is optimized. By viewing EM as block coordinate-ascent, we may also
borrow good practices typically employed when running the algorithm (Shi et
al., 2016). In particular, we employ trust-region regularization in order to
guarantee convergence in even non-convex settings (Grippo & Sciandrone, 2000).

####  3.1.1 Smoothing

At iteration $k$, we first perform smoothing by holding $\theta$ constant and
finding the ML point-estimate for $x_{1:T}$ as follows:

$x^{(k)}_{1:T}=\underset{x_{1:T}}{\operatorname*{arg\,max}}~{}J(x_{1:T},\theta^{(k-1)})-\rho_{x}\lVert
x_{1:T}-x_{1:T}^{(k-1)}\rVert^{2}_{2}$ (10)

where $\rho_{x}$ scales a soft trust-region regularizer similar to damping
terms found in Levenberg-Marquardt methods (Levenberg, 1944; Marquardt, 1963).
We note that the Hessians of this objective with respect to $x_{1:T}$ are
block-sparse, and as a result this step can be efficiently solved with a
second-order optimizer in $\mathcal{O}(n^{2}T)$.

If the cost function $J(x_{1:T},\theta)$ is a sum of quadratic terms, i.e. if
all stochasticity is Gaussian, then the smoothing can be solved with a Gauss-
Newton method, where the Hessian matrix is approximated using first-order
(Jacobian) information. In this case, the solution to the smoothing step is
equivalent to iteratively performing EKS (Aravkin et al., 2017). However, EKS
typically involves forward and backward passes along each trajectory,
performing a Riccati-like update at each time step. Though square-root
versions of the algorithm exist to improve its numerical stability (Bierman,
1981), solving the problem with batch non-linear least-squares as opposed to
iteratively performing a forward-backward method can improve stability for
long time-series (Van Dooren, 1981).

Input: observations $y_{1:T}$, control inputs $u_{1:T}$, stopping criterion
$tol$

Initialize: model parameters $\theta$

Initialize: hidden state estimates $x_{1:T}$

repeat

$J^{k}\xleftarrow{}J(x^{k}_{1:T},\theta^{k})$

$x^{k+1}_{1:T}\xleftarrow{}\underset{\tilde{x}_{1:T}}{\arg\max}~{}J(\tilde{x}_{1:T},\theta^{k})-\rho_{x}\lVert\tilde{x}_{1:T}-x^{k}_{1:T}\rVert^{2}_{2}$

$\theta^{k+1}\xleftarrow{}\underset{\tilde{\theta}}{\arg\max}~{}J(x^{k+1}_{1:T},\tilde{\theta})-\rho_{\theta}\lVert\tilde{\theta}-\theta^{k}\rVert^{2}_{2}+\log
p(\tilde{\theta})$

until $J^{k}-J(x^{k+1}_{1:T},\theta^{k+1})\leq tol$

Algorithm 1 CE-EM Implementation

####  3.1.2 Learning

In the learning step, we hold $x_{1:T}$ constant and find:

$\theta^{(k)}=\underset{\theta}{\arg\max}~{}J(x^{(k)}_{1:T},\theta)-\rho_{\theta}\lVert\theta-\theta^{(k-1)}\rVert^{2}_{2}+\log
p(\theta)$ (11)

Again, $\rho_{\theta}$ scales a soft trust-region regularizer, and specifying
$\log p(\theta)$ allows us to regularize $\theta$ toward a prior. The above
optimization problem can be solved using any non-linear optimizer. We find the
Nelder-Mead (Gao & Han, 2012) scheme well-suited for small parameter spaces,
and first-order schemes such as Adam (Kingma & Ba, 2015), or quasi-second-
order schemes such as L-BFGS (Liu & Nocedal, 1989) suited for larger parameter
spaces such as those of neural networks. The routine, summarized in Algorithm
1, iterates between the smoothing and learning steps until convergence.

It should be noted that making the certainty-equivalent approximation is
generally known to bias parameter estimates (Celeux & Govaert, 1992), but can
yield the correct solutions under the assumptions we make (Neal & Hinton,
1998).

##  4 Experiments

Readers of Goodwin & Aguero (2005) will likely be left with the incorrect
impression that CE-EM is an inferior method that would perform poorly in
practice. The objective of our experiments is to demonstrate that the CE-EM
algorithm is capable of identifying high-dimensional non-linear systems in
partially observed settings. We do so in simulation by identifying the
parameters of a system of partially observed coupled Lorenz attractors, as
well as by identifying the dynamics of a real aerobatic helicopter. In the
second experiment, we build on previous analysis of the dataset (Abbeel et
al., 2010; Punjani & Abbeel, 2015) by attempting to characterize the
interaction of the helicopter with the fluid around it, without having any
direct observation of the fluid state. Instructions for reproducing all
experiments are included in the supplementary material.

###  4.1 Identification of Coupled Lorenz Systems

In this experiment, we show that:

  1. 1.

CE-EM learns unbiased parameter estimates of systems that are close to
deterministic, and,

  2. 2.

CE-EM scales to high-dimensional problems in which particle-based methods can
be intractable.

To justify these claims, we use a system that is sufficiently non-linear and
partially observable to make particle-based smoothing methods intractable. We
choose a system of coupled Lorenz attractors for this purpose, owing to their
ability to exhibit chaotic behavior and their use in non-linear atmospheric
and fluid flow models (Bergé et al., 1984). Arbitrary increases in state
dimensionality can be achieved by coupling multiple individual attractors. The
state of a system with $K$ coupled Lorenz attractors is
$x\in\mathbb{R}^{3K}=\\{\ldots,x_{1,k},x_{2,k},x_{3,k},\ldots\\}$. The
dynamics of the system are as follows:

$\displaystyle\dot{\tilde{x}}_{1,k}$
$\displaystyle=\sigma_{k}(x_{2,k}-x_{1,k})$ (12)
$\displaystyle\dot{\tilde{x}}_{2,k}$
$\displaystyle=x_{1,k}(\rho_{k}-x_{3,k})-x_{2,k}$
$\displaystyle\dot{\tilde{x}}_{3,k}$
$\displaystyle=x_{1,k}x_{2,k}-\beta_{k}x_{3,k}$ $\displaystyle\dot{x}$
$\displaystyle=\dot{\tilde{x}}+Hx$

where $H$ is an $\mathbb{R}^{3K\times 3K}$ matrix.

We nominally set the parameters $(\sigma_{k},\rho_{k},\beta_{k})$ to the
values $(10,28,8/3)$, and randomly sample the entries of $H$ from a 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\in\mathbb{R}^{(3K-2)}$ is derived from $x$ as follows:

$y=Cx+v,\quad v\sim\mathcal{N}(0,\sigma_{v}^{2}\mathbf{I})$ (13)

where $C\in\mathbb{R}^{(3K-2)\times 3K}$ is a known matrix with full row-rank,
and $v$ is the observation noise sampled from a Gaussian with diagonal
covariance $\sigma_{v}^{2}\mathbf{I}$. The entries of $C$ are also randomly
sampled from a standard normal distribution. In the following experiments, we
simulate the system for $T=128$ timesteps at a sample rate of $\Delta
t=0.04$\mathrm{s}$$, and integrate the system using a $4$th-order Runge-Kutta
method. Initial conditions for each trajectory are sampled such that
$x_{1,k}\sim\mathcal{N}(-6,2.5^{2}),x_{2,k}\sim\mathcal{N}(-6,2.5^{2}),x_{3,k}\sim\mathcal{N}(24,2.5^{2})$.

####  4.1.1 Unbiased Estimation in Deterministic Settings

Table 1: Mean parameter estimates and standard errors for a single Lorenz
system simulated with various $\sigma_{w}$ and $\sigma_{v}$.

$\sigma_{w}$ | $\sigma_{v}$ | $\sigma$ | $\rho$ | $\beta$
---|---|---|---|---
$0.001$ | $0.01$ | $10.011(0.012)$ | $28.000(0.001)$ | $2.667(0.000)$
$0.010$ | $0.01$ | $10.017(0.012)$ | $28.000(0.001)$ | $2.668(0.001)$
$0.100$ | $0.01$ | $10.064(0.036)$ | $27.996(0.013)$ | 2.676 (0.004)
$0.001$ | $0.05$ | $10.006(0.016)$ | $27.998(0.002)$ | $2.666(0.001)$
$0.001$ | $0.10$ | $9.998(0.022)$ | $27.995(0.004)$ | $2.665(0.001)$

To test the conditions under which CE-EM learns unbiased parameter estimates,
we simulate a single Lorenz system with $H=\mathbf{0}$ and known
$C\in\mathbb{R}^{2\times 3}$. We introduce and vary the process noise
$w\sim\mathcal{N}(0,\sigma_{w}^{2}\mathbf{I})$, and vary the observation noise
coefficient $\sigma_{v}$, and then attempt to estimate the parameters
$(\sigma,\rho,\beta)$. Using initial guesses within $10\%$ of the system’s
true parameter values, we run CE-EM on a single sampled trajectory. For each
choice of $\sigma_{w}$ and $\sigma_{v}$, we repeat this process for $10$
random seeds.

Table 1 shows the mean and standard errors of parameter estimates for various
$\sigma_{w}$ and $\sigma_{v}$. We highlight in red the mean estimates that are
not within two standard errors of their true value. We see that $\sigma$ and
$\rho$ are estimated without bias for all scenarios. However, the estimate of
$\beta$ appears to become biased as the process noise is increased, but not as
the observation noise is increased. This supports the assumption that the
objective used in CE-EM is sound when systems evolve close to
deterministically, but can be biased if it is not.

####  4.1.2 Comparison to Particle Based Methods

In Section 2, we discussed methods for parameter estimation in state-space
systems that are based on particle-filtering and smoothing (Kantas et al.,
2015). Since in their E-step, these methods approximate the distribution over
unobserved state-trajectories as opposed to only their point estimate, such
methods can be asymptotically unbiased. However, for a finite number of
particles, such methods can result in high-variance estimates. In this
experiment, we compare the bias resulting from using CE-EM with the variance
of using a state-of-the-art Particle EM algorithm.

We attempt to identify the parameters $(\sigma,\rho,\beta)$ of the same single
Lorenz system as in Section 4.1.1. However, we introduce process noise
$w\sim\mathcal{N}(0,0.1^{2}\mathbf{I})$ and observation noise
$v\sim\mathcal{N}(0,0.5^{2}\mathbf{I})$. We use a training dataset of four
trajectories sampled with conditions specified in Section 4.1.

The performance of Particle EM can vary substantially depending on
implementation of the particle filter and smoother in the E-step, and
implementation of the M-step. We use a fully-adapted particle filter with
systematic and adaptive resampling, following the recommendations of Doucet et
al. (2000); Hol et al. (2006). Furthermore, we use the FFBSi algorithm
(Godsill et al., 2004; Lindsten & Schön, 2013) in order to generate iid
samples of smoothed state-trajectories, while avoiding complexity that is
quadratic in the number of particles experienced by the forward filter-
backward smoother (FFBSm) approach (Doucet et al., 2000). We then use
Stochastic Approximation EM (SAEM) (Delyon et al., 1999) to perform the
M-step.222We have found that using FFBSi with SAEM performs more reliably than
FFBSm with the M-step recommended by (Schön et al., 2011), and both
implementations can be found in the associated codebase. We use $N_{p}=100$
particles for filtering, and sample $N_{s}=10$ smoothed trajectories using
FFBSi.

Figure 2 shows the estimated parameters versus EM epoch using CE-EM and
Particle EM. We plot learning curves for 10 random seeds, each of which
initializes parameter estimates to within 10% of their true value, and uses a
different set of training trajectories. We see that CE-EM consistently
converges to accurate parameter estimates in approximately 5 epochs. Estimates
of Particle EM appear to initially diverge but in all but one case converge to
a similar accuracy in 50 epochs. Furthermore, since the complexity of FFBSi is
$O(N_{p}N_{s})$,333Variants of FFBSi algorithm improve on the $O(N_{p}N_{s})$
complexity (Lindsten & Schön, 2013). the runtime per epoch of Particle EM is
$47\times$ more than that of CE-EM.444Experiments were run on a computer with
an Intel® Core™ i7-6700K CPU @ 4.00GHz $\times 8$ processor and 32GB of
memory. Operations are parallelized across trajectories for both
implementations. Since the variance of particle-based methods generally
increases with the effective dimension of the system, the bias induced by CE-
EM may be a worthwhile trade-off for fast and reliable parameter estimation.

Figure 2: Comparison of CE-EM and Particle EM on parameter estimation of a
single Lorenz system.

####  4.1.3 Convergence of CE-EM on High-Dimensional Problems

To demonstrate that CE-EM 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, we test whether CE-EM converges to more
accurate estimates when more trajectories are provided. To test these claims,
we sample $2$, $4$, and $8$ trajectories from a deterministic system with
parameters $\theta_{\text{true}}$, and $\sigma_{v}=0.01$. We randomly
initialize each element of the parameters being optimized
($\theta=[\sigma_{1:K},\rho_{1:K},\beta_{1:K},H]$) to within 10% of the their
value in $\theta_{\text{true}}$. We then run CE-EM on each batch, tracking the
error in the estimated dynamics as training proceeds. We measure this error,
which we call $\epsilon(\theta)$, as follows:

$\epsilon(\theta)=\mathbf{E}_{x\sim p(x_{0})}\left[\lVert
f_{\theta}(x)-f_{\theta_{\text{true}}}(x)\rVert_{2}\right]$ (14)

In the learning step, we do not regularize $\theta$ to a prior and set
$\rho_{\theta}=0$.

For comparison, we also run the same Particle EM implementation as in the
previous experiment on this problem. We use the same $N_{p},N_{s}$, and
hyperparameters as before.

Figure 3: Error in estimated dynamics as CE-EM trains on a varied number of
trajectories from a system of six coupled Lorenz attractors.

Figure 3 shows the results of this experiment for four random seeds for each
batch size. We can see that, as the number of trajectories used in training
increases, the error in the estimated dynamics tends toward zero. Furthermore,
we see that CE-EM convergences monotonically to a local optimum in all cases.
In contrast, Particle EM appears to initially improve but then converges to
very poor parameter estimates.

The experiments conducted thus far have demonstrated that CE-EM can learn
unbiased parameter estimates of nearly-deterministic systems, and can scale to
high-dimensional problems for which particle-based methods are intractable. In
the next experiment, we use CE-EM to characterize the effect of unobserved
states on the dynamics of an aerobatic helicopter.

###  4.2 Characterizing Aerobatic Helicopter Dynamics

Characterizing the dynamics of a helicopter undergoing aggressive aerobatic
maneuvers is widely considered to be a challenging system-identification
problem (Abbeel et al., 2010; Punjani & Abbeel, 2015). The primary challenge
is 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 is 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 a system’s state can be reconstructed
with a finite number of lagged-observations of it (Takens, 1981). Instead of
attempting to estimate the unobserved fluid state, they directly learn a
mapping from a $0.5\,$\mathrm{s}$$ history of observed state measurements and
control commands to the forces acting on the helicopter.

This approach is equivalent to considering the past $0.5\,$\mathrm{s}$$ of
observations as the system’s state. However, it can require a very large
number of lagged observations to represent complex phenomena. In reality, the
characteristic time of unsteady flows around helicopters can easily be up to
tens of seconds. 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.

(a)

(b)

Figure 4: (a): Test performance of optimized models on various trajectories.
Error bars on SID represent the standard deviation of performance of the 10
trained models. (b): Predicted acceleration along axis $x$, $y$ and $z$ in the
body frame for one of the harder test set trajectories. Larger versions of
these plots, including rotational accelerations, can be found in Section A.3.

####  4.2.1 Objective and dataset

The objective of this learning problem is to predict $y_{t}$, the helicopter’s
acceleration at time $t$, from an input vector $u_{t}$ containing the current
measured state of the helicopter (its velocity and rotation rates) and the
control commands.

We use data collected by the Stanford Autonomous Helicopter Project (Abbeel et
al., 2010). Trajectories are split into $10\,$\mathrm{s}$$ long chunks and
then randomly distributed into train, test, and validation sets according to
the established protocol (Abbeel et al., 2010; Punjani & Abbeel, 2015) and
summarized in Section A.1. The train, test, and validation sets respectively
contain 466, 100, and 101 trajectories of 500 time-steps each.

A simple success metric on a given trajectory is the root mean squared
prediction error,

$\text{RMSE}=\sqrt{\frac{1}{T}\sum_{t=1}^{T}\left\lVert
y_{t}^{\text{(measured)}}-y_{t}^{\text{(pred)}}\right\rVert^{2}_{2}}$ (15)

where $y_{t}^{\text{(measured)}}$ is the measured force from the dataset,
$y_{t}^{\text{(pred)}}$ is the force predicted by the model, and $T$ is the
number of time-steps in each trajectory.

####  4.2.2 Previous work and baselines

Naive: We first consider a 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_{t}=\text{NN}_{\theta_{n}}(u_{t})$, where $\text{NN}_{\theta_{n}}$ is a
neural-network with parameters $\theta_{n}$. We refer to this model as the
naive model.

H25: We also compare to the work of Punjani & Abbeel (2015). They predict
$y_{t}$ using a time-history $u_{t-H:t}$ of $H=25$ lagged observations of the
helicopter’s measured state and control commands. This input is passed through
a ReLU-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 $y_{t}=\text{NN}_{\theta_{h}}(u_{t-H:t})$
with parameters $\theta_{h}$. We call this neural network model the H25 model.
Both of these models can be trained via stochastic gradient descent to
minimize the Mean-Squared-Error (MSE) of their predictions for $y$. The
optimization methodology for these models is described in Section A.2.

SID: As a third baseline, we compare with subspace-identification methods (Van
Overschee & De Moor, 1994). We let
$\tilde{y}_{t}=y_{t}-\text{NN}_{\theta_{n}}(u_{t})$ be the prediction errors
of the trained naive model. We use the MATLAB command n4sid to fit a linear
dynamical system of the following form:

$\displaystyle x_{t+1}$ $\displaystyle=A_{\text{s}}x_{t}+B_{\text{s}}u_{t};$
$\displaystyle\quad\tilde{y}_{t}$
$\displaystyle=C_{\text{s}}x_{t}+D_{\text{s}}u_{t}$ (16)

Here, $x\in R^{d}$ is the unobserved state with arbitrary dimension $d$. The
learned parameters are
$\theta_{\text{s}}=[A_{\text{s}},B_{\text{s}},C_{\text{s}},D_{\text{s}}]$. We
use a state dimension of $10$ and call this model the SID model. The n4sid
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. This approach fits a linear system
to the residual error of the naive model, therefore the prediction of $y_{t}$
obtained from the SID model is a nonlinear function of the states $x_{t}$.

LSTM: We also train an LSTM (Hochreiter & Schmidhuber, 1997) on the residual
time-series $\tilde{y}_{1:T}$ and $u_{1:T}$.

####  4.2.3 Non-linear Unobserved State Model

We use CE-EM to train a non-linear SSM. Similar to the parameterization used
for subpace-identification, we fit the prediction errors of the naive model
using the following dynamical system:

$\displaystyle x_{t+1}$ $\displaystyle=A_{\text{\tiny NL}}x_{t}+B_{\text{\tiny
NL}}u_{t}$ (17) $\displaystyle\tilde{y}_{t}$ $\displaystyle=C_{\text{\tiny
NL}}x_{t}+D_{\text{\tiny NL}}u_{t}+\text{NN}_{\bar{\theta}_{\text{\tiny
NL}}}(x_{t},u_{t})$

where $\text{NN}_{\bar{\theta}_{\text{\tiny NL}}}$ is a neural network, and
$\theta_{\text{\tiny NL}}=[A_{\text{\tiny NL}},B_{\text{\tiny
NL}},C_{\text{\tiny NL}},D_{\text{\tiny NL}},\bar{\theta}_{\text{\tiny NL}}]$
are the learned parameters. We introduce non-linearity only in the observation
function because it is known from Koopman theory that a non-linear system can
be approximated by a high-dimensional linear system provided the correct non-
linear mapping between them (Brunton et al., 2016).

While learning, we assume that both process and observation noise are
distributed with diagonal Gaussian covariance matrices $\sigma_{w}\mathbf{I}$
and $\sigma_{v}\mathbf{I}$ respectively. The values of $\sigma_{w}$ and
$\sigma_{v}$ are treated as hyperparmeters of CE-EM, and are both set to 1.
Here as well, we use a state dimension of $10$ and call this model the NL
model. The optimization methodology for this model is described in Section
A.2.

It should be noted that the system we learn need not actually correspond to 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 a gray-
box model trained with CE-EM.

####  4.2.4 Evaluation methodology

The test RMSE of the naive, H25, and LSTM models can be evaluated directly on
the test trajectories using next-step prediction. However, the SID and NL
models require an estimate of the unobserved state 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 a given time-step, a prediction of $\tilde{y}_{t}$ is made
using the current estimate of $x_{t}$, and is used in the computation of RMSE.
The state-estimate is then corrected with the measured $\tilde{y}_{t}$.

####  4.2.5 Results

Figure 4(a) shows the RMSE of the compared models on trajectories in the test-
set. We see that the NL model is able to 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 Section A.3. The LSTM model also
performs poorly, on average 1.4 times worse than the H25 model. The SID model
notably outperforms the state-of-the-art 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 4(b) depicts the errors in prediction over a sample trajectory in the
test-set. Here, we also see that the NL model is able to attenuate the time-
varying error present in predictions made by the H25, suggesting that it has
accurately characterized the dynamics of unobserved, time-varying states.

This experiment validates the effectiveness of CE-EM to identify a non-linear
dynamical model of unobserved states that affect the forces acting an
aerobatic helicopter.

##  5 Conclusions

This paper presented an 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 a partially observed system of coupled Lorenz
attractors, we showed that CE-EM can perform identification on a problem that
particle-based EM methods are ill-suited for. However, we also find that CE-EM
yields biased parameter estimates in the presence of large process noise. This
bias can be partially mitigated by locally approximating the posterior state-
marginals as Gaussian, as is done by MAP-EM (Goodwin & Aguero, 2005). We then
used the algorithm to model the time-varying hidden-states that affect the
dynamics of an aerobatic helicopter. The model trained with CE-EM outperforms
state-of-the-art methods because it is able to fit large non-linear models to
unobserved states.

Numerous system-identification problems can be studied using CE-EM. Recently,
there have been tremendous efforts to characterize predictive models for the
spread of COVID-19 (Fanelli & Piazza, 2020). Limited capacity for testing the
prevalence of the disease makes relevant states partially observed, and thus
CE-EM may be useful for its modeling. We also hope to apply CE-EM to very
high-dimensional systems with sparsely coupled dynamics using general-form
consensus optimization (Boyd et al., 2011).

## Acknowledgments

This work is supported in part by DARPA under agreement number D17AP00032, and
by the King Abdulaziz City for Science and Technology (KACST) through the
Center of Excellence in Aeronautics and Astronautics. The content is solely
the responsibility of the authors and does not necessarily represent the
official views of DARPA or KACST.

## References

  * Abbeel et al. (2010) Abbeel, P., Coates, A., and Ng, A. Y.  Autonomous helicopter aerobatics through apprenticeship learning.  _The International Journal of Robotics Research_ , 29(13):1608–1639, 2010.
  * Aravkin et al. (2017) Aravkin, A., Burke, J. V., Ljung, L., Lozano, A., and Pillonetto, G.  Generalized Kalman smoothing: Modeling and algorithms.  _Automatica_ , 86:63–86, 2017.
  * Bailer-Jones et al. (1998) Bailer-Jones, C. A., MacKay, D. J., and Withers, P. J.  A recurrent neural network for modelling dynamical systems.  _Network: Computation in Neural Systems_ , 9(4):531–547, 1998.
  * Bergé et al. (1984) Bergé, P., Pomeau, Y., and Vidal, C.  _Order within Chaos: Towards a Deterministic Approach to Turbulence_.  Wiley & Sons, 1984.
  * Bierman (1981) Bierman, G. J.  A new computationally efficient fixed-interval, discrete-time smoother.  In _IEEE Conference on Decision and Control (CDC)_ , pp. 1054–1060, 1981.
  * Billings (2013) Billings, S. A.  _Nonlinear System Identification: NARMAX Methods in the Time, Frequency, and Spatio-Temporal Domains_.  Wiley, 2013.
  * Boyd et al. (2011) Boyd, S., Parikh, N., Chu, E., Peleato, B., and Eckstein, J.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  _Foundations and Trends in Machine Learning_ , 3(1):1–122, 2011.
  * Brunton et al. (2016) Brunton, S. L., Brunton, B. W., Proctor, J. L., and Kutz, J. N.  Koopman invariant subspaces and finite linear representations of nonlinear dynamical systems for control.  _PloS one_ , 11(2), 2016.
  * Celeux & Govaert (1992) Celeux, G. and Govaert, G.  A classification EM algorithm for clustering and two stochastic versions.  _Computational Statistics & Data Analysis_, 14(3):315–332, 1992.
  * Cory & Tedrake (2008) Cory, R. and Tedrake, R.  Experiments in fixed-wing UAV perching.  In _AIAA Guidance, Navigation and Control Conference_ , 2008.
  * Delyon et al. (1999) Delyon, B., Lavielle, M., Moulines, E., et al.  Convergence of a stochastic approximation version of the EM algorithm.  _The Annals of Statistics_ , 27(1):94–128, 1999\.
  * Dempster et al. (1977) Dempster, A. P., Laird, N. M., and Rubin, D. B.  Maximum likelihood from incomplete data via the EM algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 39(1):1–22, 1977.
  * Doucet et al. (2000) Doucet, A., Godsill, S., and Andrieu, C.  On sequential Monte Carlo sampling methods for Bayesian filtering.  _Statistics and Computing_ , 10(3):197–208, 2000\.
  * Fanelli & Piazza (2020) Fanelli, D. and Piazza, F.  Analysis and forecast of COVID-19 spreading in China, Italy and France.  _Chaos, Solitons & Fractals_, 134:109761, 2020.
  * Gao & Han (2012) Gao, F. and Han, L.  Implementing the Nelder-Mead simplex algorithm with adaptive parameters.  _Computational Optimization and Applications_ , 51(1):259–277, 2012.  doi: 10.1007/s10589-010-9329-3.
  * Ghahramani & Hinton (1996) Ghahramani, Z. and Hinton, G. E.  Parameter estimation for linear dynamical systems.  Technical report, CRG-TR-96-2, University of Totronto, Dept. of Computer Science, 1996.
  * Ghahramani & Roweis (1999) Ghahramani, Z. and Roweis, S. T.  Learning nonlinear dynamical systems using an EM algorithm.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pp. 431–437, 1999.
  * Godsill et al. (2004) Godsill, S. J., Doucet, A., and West, M.  Monte Carlo smoothing for nonlinear time series.  _Journal of the American Statistical Association_ , 99(465):156–168, 2004.
  * Goodwin & Aguero (2005) Goodwin, G. C. and Aguero, J. C.  Approximate em algorithms for parameter and state estimation in nonlinear stochastic models.  In _IEEE Conference on Decision and Control (CDC)_ , pp. 368–373, 2005.
  * Grippo & Sciandrone (2000) Grippo, L. and Sciandrone, M.  On the convergence of the block nonlinear Gauss–Seidel method under convex constraints.  _Operations research letters_ , 26(3):127–136, 2000.
  * Gupta et al. (2019) Gupta, J. K., Menda, K., Manchester, Z., and Kochenderfer, M. J.  A general framework for structured learning of mechanical systems.  _arXiv preprint arXiv:1902.08705_ , 2019.
  * Gupta et al. (2020) Gupta, J. K., Menda, K., Manchester, Z., and Kochenderfer, M. J.  Structured mechanical models for robot learning and control.  In _Conference on Learning for Dynamics and Control (L4DC)_ , 2020\.
  * Gustavsson (1975) Gustavsson, I.  Survey of applications of identification in chemical and physical processes.  _Automatica_ , 11(1):3–24, 1975.
  * Hochreiter & Schmidhuber (1997) Hochreiter, S. and Schmidhuber, J.  LSTM can solve hard long time lag problems.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pp. 473–479, 1997.
  * Hol et al. (2006) Hol, J. D., Schon, T. B., and Gustafsson, F.  On resampling algorithms for particle filters.  In _IEEE Nonlinear Statistical Signal Processing Workshop_ , pp. 79–82, 2006.
  * Hou & Wang (2013) Hou, Z.-S. and Wang, Z.  From model-based control to data-driven control: Survey, classification and perspective.  _Information Sciences_ , 235:3–35, 2013.
  * Jones et al. (2001–) Jones, E., Oliphant, T., Peterson, P., et al.  SciPy: Open source scientific tools for Python, 2001–.  URL http://www.scipy.org/.
  * Kantas et al. (2015) Kantas, N., Doucet, A., Singh, S. S., Maciejowski, J., and Chopin, N.  On particle methods for parameter estimation in state-space models.  _Statistical Science_ , 30(3):328–351, 2015.
  * Kingma & Ba (2015) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  In _International Conference on Learning Representations_ , 2015.
  * Korda & Mezić (2018) Korda, M. and Mezić, I.  Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control.  _Automatica_ , 93:149–160, July 2018.  doi: 10.1016/j.automatica.2018.03.046.
  * Levenberg (1944) Levenberg, K.  A method for the solution of certain non-linear problems in least squares.  _Quarterly of Applied Mathematics_ , 2(2):164–168, July 1944.  doi: 10.1090/qam/10666.
  * Lindsten (2013) Lindsten, F.  An efficient stochastic approximation em algorithm using conditional particle filters.  In _IEEE International Conference on Acoustics, Speech and Signal Processing_ , pp. 6274–6278, 2013.
  * Lindsten & Schön (2013) Lindsten, F. and Schön, T. B.  Backward simulation methods for Monte Carlo statistical inference.  _Foundations and Trends® in Machine Learning_ , 6(1):1–143, 2013.
  * Liu & Nocedal (1989) Liu, D. C. and Nocedal, J.  On the limited memory bfgs method for large scale optimization.  _Mathematical Programming_ , 45(1-3):503–528, 1989.
  * Ljung (1999) Ljung, L.  System identification.  _Wiley Encyclopedia of Electrical and Electronics Engineering_ , pp. 1–19, 1999.
  * Ljung (2010) Ljung, L.  Perspectives on system identification.  _Annual Reviews in Control_ , 34:1–12, 2010.
  * Lutter et al. (2019a) Lutter, M., Listmann, K., and Peters, J.  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.
  * Lutter et al. (2019b) Lutter, M., Ritter, C., and Peters, J.  Deep Lagrangian networks: Using physics as model prior for deep learning.  In _International Conference on Learning Representations (ICLR)_ , 2019b.
  * Marquardt (1963) Marquardt, D. W.  An algorithm for least-squares estimation of nonlinear parameters.  _Journal of the Society for Industrial and Applied Mathematics_ , 11(2):431–441, June 1963.  doi: 10.1137/0111030.
  * Menda et al. (2020) Menda, K., de Becdelièvre, J., Gupta, J. K., Kroo, I., Kochenderfer, M. J., and Manchester, Z.  Normalized Stanford helicopter dataset, February 2020.  URL https://doi.org/10.5281/zenodo.3662987.
  * Naesseth et al. (2015) Naesseth, C. A., Lindsten, F., and Schön, T. B.  Nested sequential Monte Carlo methods.  _International Conference on Machine Learning (ICML)_ , 2015.
  * Neal & Hinton (1998) Neal, R. M. and Hinton, G. E.  A view of the EM algorithm that justifies incremental, sparse, and other variants.  In _Learning in Graphical Models_ , pp. 355–368. Springer, 1998\.
  * Ordonez et al. (2017) Ordonez, C., Gupta, N., Reese, B., Seegmiller, N., Kelly, A., and Collins Jr, E. G.  Learning of skid-steered kinematic and dynamic models for motion planning.  _Robotics and Autonomous Systems_ , 95:207–221, 2017.
  * Punjani & Abbeel (2015) Punjani, A. and Abbeel, P.  Deep learning helicopter dynamics models.  In _IEEE International Conference on Robotics and Automation (ICRA)_ , pp. 3223–3230, 2015.
  * Rauch et al. (1965) Rauch, H. E., Tung, F., and Striebel, C. T.  Maximum likelihood estimates of linear dynamic systems.  _AIAA Journal_ , 3(8):1445–1450, 1965.
  * Schön et al. (2011) Schön, T. B., Wills, A., and Ninness, B.  System identification of nonlinear state-space models.  _Automatica_ , 47(1):39–49, 2011.
  * Shi et al. (2016) Shi, H.-J. M., Tu, S., Xu, Y., and Yin, W.  A primer on coordinate descent algorithms.  _arXiv preprint arXiv:1610.00040_ , 2016.
  * Snyder et al. (2008) Snyder, C., Bengtsson, T., Bickel, P., and Anderson, J.  Obstacles to high-dimensional particle filtering.  _Monthly Weather Review_ , 136(12):4629–4640, 2008.
  * Sun et al. (2008) Sun, X., Jin, L., and Xiong, M.  Extended kalman filter for estimation of parameters in nonlinear state-space models of biochemical networks.  _PLOS One_ , 3(11):e3758, 2008.
  * Svensson & Lindsten (2018) Svensson, A. and Lindsten, F.  Learning dynamical systems with particle stochastic approximation em.  _arXiv preprint arXiv:1806.09548_ , 2018.
  * Svensson et al. (2018) Svensson, A., Lindsten, F., and Schön, T. B.  Learning nonlinear state-space models using smooth particle-filter-based likelihood approximations.  _IFAC-PapersOnLine_ , 51(15):652–657, 2018.
  * Takens (1981) Takens, F.  Detecting strange attractors in turbulence.  In _Dynamical Systems and Turbulence_ , pp. 366–381. Springer, 1981\.
  * Van Dooren (1981) Van Dooren, P.  A generalized eigenvalue approach for solving riccati equations.  _SIAM Journal on Scientific and Statistical Computing_ , 2(2):121–135, 1981.
  * Van Overschee & De Moor (1994) Van Overschee, P. and De Moor, B.  N4SID: Subspace algorithms for the identification of combined deterministic-stochastic systems.  _Automatica_ , 30(1):75–93, 1994.
  * Zimmermann & Neuneier (2000) Zimmermann, H. and Neuneier, R.  Modeling dynamical systems by recurrent neural networks.  _WIT Transactions on Information and Communication Technologies_ , 25:557–566, 2000.

##  Appendix 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://heli.stanford.edu/. A gas-powered helicopter was flown by a
professional pilot to collect a large dataset of $6290$\mathrm{s}$$ 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 to an
accelerometer, a gyroscope, a magnetometer and vision cameras. Abbeel et al.
(2010) provide the raw data, as well as states estimates in the Earth
reference frame obtained with extended Kalman smoothing. Following Punjani &
Abbeel (2015), 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:

$\dot{s}=F(s,\delta)=\begin{bmatrix}C_{12}v\\\ \frac{1}{2}\hat{\omega}q\\\
C_{12}^{\top}g-\omega\times v+f_{v}(s,\delta)\\\
f_{\omega}(s,\delta)\end{bmatrix}$ (18)

where $s\in\mathbb{R}^{13}$ is the helicopter state consisting of its position
$r$, quaternion-attitude $q$, linear velocity $v$, angular velocity $\omega$,
and $\delta\in\mathbb{R}^{4}$ to be the control command. $C_{12}$ is the
rotation-matrix from the body to Earth reference frame, and $f_{v}$ and
$f_{\omega}$ are the linear and angular accelerations caused by aerodynamic
forces, and are what we aim to predict.

The above notation is related to that used in Section 4.2 as follows:

  * •

We define $u$ as the concatenation of all inputs to the model, including the
relevant state variables $v$ and $\omega$ and control commands $\delta$.

  * •

We define $y$ as the output predicted, which would correspond to a
concatenation of $f_{v}$ and $f_{\omega}$.

  * •

We define $x$ as the vector of unobserved flow states to be estimated and is
not present in their model.

The processed dataset used in our experiments can be found at
https://doi.org/10.5281/zenodo.3662987 (Menda et al., 2020).

###  A.2 Training Helicopter Models

Neural networks in the naive and H25 models have eight hidden layers of size
$32$ each, and tanh non-linearities. We optimize these models using an Adam
optimizer (Kingma & Ba, 2015) with a harmonic learning rate decay, and mini-
batch 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 CE-EM, using
$\rho_{x}=\rho_{\theta}=0.5$, $\sigma_{w}=\sigma_{v}=1.0$, and use an Adam
optimizer to optimize Equation 11 in the learning step. The learning rate for
dynamics parameters in $\theta_{\text{\tiny NL}}$ is $5.0\times 10^{-4}$ and
observation parameters in $\theta_{\text{\tiny NL}}$ is $1.0\times 10^{-3}$.
For its relative robustness, we optimize Equation 10 using a non-linear least
squares optimizer with a Trust-Region Reflective algorithm (Jones et 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 an EKF
depends on several user-provided parameters:

  * •

$x_{0}$: value of the initial state

  * •

$\Sigma_{0}$: 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 $\Sigma_{0}$ are all set to the
identity matrix. $x_{0}$ is assumed to be 0 on all dimensions.

A well-tuned EKF with an inaccurate initial state value converges to accurate
estimations in only a few time steps of transient behavior. Since the H25
model needs 25 past inputs to predict its first output prediction, we drop the
first 25 predictions from the EKF when computing RMSE, thereby omitting some
of the transient regime.

###  A.3 Figures

Figure 5: Test performance of optimized models on various trajectories

Figure 6: Predicted acceleration along axis $x$, $y$ and $z$ in the body
frame. For subspace identification and models trained by CE-EM, this plot
requires running an extended Kalman filter. These figures can be reproduced
for any other trajectory with the included code.

Figure 7: Predicted circular accelerations around axis $x$, $y$ and $z$ in the
body frame. For subspace identification and models trained by CE-EM, this plot
requires running an extended Kalman filter. These figures can be reproduced
for any other trajectory with the included code.