Identification of LTV Dynamical Models with Smooth or Discontinuous Time Evolution by means of Convex Optimization
Fredrik Bagge Carlson Anders Robertsson Rolf Johansson

propositionProposition
corollaryCorollary
theoremTheorem




Slowdown of the surface diffusion during early stages of bacterial colonization
    C. Cassier-Chauvat
    December 30, 2023
===============================================================================

empty
empty


    We establish a connection between trend filtering and system identification which results in a family of new identification methods for linear, time-varying (LTV) dynamical models based on convex optimization. We demonstrate how the design of the cost function promotes a model with either a continuous change in dynamics over time, or causes discontinuous changes in model coefficients occurring at a finite (sparse) set of time instances. We further discuss the introduction of priors on the model parameters for situations where excitation is insufficient for identification. The identification problems are cast as convex optimization problems and are applicable to, e.g., ARX models and state-space models with time-varying parameters. We illustrate usage of the methods in simulations of jump-linear systems, a nonlinear robot arm with non-smooth friction and stiff contacts as well as in model-based, trajectory centric reinforcement learning on a smooth nonlinear system.
    











§ INTRODUCTION


The difficulty of the task of identifying time-varying dynamical models of systems varies greatly with the model considered and the availability of measurements of the state sequence. For smoothly changing dynamics, linear in the parameters, the recursive least-squares algorithm with exponential forgetting (RLSλ) is a common option. If a Gaussian random-walk model for the parameters is assumed, a Kalman filtering/smoothing algorithm <cit.> gives the filtering/smoothing densities of the parameters in closed form. The assumption of smoothly (Gaussian) varying dynamics is often restrictive. Discontinuous dynamics changes occur, for instance, when an external controller changes operation mode, when a sudden contact between a robot and its environment is established, an unmodeled disturbance enters the system or when a system is suddenly damaged.

Identification of systems with non-smooth dynamics evolution has been studied extensively. The book <cit.> treats the case where the dynamics are known, but the state sequence unknown, i.e., state estimation. In <cit.>, the authors examine the residuals from an initial constant dynamics fit to determine regions in time where improved fit is needed by the introduction of additional constant dynamics models. Results on identifiability and observability in jump-linear systems in the non-controlled (autonomous) setting are available in <cit.>. The main result on identifiability in <cit.> was a rank condition on a Hankel matrix constructed from the collected output data, similar to classical results on the least-squares identification of ARX models which appears as rank constraints on the, typically Toeplitz or block-Toeplitz, regressor matrix. Identifiability of the methods proposed in this article are discussed in <ref>.

An LTV model can be seen as a first-order approximation of the dynamics of a nonlinear system around a trajectory. We emphasize that such an approximation will in general fail to generalize far from the this trajectory, but many methods in reinforcement learning and control make efficient use of the linearized dynamics for optimization, while ensuring validity of the approximation by constraints or penalty terms. An example provided in <ref> highlights such a method.

An important class of identification methods that has been popularized lately is trend filtering methods <cit.>. Trend filtering methods work by specifying a fitness criterion that determines the goodness of fit, as well as a regularization term, often chosen with sparsity promoting qualities. As a simple example, consider the reconstruction ŷ of a noisy signal y = {y_t}_t=1^T with piecewise constant segments. To this end, we may formulate and solve the convex optimization problem

    ŷy-ŷ^2 + λ∑_t |ŷ_t+1 - ŷ_t|

The first term is the fitness criterion or loss function, whereas the second term is a sparsity-promoting regularizer which promotes small changes between consecutive samples in the reconstructed signal. The sparsity promoting effect of the 1-norm regularizer is well known, and stems from the constant length of the gradient whenever the argument is non-zero <cit.>. Compare this to the squared difference, for which the gradient rapidly vanishes as the argument approaches zero. The squared difference will thus promote small arguments, whereas the 1-norm promotes sparse arguments.

In this work, we will draw inspiration from the trend-filtering literature to develop new system identification methods for LTV models, with interesting properties. In trend filtering, we decompose a curve as a set of polynomial segments. In the identification methods proposed in this work, we instead decompose a multivariable state sequence as the output of a set of LTV models, where the model coefficients evolve as polynomial functions of time. We start by defining a set of optimization problems with a least-squares loss function and carefully chosen regularization terms. We further discuss how prior information can be utilized to increase the accuracy of the identification and end the article with identification of a nonlinear system with non-smooth friction and an example of model-based reinforcement learning followed by a discussion.



