System identification is the process of building a mathematical model of an unknown system from measurements of its inputs and outputs. It is a key step for model-based control, estimator design, and output prediction. This work presents an algorithm for non-linear offline system identification from partial observations, i.e. situations in which the system's full-state is not directly observable. The algorithm presented, called SISL, iteratively infers the system's full state through non-linear optimization and then updates the model parameters. We test our algorithm on a simulated system of coupled Lorenz attractors, showing our algorithm's ability to identify high-dimensional systems that prove intractable for particle-based approaches. We also use SISL to identify the dynamics of an aerobatic helicopter. By augmenting the state with unobserved fluid states, we learn a model that predicts the acceleration of the helicopter better than state-of-the-art approaches. The performance of controllers and state-estimators for non-linear systems depends heavily on the quality of the model of system dynamics . System-identification tackles the problem of learning or calibrating dynamics models from data , which is often a timehistory of observations of the system and control inputs. In this work, we address the problem of learning dynamics models of partially observed, high-dimensional non-linear systems. That is, we consider situations in which the system's state cannot be inferred from a single observation, but instead requires inference over a time-series of observations. The problem of identifying systems from partial observations arises in many robotics domains (; ;). Though we often have direct measurements of a robot's pose and velocity, in many cases we cannot directly observe relevant quantities such as the temperature of actuators or the state the environment around the robot. Consider learning a dynamics model for an aerobatic helicopter. attempted to map only the helicopter's pose and velocity to its acceleration and they found their model to be inaccurate when predicting aggressive maneuvers. They posited that the substantial airflow generated by the helicopter affected the dynamics. Since it is often impossible to directly measure the state of the airflow around a vehicle, identification must be performed in a partially observed setting. System-identification is a mature field with a rich history (; 2010). Various techniques can be classified by whether they apply to linear or non-linear systems, with partially or fully observed states. Additionally, techniques are applied in an online or batch-offline setting. This work presents an approach to offline identification of non-linear and partially observed systems. When a system is fully observed, i.e. its full state is observed but corrupted by noise, a set of techniques called equation-error methods are typically employed (Åström &). In such cases, we can consider observations as independent, and minimize the error between the observed statederivatives and those predicted by the model given the control input and observed states. In partially observed settings, merely knowing the current input is insufficient to accurately predict the observation. Several black-box approaches exist to predict observations from time-series of inputs. Autoregressive approaches directly map a time-history of past inputs to observations . Recurrent neural networks and subspace-identification methods can also be used to learn blackbox dynamical systems from this data. However, in many cases prior knowledge can be used to specify structured, parameterized models of the system . Such models can be trained with less data and used with a wider array of control and state-estimation techniques than non-linear black-box models (; b; a). Techniques used to identify partially observed structured models are often based on Expectation-Maximization (EM) (; Schön et al., 2011; ;). An alternating procedure is performed in which a smoothing step uses the current system dynamics estimate to infer the distribution over state-trajectories, and is followed by a learning step that uses this distribution to update the system dynamics estimate. In the non-linear or non-Gaussian case, it is typically not possible to analytically characterize the distribution over trajectories, and thus methods based on Sequential Monte-Carlo such as Particle Smoothing (PS) (Schön et al., 2011;), or Extended Kalman Smoothing (EKS) are employed in the E-step. Though considered state-of-theart for this problem, both methods become intractable in high-dimensional state spaces. PS suffers from the curse of dimensionality, requiring an intractably large number of particles if the state space is high-dimensional , and an M-step that can be quadratic in complexity with respect to the number of particles (Schön et al., 2011). EKS-based methods are fast during the E-step, but the M-step requires approximations to integrate out state uncertainty, such as fitting non-linearities with Radial Basis Function approximators, and scales poorly with the dimension of the state-space . In this work, we present a system-identification algorithm that is suited for high-dimensional, nonlinear, and partially observed systems. By assuming that the systems are close to deterministic, as is often the case in robotics, we approximate the distribution over unobserved states using only their maximum-likelihood (ML) point-estimate. Our algorithm, called SISL (System-identification via Iterative Smoothing and Learning) performs the following two steps until convergence: • In the smoothing or E-step, we use non-linear programming to tractably find the ML pointestimate of the unobserved states. • In the learning or M-step, we use the estimate of unobserved states to improve the estimate of system parameters. The idea to use an ML point-estimate in lieu of the distribution over unobserved states in the EM procedure's E-step is not new, and, in general, does not guarantee monotonic convergence to a local optimum . However, such an approximation is equivalent to regular EM if the ML point-estimate is the only instance of unobserved variables with non-negligible probability . We apply this idea to the problem of system-identification for nearly deterministic systems, in which ML point-estimates can serve as surrogates for the true distribution over unobserved state-trajectories. The primary contribution of this work is an algorithm for identifying non-linear, partially observed systems that is able to scale to high-dimensional problems. In Section 2, we specify the assumptions underpinning our algorithm and discuss the computational methodology for using it. In Section 3, we empirically demonstrate that it is able to identify the parameters of a high-dimensional system of coupled Lorenz attractors, a problem that proves intractable for particle-based methods. We also demonstrate our algorithm on the problem of identifying the dynamics of an aerobatic helicopter, and compare against various approaches including the state-of-the-art approach . In this work, we assume that we are given a batch of trajectories containing observations y 1:T ∈ R m×T of a dynamical system as it evolves over a time horizon T, possibly forced by some known input sequence u 1:T −1. We assume that this dynamical system has a state x ∈ R n that evolves and generates observations according to the following equations, where w t is referred to as the process noise and v t as the observation noise. w t and v t are both assumed to be additive for notational simplicity, but this is not a required assumption. 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 θ (x, t) and g θ (x, t) for θ ∈ Θ that approximate the dynamical system's evolution and observation processes. The goal of our algorithm is to find the parameters θ that maximize the likelihood of the observations. That is, we seek to find: Assuming, the system is Markovian in x t, we can factorize the distributions as: In order to tractably solve the maximization problem in Equation, particle-based EM techniques typically approximate the integral as an expectation over a particle set (Schön et al., 2011;). However, since particle-based methods can struggle with high-dimensional spaces , this work seeks point estimates for the likelihoodmaximizing state sequence, which can be found using non-linear programming. Instead of solving the full maximum likelihood problem shown in Equation, we solve a surrogate problem: The maximized surrogate objective is equivalent to the original objective if the distribution over x 1:T given y 1:T and θ is well-approximated by its ML point-estimate, i.e. the distribution is a Dirac delta function. We expect this assumption to hold if the trajectories are sufficiently long, the dynamics are close to deterministic, and all relevant dynamic modes are persistently excited. Making these assumptions, we now present an algorithm for solving Equation. By taking the logarithm of the likelihood-objective in Equation, and using the factorization from Equation, we get the SISL objective J(x 1:T, θ) as follows: Jointly maximizing this objective over x 1:T and θ yields θ SISL. Though this objective can be optimized as is by a non-linear optimizer, it is not necessarily efficient to do so since θ and x 1:T are highly coupled, leading to inefficient and potentially unstable updates. For this reason, we take an iterative approach akin to EM, which performs block-coordinate ascent on the J(x 1:T, θ). At iteration k, we first perform smoothing by holding θ constant, and finding the ML point-estimate for x 1:T as follows: Here ρ x scales a soft trust-region regularizer similar to damping terms found in LevenbergMarquardt methods . The smoothing problem can be solved in O(n 2 T) by taking advantage of sparsity in the smoothing objective's Hessian matrix. In the learning step, we hold x 1:T constant and find: Again, ρ θ scales a soft trust-region regularizer, and specifying log p(θ) allows us to regularize θ toward a prior. The above optimization problem can be solved using any non-linear optimizer, such as a first or second-order gradient descent scheme. The SISL algorithm iterates between the smoothing and learning steps until convergence. The objective of our experiments is to demonstrate that the SISL 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 by attempting to characterize the interaction of the helicopter with the fluid around it, without having any direct observation of the fluid state. In this didactic experiment, we will show that: 1. SISL learns unbiased parameter estimates of systems that are close to deterministic, and, 2. SISL scales to high-dimensional problems in which particle-smoothing is 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 increase in state dimensionality can be achieved by coupling multiple individual attractors. The state of a system with K coupled Lorenz attractors is x ∈ R 3K = {. . ., x 1,k, x 2,k, x 3,k, . . .}. The dynamics of the system are as follows: where H is an R 3K×3K matrix. We nominally set the parameters (σ k, ρ k, β 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 ∈ R (3K−2) is found from x as follows: where C ∈ R (3K−2)×3K is a known matrix with full row-rank, and v is the observation noise sampled from a Gaussian with diagonal covariance σ 2 v 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 ∆t = 0.04s, and integrate the system using a 4 th -order Runge-Kutta method. Initial conditions for each trajectory are sampled such that x 1,k ∼ N (−6, 2.5 2), x 2,k ∼ N (−6, 2.5 2), x 3,k ∼ N (24, 2.5 2). To test the conditions under which SISL learns unbiased parameter estimates, we simulate a single Lorenz system with H = 0, and known C ∈ R 2×3. We introduce and vary the process noise w ∼ N (0, σ 2 w I), and vary the observation noise coefficient σ v, and then attempt to estimate the parameters (σ, ρ, β). Using initial guesses within 10% of the system's true parameter values, we run SISL on a single sampled trajectory. For each choice of σ w and σ v, we repeat this process for 10 random seeds. Table 1 shows the mean and standard errors of parameter estimates for various σ w and σ v. We highlight in red the mean estimates that are not within two standard errors of their true value. We see that σ and ρ are estimated without bias for all scenarios. However, the estimate of β 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 SISL is sound when systems evolve close to deterministically, but can be biased if it is not. State-of-the-art methods for parameter identification of partially observed systems rely on particlesmoothing (PS) to estimate a distribution over states x 1:T given a trajectory of observations y 1:T. We will experimentally demonstrate that the performance of PS does not scale to high-dimensional systems. To do so, we compare the number of particles required for PS to reliably characterize the distribution over x 1:T for a system of two and six coupled Lorenz systems. The systems are simulated with observation noise σ v = 0.01 but without process noise. We implement a PS as specified by Schön et al.. Letx t,n, w t,n be the estimated system-state and particle weight corresponding to the nth particle at the tth timestep. In order to test whether PS can reliably characterize the posterior distribution over x 1:T, we measure the weighted root mean square error (RMSE), ξ as follows: If PS is reliable using N particles, we should see that ξ N is tightly distributed across random seeds. These experimentally demonstrate that, for a sufficiently non-linear and partially observable system, the number of particles required to reliably characterize the posterior distribution over hidden states grows intractably with the dimension of the system. Since particle-based EM methods for system-identification are typically super-linear in complexity with respect to the number of particles (Schön et al., 2011), these methods are ill-suited to high-dimensional problems. To demonstrate that SISL is capable of identifying high-dimensional systems, we show that we can estimate the dynamics of an 18 dimensional system of six coupled Lorenz attractors. Moreover, as the number of trajectories provided to SISL increases, it converges to more accurate estimates. To test this claim, we sample 2, 4, and 8 trajectories from a system with parameters θ true, and σ v = 0.01. We randomly initialize each element of the parameters being optimized (θ = [σ 1:K, ρ 1:K, β 1:K, H]) to within 10% of the their value in θ true. We then run SISL on each batch, tracking the error in the estimated dynamics as training proceeds. We measure this error, which we call (θ), as follows: In the learning step, we do not regularize θ to a prior and set ρ θ = 0. In Figure 2, we see the 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 SISL convergences monotonically to a local optimum in all cases. This experiment supports our claim that SISL is able to identify the parameters of a high-dimensional, non-linear, and partially observable system that is intractable for particle-based methods. The experiments conducted thus far have demonstrated that SISL 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 SISL to characterize the effect of unobserved states on the dynamics of an aerobatic helicopter. Characterizing the dynamics of a helicopter undergoing aggressive aerobatic maneuvers is widely considered to be a challenging system-identification problem . 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, 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 . Instead of attempting to estimate the unobserved fluid state, they directly learn a mapping from a 0.5 s long history of observed state measurements and control commands to the forces acting on the helicopter. This approach is sensible, and is equivalent to considering the past 0.5 s of observations as the system's state. However, it can require a very large number of lagged observations to represent complex phenomena. Having such a high dimensional state can make the control design and stateestimation 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. 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 . Trajectories are split into 10 s long chunks and then randomly distributed into train, test, and validation sets according to the protocol established by; and summarized in Appendix 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, RMSE =, where y (measured) t is the measured force from the dataset, y (pred) t is the force predicted by the model, and T is the number of time-steps in each trajectory. 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 = NN θn (u t), where NN θn is a neural-network with parameters θ n. We refer to this model as the naive model. We also compare to the work of. They predict y t using a time-history u t−H:t of H 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 = NN θ h (u t−H:t) with parameters θ 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 Appendix A.2. As a third baseline, we compare with subspace-identification methods . We letỹ t = y t − NN θ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: Here, x ∈ R d is the unobserved state with arbitrary dimension d. The learned parameters are 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. Particle-based EM methods are not presented as baselines because they are intractable on problems with large state-spaces, as shown in Section 3.1.2. Similar to the parameterization used for subpace-identification, we fit the prediction errors of the naive model using the following dynamical system: where NNθ NL is a neural network, and θ NL = [A NL, B NL, C NL, D NL,θ NL] are the learned parameters. While learning, we assume that both process and observation noise are distributed with diagonal Gaussian covariance matrices σ w I and σ v I respectively. The values of σ w and σ v are treated as hyperparmeters of SISL. 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 Appendix 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 SISL. The test RMSE of the naive and H25 models can be evaluated directly on the test trajectories using next-step prediction. However, the SID and 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ỹ 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ỹ t. Figure 3a 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 can be found in Appendix A.3. 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 3b 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 SISL to identify a non-linear dynamical model of unobserved states that affect the forces acting an aerobatic helicopter. 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 our algorithm can perform identification on a problem that particle-based EM methods are fundamentally ill-suited for. We also validated that our algorithm is an effective replacement for identification methods based on EM if the system is close to deterministic, but can yield biased parameter estimates if it is not. We then used our algorithm to model the time-varying hiddenstates that affect the dynamics of an aerobatic helicopter. Our approach outperforms state-of-the-art methods because it is able to fit large non-linear models to unobserved states. We aim to apply our algorithm to system identification problems in a number of domains. There has recently been interest in characterizing the dynamics of aircraft with high aspect ratios, for which the difficult-to-observe bending modes substantially impact dynamics. Additionally, the inability to measure friction forces in dynamic interactions involving contact typically stands in the way of system identification, and thus requires algorithms that are capable of identification under partial observation. A APPENDIX In this work we use the dataset gathered by and available at http://heli. stanford.edu/. A gas-powered helicopter was flown by a professional pilot to collect a large dataset of 6290s of flight. There are four controls: the longitudinal and lateral cyclic pitch, the tail rotor pitch and the collective pitch. The state is measured thanks to an accelerometer, a gyroscope, a magnetometer and vision cameras. provide the raw data, as well as states estimates in the Earth reference frame obtained with extended Kalman smoothing.'s protocol, we use the fused sensor data and downsample it from 100Hz to 50Hz. From the Earth frame accelerations provided in the dataset, we compute body frame accelerations (minus gyroscopic terms) which are the prediction targets for our training. Using the notations from , we can write the helicopter dynamics in the following form: where s ∈ R 13 is the helicopter state consisting of its position r, quaternion-attitude q, linear velocity v, angular velocity ω, and δ ∈ 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 ω are the linear and angular accelerations caused by aerodynamic forces, and are what we aim to predict. This notation connects with the one used in our paper in the following way: • We define u as the concatenation of all inputs to the model, including the relevant state variables v and ω and control commands δ. • We define y as the output predicted, which would correspond to a concatenation of f v and f ω. • We define x as the vector of unobserved flow states to be estimated and is not present in their model. 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 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 SISL, using ρ x = ρ θ = 0.5, σ w = σ v = 1.0, and use an Adam optimizer to optimize Equation in the learning step. The learning rate for dynamics parameters in θ NL is 5.0 × 10 −4 and observation parameters in θ NL is 1.0 × 10 −3. For its relative robustness, we optimize Equation using a non-linear least squares optimizer with a Trust-Region Reflective algorithm (-) 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 • Σ 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 Σ 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.
This work presents a scalable algorithm for non-linear offline system identification from partial observations.