§ LTI IDENTIFICATION


We start by considering the case of identification of the parameters in an LTI model on the form

    x_t+1 = A x_t + B u_t + v_t,    t ∈ [1,T]

where xn, um are the state and input respectively. A discussion around the noise term v_t is deferred until <ref>, where we indicate how statistical assumptions on v_t influence the cost function and the properties of the estimate. If the state and input sequences are known, a plethora of methods for estimating the parameters exists. A common method for systems that are linear in the parameters is the least-squares (LS) method, which in case of Gaussian noise, v, coincides with the maximum likelihood (ML) estimate. To facilitate estimation using the LS method, we write the model on the form =, and arrange the data according to

    =
       [ x_1;   ⋮; x_T ]      Tn
       = A⃗ ⃗ ⃗ ⃗ ⃗B⃗      K
       =
       [   I_n ⊗ x_0   I_n ⊗ u_0;           ⋮           ⋮; I_n ⊗ x_T-1 I_n ⊗ u_T-1 ]      ∈ℝ^Tn× K

where ⊗ denotes the Kronecker product and K=n^2+nm is the number of model parameters, and solve the optimization problem eq:lscost with closed-form solution eq:ls.

    ^*    = _ - ^2 
        =






§ TIME-VARYING DYNAMICS

We now move on to the contribution of this work, and extend our view to systems where the dynamics change with time. We limit the scope of this article to models on the form

    x_t+1   = A_t x_t + B_t u_t + v_t
    _t    = A⃗_⃗t⃗ ⃗ ⃗ ⃗ ⃗B⃗_⃗t⃗

where the parameters  are assumed to evolve according to the dynamical system

    k_t+1   = H_t k_t + w_t
    
       y_t    = (I_n ⊗x_t    u_t) _t

where, if no prior knowledge is available, the dynamics matrix H_t can be taken as the identity matrix; H = I implies that the model coefficients follow a random walk dictated by the properties of w_t, i.e., the state transition density function p_w(k_t+1|k_t). The emission density function p_v(x_t+1 | x_t, u_t, k_t) is determining the drift of the state, which for the parameter estimation problem can be seen as the distribution of measurements, given the current state of the system.
We emphasize here that the state in the parameter evolution model refers to the current parameters k_t and not the system state x_t, hence, p_v is called the emission density and not the transition density. Particular choices of p_v and p_w emit data likelihoods concave in the parameters and hence amenable to convex optimization.


The following sections will introduce a number of optimization problems with different regularization functions, corresponding to different choices of p_w, and different regularization arguments, corresponding to different choices of H. We also discuss the quality of the identification resulting from the different modeling choices.



 §.§ Low frequency time evolution

A slowly varying signal is characterized by small first-order time differences.
To identify slowly varying dynamics parameters, we thus penalize the squared 2-norm of the first-order time difference of the model parameters, and solve the optimization problem

    -^2 + λ^2∑_t _t+1 - _t^2

where ∑_t denotes the sum over relevant indices t, in this case t∈ [1,T-1].
This optimization problem has a closed form solution given by

    ^*    = ( + λ^2 D_1 D_1)^-1Ỹ
       = vec(_1, ... , _T)

where  and Ỹ are appropriately constructed matrices and the first-order differentiation operator matrix D_1 is constructed such that λ^2D_1 ^2 equals the second term in eq:slow.
The computational complexity 𝒪((TK)^3) of computing k^* using the closed-form solution eq:closedform becomes prohibitive for all but toy problems. We note that the cost function in eq:slow is the negative data log-likelihood of a Brownian random-walk parameter model with H=I, which motivates us to develop a dynamic programming algorithm based on a Kalman smoother, detailed in <ref>.



 §.§ Smooth time evolution

A smoothly varying signal is characterized by small second-order time differences.
To identify smoothly time-varying dynamics parameters, we thus penalize the squared 2-norm of the second-order time difference of the model parameters, and solve the optimization problem

    -^2 + λ^2∑_t _t+2 -2 _t+1 + _t^2

Also this optimization problem has a closed form solution on the form eq:closedform with the corresponding second-order differentiation operator D_2. <Ref> is the negative data log-likelihood of a Brownian random-walk parameter model with added momentum and H derived in <ref>, where a Kalman smoother with augmented state is developed to find the optimal solution. We also extend problem eq:smooth to more general regularization terms in <ref>.





 §.§ Piecewise constant time evolution

In the presence of discontinuous or abrupt changes in the dynamics, estimation method eq:smooth might perform poorly. A signal which is mostly flat, with a small number of distinct level changes, is characterized by a sparse first-order time difference. To detect sudden changes in dynamics, we thus formulate and solve the problem

    -^2 + λ∑_t _t+1 - _t

We can give eq:pwconstant an interpretation as a grouped-lasso cost function, where instead of groups being formed out of variables, our groups are defined by differences between variables. We thus have a penalty on the 1-norm on the length of the difference vectors _t+1 - _t since ·_1 = ·.
The 1-norm is a sparsity-promoting penalty, hence a solution in which only a small number of non-zero first-order time differences in the model parameters is favored, i.e., a piecewise constant dynamics evolution.
At a first glance, one might consider the formulation

    -^2 + λ∑_t _t+1 - _t_1

which results in a dynamics evolution with sparse changes in the coefficients, but changes to different entries of _t are not necessarily occurring at the same time instants. The formulation eq:pwconstant, however, promotes a solution in which the change occurs at the same time instants for all coefficients in A and B, i.e., _t+1 = _t for most t.




  §.§.§ Implementation

Due to the non-squared norm penalty ∑_t _t+1 - _t, problem eq:pwconstant is significantly harder to solve than eq:smooth. An efficient implementation using the linearized ADMM algorithm <cit.> is made available in the accompanying repository.




 §.§ Piecewise constant time evolution with known number of steps


If the number of switches in dynamics parameters, M, is known in advance, the optimal problem to solve is

    -^2
             ∑_t 1{_t+1≠_t}≤ M

where 1{·} is the indicator function.
This problem is non-convex and we propose solving it using dynamic programming (DP). For this purpose we modify the algorithm developed in <cit.>, an algorithm frequently referred to as segmented least-squares <cit.>, where a curve is approximated by piecewise linear segments. The modification lies in the association of each segment (set of consecutive time indices during which the parameters are constant) with a dynamics model, as opposed to a simple straight line.[Indeed, if a simple integrator is chosen as dynamics model and a constant input is assumed, the result of our extended algorithm reduces to the segmented least-squares solution.] Unfortunately, the computational complexity of the dynamic programming solution, 𝒪(T^2K^3), becomes prohibitive for large T.[For details regarding the DP algorithm and implementation, the reader is referred to the source-code repository accompanying this article.]



 §.§ Piecewise linear time evolution

A piecewise linear signal is characterized by a sparse second-order time difference, i.e., it has a small number of changes in the slope. A piecewise linear time-evolution of the dynamics parameters is hence obtained if we solve the optimization problem.

    -^2 + λ∑_t _t+2 -2 _t+1 + _t




 §.§ Summary

The proposed optimization problems are summarized in <ref>. The table illustrates how the choice of regularizer and order of time-differentiation of the parameter vector affects the quality of the resulting solution.






 §.§ Two-step refinement


Since many of the proposed formulations of the optimization problem penalize the size of the changes to the parameters, solutions in which the changes are slightly underestimated are favored. To mitigate this issue, a two-step procedure can be implemented where in the first step, change points (knots) are identified. In the second step, the penalty on the one-norm is removed and equality constraints are introduced between consecutive time-indices for which no change in dynamics was indicated by the first step.

The second step can be computed very efficiently by noticing that the problem can be split into several identical sub-problems at the knots identified in the first step. The sub-problems have closed-form solutions if the problem in <ref> is considered.

To identify the points at which the dynamics change, we observe the argument inside the sum of the regularization term, i.e., a_t1 = _t+1 - _t or a_t2 = _t+2 -2 _t+1 + _t. Time instances where a_t is taking non-zero values indicate change points.




§ DYNAMICS PRIOR AND KALMAN FILTERING

Reviewer found this section rushed
The identifiability of the parameters in a dynamical model hinges on the observability of the dynamics system eq:dynsys, or more explicitly, only modes excited by the input u will be satisfactorily identified. If the identification is part of an iterative learning and control scheme, e.g., ILC or reinforcement learning, it might be undesirable to introduce additional noise in the input to improve excitation for identification. This section will introduce prior information about the dynamics which mitigates the issue of poor excitation of the system modes. The prior information might come from, e.g., a nominal model known to be inaccurate, or an estimated global model such as a Gaussian mixture model (GMM). A statistical model of the joint density p(x_t+1,x_t,u_t) constructed from previously collected tuples (x_t+1,x_t,u_t) provides a dynamical model of the system through the conditional pdf p(x_t+1|x_t,u_t).

We will see that for priors from certain families, the resulting optimization problem remains convex. For the special case of a Gaussian prior over the dynamics parameters or the output, the posterior mean of the parameter vector is conveniently obtained from a Kalman-smoothing algorithm, modified to include the prior.



 §.§ General case


If we introduce a parameter state  (c.f., eq:tvk) and a prior over all parameter-state variables p(_t|z_t), where the variable z_t might be, for instance, the time index t or state x_t, we have the data log-likelihood



 

    log p(,y|x,z)_1:T   = ∑_t=1^T log p(y_t|_t,x_t) 
    
        + ∑_t=1^T-1log p(_t+1|_t)    + ∑_t=1^Tlog p(_t|z_t)

which factors conveniently due to the Markov property of a state-space model.
For particular choices of density functions in eq:lcdll, notably Gaussian and Laplacian, the negative likelihood function becomes convex. The next section will elaborate on the Gaussian case and introduce a recursive algorithm that solves for the full posterior efficiently. The Laplacian case, while convex, does not admit an equally efficient algorithm, but is more robust to outliers in the data.



 §.§ Gaussian case
 
If all densities in eq:lcdll are Gaussian and  is modeled with the Brownian random walk model eq:dynsys (Gaussian v_t), eq:lcdll can be written on the form (scaling constants omitted)

    -log p(,y|x,z)_1:T   =
        ∑_t=1^T y_t - ŷ(_t,x_t)^2_Σ^-1_y
       + ∑_t=1^T-1_t+1 - _t^2_Σ^-1_
       + ∑_t=1^Tμ_0(z_t) - _t^2_Σ^-1_0(z_t)

for some function μ_0(z_t) which produces the prior mean of  given z_t. Σ_y, Σ_, Σ_0(z_t) are the covariance matrices of the state-drift, parameter drift and prior respectively and x_Σ^-1^2 = xΣ^-1x.

In this special case, we introduce a recursive solution given by a modified Kalman smoothing algorithm, where the conditional mean of the state is updated with the prior. Consider the standard Kalman filtering equations, reproduced here to establish the notation
    
    x̂_t|t-1   = A x̂_t-1|t-1 + B u_t-1
    
            P_t|t-1   = A P_t-1|t-1A + R_1
    
            K_t    = P_t|t-1C(CP_t|t-1C+R_2)^-1
    x̂_t|t   = x̂_t|t-1 + K_t(y_t-Cx̂_t|t-1)
    
            P_t|t   = P_t|t-1 - K_t C P_t|t-1

    where x is the state vector, with state-drift covariance R_1 and C is a matrix that relates x to a measurement y=Cx with covariance R_2. The first two equations constitute the prediction step, and the last two equations incorporate the measurement y_t in the correction step.
    The modification required to incorporate a Gaussian prior on the state variable p(x_t|v_t) = (μ_0(v_t), Σ_0(v_t)) involves a repeated correction step and takes the form
    
    K̅_t    = P_t|t(P_t|t+Σ_0(v_t))^-1
    x̅_t|t   = x̂_t|t + K̅_t(μ_0(v_t)-x̂_t|t)
    P̅_t|t   = P_t|t - K̅_t P_t|t

    where ·̅ denotes the posterior value. This additional correction can be interpreted as receiving a second measurement μ_0(v_t) with covariance Σ_0(v_t).
    For the Kalman-smoothing algorithm, x̂_t|t and P_t|t in eq:postx,eq:postcov are replaced with x̂_t|T and P_t|T.

A prior over the output of the system, or a subset thereof, is straight forward to include in the estimation by means of an extra update step, with C,R_2 and y being replaced with their corresponding values according to the prior.



 §.§ Kalman filter for identification

We can employ the Kalman-based algorithm to solve two of the proposed optimization problems:


  §.§.§ Low frequency

The Kalman smoother can be used for solving identification problems like eq:slow by noting that eq:slow is the negative log-likelihood of the dynamics model eq:dynsys. The identification problem is thus reduced to a standard state-estimation problem.


  §.§.§ Smooth


To develop a Kalman-filter based algorithm for solving eq:smooth, we augment the model eq:dynsys with the state variable k^'_t = k_t - k_t-1 and note that k^'_t+1 - k^'_t = k_t+1 - 2k_t + k_t-1. We thus introduce the augmented-state model

    [   k_t+1; k^'_t+1 ]   = [ I_K I_K; 0_K I_K ][   k_t; k^'_t ] + [   0; w_t ]
    
       y_t    =  [ (I_n ⊗x_t    u_t)                 0 ][   k_t; k^'_t ]

which is on a form suitable for filtering/smoothing with the machinery developed above.



  §.§.§ General case

The Kalman-filter based identification method can be generalized to solving optimization problems where the argument in the regularizer appearing in eq:smooth is replaced by a general linear operation on the parameter vector, P(z)k, and we have the following proposition

    Any optimization problem on the form
    
    -^2 + λ^2∑_t P(z)_t^2

    where P(z) is a polynomial of degree n>0 in the time difference operator z with z^-n P(1) = 0, can be solved with a Kalman smoother employed to an autonomous state-space system.


    Let P^*(z^-1) = z^-n P(z). We assume without loss of generality that P^*(0) = 1 since any constant P^*(0) can be factored out of the polynomial. Q(z^-1) = P^*(z^-1)^-1 is a strictly proper transfer function and has a realization as a linear, Gaussian state-space system of degree n. Since Q(z^-1) is strictly proper, the realization has no direct term.
    The negative data log-likelihood of Q(z^-1) is equal, up to constants idenpendent of , to the cost function in <ref>, hence the Kalman smoother applied to Q optimizes <ref>.

For eq:smooth P(z) equals z^2 - 2z + 1 and Q(z^-1) has a realization on the form eq:dynsys2.





§ WELL-POSEDNESS AND IDENTIFIABILITY

To assess the well-posedness of the proposed identification methods, we start
by noting that the problem of finding A in x_t+1 = Ax_t given a pair
(x_t+1,x_t) is an ill-posed problem in the sense that the solution is non
unique. If we are given several pairs (x_t+1,x_t), for different t, while
A remains constant, the problem becomes over-determined and well-posed in the
least-squares sense, provided that the vectors {x_t^(i)}_t=1^T span ℝ^n. The LTI-case in <ref> is well posed according to
classical results, when Φ has full column rank.

When we extend our view to LTV models, the number of free
parameters is increased significantly, and the corresponding
regressor matrix  will never have full column
rank and the introduction of a regularization term is necessary.
Informally, for every n measurements, we have K=n^2+nm
free parameters. If we consider the identification problem
of <ref> and let
λ→∞, the regularizer terms essentially becomes equality constraints.
This will enforce a solution in which all parameters
in k are constant over time, and the problem reduces
to the LTI-problem. As λ decreases,
the effective number of free parameters increases until the
problem gets ill-posed for λ = 0.
We formalize the above arguments as

    Optimization problems eq:slow,eq:pwconstant have unique global minima for λ > 0 if and only if the corresponding LTI optimization problem has a unique solution.


    The cost function is a sum of two convex terms. For a global minimum to be non-unique, the Hessians of the two terms must have intersecting nullspaces.
    In the limit λ→∞ the problem reduces to the LTI problem. The nullspace of the regularization Hessian, which is invariant to λ, does thus not share any directions with the nullspace of  which establishes the equivalence of identifiability between the LTI problem and the LTV problems.


    Optimization problems eq:smooth,eq:pwlinear with higher order differentiation in the regularization term have unique global minima for λ > 0 if and only if there exists no vector v n+m such that
    
    C^xu_t v = x_tx_t    x_tu_t
     u_tx_t    u_tu_tv = 0  ∀ t



    Again, the cost function is a sum of two convex terms and for a global minimum to be non-unique, the Hessians of the two terms must have intersecting nullspaces.
    In the limit λ→∞ the regularization term reduces to a linear constraint set, allowing only parameter vectors that lie along a line through time. Let ṽ≠ 0 be such a vector, parametrized by t as ṽ = v̅    2v̅   ⋯    Tv̅TK where v̅ = vec({v}_1^N) K and v is an arbitrary vector n+m. ṽ∈null() implies that the loss is invariant to the pertubation αṽ to k̃ for an arbitrary α.  () is given by blkdiag({I_n ⊗ C^xu_t }_1^T) which means that ṽ∈null()⟺α t (I_n ⊗ C^xu_t) v̅ = 0  ∀ (α,t) ⟺v̅∈null(I_n ⊗ C^xu_t) ∀ t, which implies v ∈nullC^xu_t due to the block-diagonal nature of I_n ⊗ C^xu_t





For the LTI problem to be well-posed, the system must be identifiable and the input u must be persistently exciting of sufficient order <cit.>.




















§ EXAMPLE – JUMP-LINEAR SYSTEM

Add two link example where a stiff contact is established and there is a sign change in control signal and hence friction. Apply pwconstant method and maybe scrap linear system example. Fig 3 can be reproduced with analytical jacobian as comparison. Mention in abstract that we then evaluate on one non-linear but smooth system (pendcart) and one non-smooth system (twolink).
Only linear systems considered (pend-cart actually not linear. Highligt that pendcart is nonlinear or introduce simple robotic example. Mention early that an LTV system along a trajectory does not generalize far for a nonlinear system, hence the use of KL constraint.

We now consider a simulated example. We generate a state sequence from the following LTV system, where the change in dynamics, from

    A_t = [
    [ 0.95  0.1;  0.0 0.95;      ]],    B_t = [
    [ 0.2; 1.0;     ]]

to 
    A_t = [
    [  0.5 0.05;  0.0  0.5;      ]],    B_t = [
    [ 0.2; 1.0;     ]]

occurred at t=200.
The input was Gaussian noise of zero mean and unit variance, state transition noise and measuremet noise of zero mean and σ = 0.2 were added.
<Ref> depicts the estimated coefficients in the dynamics matrices for a value of λ chosen using the L-curve method <cit.>.





§ EXAMPLE – NON-SMOOTH ROBOT ARM WITH STIFF CONTACT

To illustrate the ability of the proposed models to represent the non-smooth dynamics along a trajectory of a robot arm, we simulate a two-link robot with discontinuous Coulomb friction. We also let the robot establish a stiff contact with the environment to illustrate both strengths and weaknesses of the modeling approach.

The state of the robot arm consists of two joint coordinates, q, and their time derivatives, q̇. <Ref> illustrates the state trajectories, control torques and simulations of a model estimated by solving eq:pwconstant. The figure clearly illustrates that the model is able to capture the dynamics both during the non-smooth sign change of the velocity, but also during establishment of the stiff contact. The learned dynamics of the contact is however time-dependent, which is illustrated in <Ref>, where the model is used on a validation trajectory where a different noise sequence was added to the control torque. Due to the novel input signal, the contact is established at a different time-instant and as a consequence, there is an error transient in the simulated data.





§ EXAMPLE – REINFORCEMENT LEARNING
 
More details needed
In this example, we use the proposed methods to identify LTV dynamics models for
reinforcement learning. The goal of the task is to dampen oscillations of a pendulum attached to a moving
cart by means of moving the cart, with bounds on the control signal and a quadratic
cost on states and control. Due to the nonlinear nature of the pendulum dynamics,
linear expansions of the dynamics in the upward (initial) position and downward (final) position have poles
on opposite sides of the imaginary axis. To this end, we employ a reinforcement-learning framework inspired by <cit.>, where we perform a series of rollouts whereafter each we
1) fit a dynamics model along the last obtained trajectory,
2) optimize the cost function under
the model using iterative LQG (differential dynamic programming),[Implementation made available at
github.com/baggepinnen/DifferentialDynamicProgramming.jlgithub.com/baggepinnen/DifferentialDynamicProgramming.jl] an algorithm that calculates the value function exactly under the LTV dynamics and a quadratic expansion of the cost function.
In order to stay close to the validity region of the linear model, we put bounds on the deviation between each new trajectory and the last trajectory.

We compare three different models;
the ground truth system model, an LTV model (obtained by solving eq:smooth) and an LTI model.
The total cost over T=500 time steps is shown as a function of learning iteration
in <ref>. The figure illustrates how the learning procedure reaches the
optimal cost of the ground truth model when an LTV model is used, whereas when
using an LTI model, the learning diverges. The figure further illustrates that if
the LTV model is fit using a prior (<ref>), the learning speed is increased. The prior in this case was constructed from the true system model, linearized around the last trajectory. This strategy is unavailable in a real application, but the experiment serves as an indication of the effectiveness of inclusion of a prior in this example. Future work is targeting the incremental estimation of these priors.
































§ DISCUSSION
















This article presents methods for estimation of linear, time-varying models. The methods presented extend directly to nonlinear models that remain linear in the parameters. When estimating an LTV model from a trajectory obtained from a nonlinear system, one is effectively estimating the linearization of the system around that trajectory. A first-order approximation to a nonlinear system is not guaranteed to generalize well as deviations from the trajectory become large. Many non-linear systems are, however, approximately locally linear, such that they are well described by a linear model in a small neighborhood around the linearization/operating point. For certain methods, such as iterative learning control and trajectory centric reinforcement learning, a first-order approximation to the dynamics is used for efficient optimization, while the validity of the approximation is ensured by incorporating penalties or constraints between two consecutive trajectories.

The methods presented allow very efficient learning of this first-order approximation due to the prior belief over the nature of the change in dynamics parameters, encoded by the regularization terms. By postulating a prior belief that the dynamics parameters change in a certain way, less demand is put on the data required for identification. The identification process will thus not interfere with normal operation in the same way as if excessive noise would be added to the input for identification purposes. This allows learning of flexible, over-parametrized models that fit available data well. This makes the proposed identification methods attractive in applications such as guided policy search (GPS) <cit.> and non-linear iterative learning control (ILC) <cit.>, where they can lead to dramatically decreased sample complexity.


When faced with a system where time-varying dynamics is suspected and no particular
knowledge regarding the dynamics evolution is available, or when the dynamics are known to vary slowly,
a reasonable first choice of algorithm is eq:smooth. It is also by
far the fastest of the proposed methods due to the Kalman-filter implementation
of <ref>.[The Kalman-filter implementation is often
several orders of magnitude faster than solving the optimization problems with
an iterative solver.] Example use cases include when dynamics are changing with a
continuous auxiliary variable, such as temperature, altitude or velocity.
If a smooth parameter drift is found to correlate with an auxiliary variable,
LPV-methodology can be employed to model the dependence explicitly.

Dynamics may change abruptly as a result of, e.g., system failure, change of operating mode,
or when a sudden disturbance enters the system, such as a policy change affecting a market or
a window opening affecting the indoor temperature. The identification method eq:pwconstant can be employed
to identify when such changes occur, without specifying a priori how many changes are expected.




For simplicity, the regularization weights were kept as simple scalars in this article.
However, all terms λΔ k^2 = (Δ k) (λ I) (Δ k)
can be generalized to (Δ k)Λ (Δ k), where Λ is an arbitrary positive definite matrix.
This allows incorporation of different scales for different variables with little added implementation complexity.




§ CONCLUSIONS

We have proposed a framework for identification of linear, time-varying models along trajectories of nonlinear dynamical systems using convex optimization. We showed how a Kalman smoother can be used to estimate the dynamics efficiently in a few special cases, and demonstrated the use of the proposed LTV models on three examples, highlighting their efficiency for trajectory-centric, model-based reinforcement learning, iterative learning control (ILC), and jump-linear system identification. We have also demonstrated the ability of the models to handle non-smooth friction dynamics as well as analyzed the identifiability of the models.









IEEEtran
