# Direct identification of continuous-time linear switched state-space models

Manas Mejari  Dario Piga  IDSIA Dalle Molle Institute for Artificial
Intelligence, USI-SUPSI, Via la Santa 1, CH-6962 Lugano-Viganello,
Switzerland. (e-mail: {manas.mejari, dario.piga}@supsi.ch).

###### Abstract

This paper presents an algorithm for direct _continuous-time_ (CT)
identification of _linear switched state-space_ (LSS) models. The key idea for
direct CT identification is based on an integral architecture consisting of an
LSS model followed by an integral block. This architecture is used to
approximate the continuous-time state map of a switched system. A properly
constructed objective criterion is proposed based on the integral architecture
in order to estimate the unknown parameters and signals of the LSS model. A
coordinate descent algorithm is employed to optimize this objective, which
alternates between computing the unknown model matrices, switching sequence
and estimating the state variables. The effectiveness of the proposed
algorithm is shown via a simulation case study.

###### keywords:

Continuous-time system estimation, Hybrid and switched systems modeling.

††thanks: This work has been submitted to IFAC World Congress’23 for possible
publication.

##  1 Introduction

###  1.1 Linear switched systems

Switched linear models belong to a class of hybrid systems, which consists of
multiple linear subsystems and a switching signal dictating the _active_
linear subsystem at each time instance. Such model class is widely used to
describe the behavior of dynamical systems subject to abrupt changes,
exhibiting both continuous and discrete dynamics. These changes can occur, for
instance, due to sensor/actuator failures, external disturbances or a change
in the operating point of a non-linear system.

Over the past few decades, switched linear models have found several
applications in a variety of fields including, mobile communication (Abdollahi
and Khorasani, 2011), signal processing (Doucet et al., 2001), computer vision
and bio-tracking (Oh et al., 2008), energy disaggregation (Mejari et al.,
2018), modeling human motion dynamics (Pavlovic et al., 2000), among many
others.

###  1.2 On direct continuous-time identification

Concerning the identification of linear switched systems, majority of the
approaches proposed in the literature have been developed for the
identification of _discrete-time_ (DT) models. Among these, we mention
optimization based algorithms (Bako, 2011; Ohlsson and Ljung, 2013), recursive
clustering-based approaches (Breschi et al., 2016; Mejari et al., 2020a),
mixed-integer programming algorithms (Mejari et al., 2020b), Bayesian
inference (Piga et al., 2020), algebraic-geometric approach (Vidal, 2008),
which identify DT switched linear models in _input-output_ (IO) form. Although
IO models are able to describe the behavior of the underlying system, often it
is desirable to estimate state-space representations, as they are more
convenient for stability analysis and controller synthesis of multi-input
multi-output plants. To this end, realization theory and subspace based
algorithms have been developed for DT switched linear state-space models, see
(Bako et al., 2009; Petreczky et al., 2013; Verdult and Verhaegen, 2004).

In comparison with the large number of contributions dedicated to the
estimation of DT switched models, very few works have addressed direct
_continuous-time_ (CT) identification of switched models. However, as
discussed in (Garnier, 2015; Garnier and Wang, 2008; Piga, 2018) for _linear
time-invariant_ (LTI) models, direct identification of CT model from sampled
data offers multiple advantages over the discrete-time case. Most of the
physical systems are naturally modelled in continuous-time, and thus, the
estimated parameters of CT models usually have a physical interpretation.
Direct CT identification methods can also deal with non-uniformly sampled
data, while discrete-time models implicitly rely on a fixed sampling time.
Moreover, CT identification methods are generally more robust to numerical
issues that may arise when using discrete-time methods in the case of high-
frequency sampled data.

Motivated by these advantages, our goal in this paper is to develop an
algorithm for direct CT identification for switched linear models. The core
idea is based on the concept of _integral architecture_ , recently introduced
by the authors in Mavkov et al. (2020); Mejari et al. (2022) for
identification of CT non-linear and LPV systems. In this work, we extend that
methodology for _linear switched state-space_ (LSS) model class.

###  1.3 Paper contributions and related works

We consider the problem of direct CT identification of LSS models which
involves: estimating the matrices of each LTI submodel, computing a discrete
mode sequence which indicates the active submodel at a given time, and
estimating the continuous state-sequence, from a given sampled IO data. The
proposed solution is based on an integral architecture consisting of LSS model
followed by an integral block, which is used to approximate the continuous
state dynamics of an LSS system. A block coordinate descent algorithm is
employed to optimize a properly constructed dual-objective criterion, which
alternates between computing the unknown matrices, discrete mode sequence and
estimating the states.

To the best of our knowledge, direct CT identification of switched linear
models has been addressed very recently only in (Goudjil et al., 2020;
Kersting and Buss, 2019; Du et al., 2021). These approaches, however, rely on
strong assumptions imposed on the system’s signals. In particular, the CT
identification method proposed in Goudjil et al. (2020) requires that the
input signal exciting the system is sinusoidal. Then, by exploiting the
linearity of subsystems, outputs of individual subsystems are estimated using
a DT switched IO method. In the next stage, CT identification approaches
developed for LTI models are employed to estimate model parameters based on
estimated subsystem outputs. In Kersting and Buss (2019), parameter
identifiers and concurrent learning is proposed based on the assumptions that
discrete mode sequence as well as continuous states are measured. Integral
concurrent learning is proposed in Du et al. (2021) relaxing the assumption of
known switching sequence. However, the continuous state is assumed to be
measured. To position our work w.r.t. these contributions, we do not impose
any of the aforementioned assumptions required in (Goudjil et al., 2020;
Kersting and Buss, 2019; Du et al., 2021), which are quite restrictive in
practice. In particular, in our contribution the input signals used to excite
the system are not restricted to sinusoidal inputs. On the contrary, any class
of input signals which excite all modes of the system can be used.
Furthermore, neither the discrete mode sequence nor the continuous state are
assumed to be known. The proposed algorithm estimates both these signals along
with the model matrices through a block coordinate-descent approach tailored
to the considered identification problem.

The paper is organized as follows. The identification problem for LSS models
is formalized in Section 2. The description of the integral architecture and
details of proposed proposed identification algorithm are provided in Section
3. A simulation example is reported in Section 4.

##  2 PROBLEM FORMULATION

We consider a data-generating system $\mathcal{S}$ governed by the following
CT linear switched state-space representation:

$\displaystyle\dot{\bf{x}}(t)$ $\displaystyle=A_{{\bf
s}(t)}{\bf{x}}(t)+B_{{\bf s}(t)}{{\bf u}}(t),$ (1a) $\displaystyle{\bf{x}}(0)$
$\displaystyle={\bf{x}}_{0},$ (1b) $\displaystyle{{\bf y}^{\rm o}}(t)$
$\displaystyle=C_{{\bf s}(t)}{\bf{x}}(t)+D_{{\bf s}(t)}{\bf u}(t),$ (1c)

where ${\bf{x}}(t)\in\mathbb{R}^{n_{x}}$ and
$\dot{\bf{x}}(t)\in\mathbb{R}^{n_{x}}$ are the state vector and its time
derivative, respectively; ${\bf{x}}_{0}\in\mathbb{R}^{n_{x}}$ is the initial
condition; ${\bf{u}}(t)\in\mathbb{R}^{n_{u}}$ is the system input; ${\bf
s}(t)\in\\{1,\ldots,K\\}$ is the switching signal and ${{\bf y}^{\rm
o}}(t)\in\mathbb{R}^{n_{y}}$ is the (noise-free) system output at time
$t\in\mathbb{R}$. Note that, the system is assumed to have $K$ operating
_modes_ , each corresponds to an LTI state-space system with real-valued
matrices $\\{A_{i},B_{i},C_{i},D_{i}\\}_{i=1}^{K}$ of appropriate dimensions.

A training dataset $\mathcal{D}$ of length $N$ is gathered from the linear
switched system $\mathcal{S}$ defined in (1) at time instants
$\\{t_{0}=0,t_{1},\dots,\;t_{N-1}\\}$. The dataset consists of input and
_noisy_ output samples: $\mathcal{D}=\\{{\bf u}(t_{k}),{\bf
y}(t_{k})\\}_{k=0}^{N-1}$ with sampling time $\Delta t$. The measured output
is corrupted by a zero-mean white Gaussian noise
$\eta\sim\mathcal{N}(0,\sigma^{2}_{\eta}I)$, _i.e._ , ${\bf y}(t_{k})={\bf
y}^{\rm o}(t_{k})+\eta(t_{k})$.

###### Problem 1

Given a training dataset $\mathcal{D}$, our goal is to identify a continuous-
time linear switched state-space (LSS) model, such that the model output
matches closely with the measured system output ${\bf y}(t)$.

The following assumptions are imposed on the system and signals in (1):

  1. 1.

input signal ${\bf u}(t)$ can be reconstructed (or reasonably approximated)
for all time instants $t\in[0\;\;t_{N-1}]\subset\mathbb{R}$ from the measured
samples $\\{{\bf u}(t_{k})\\}_{k=0}^{N-1}$.

  2. 2.

time evolution of the switching signal ${\bf s}(t)$ is stochastic, but it’s
value does not change during the sampling interval $\Delta t$.

  3. 3.

all $K$ operating modes are distinguishable and sufficiently excited.

In the paper, we consider a constant sampling time $\Delta t$ only to ease the
notation. Nevertheless, the approach is valid also for varying sampling times.

In the next paragraph, we introduce an _integral architecture_ for the
identification of continuous-time LSS model.

##  3 CONTINUOUS-TIME IDENTIFICATION OF LSS MODELS

###  3.1 Integral architecture

In order to describe the continuous-time state dynamics in (1a), we define an
LSS block $\mathcal{M}_{x}(\hat{\bf{x}},{\bf u},\hat{\bf s})$, which is fed by
the system input ${\bf u}(t)$, (estimated) switching signal $\hat{\bf s}(t)$
and (estimated) state $\hat{\bf{x}}(t)$ at time $t$, and returns the estimated
state time-derivative $\dot{\hat{{\bf{x}}}}(t)$, _i.e._ ,

$\displaystyle\mathcal{M}_{x}(\hat{\bf{x}},{\bf u},\hat{\bf s}):\ \
\dot{\hat{{\bf{x}}}}(t)\\!=\\!\hat{A}_{\hat{\bf
s}(t)}\hat{\bf{x}}(t)+\hat{B}_{\hat{\bf s}(t)}{\bf u}(t),$ (2)

where $\hat{A}_{i}\in\mathbb{R}^{n_{x}\times n_{x}}$ and
$\hat{B}_{i}\in\mathbb{R}^{n_{x}\times n_{u}}$ (for $i=1,\ldots,K$) are the
model matrices to be identified. Similarly, the output equation in (1c) is
represented by another block $\mathcal{M}_{y}(\hat{\bf{x}},{\bf u},\hat{\bf
s})$, which is fed by the estimated state $\hat{\bf{x}}(t)$, input ${\bf
u}(t)$ and estimated switching signal $\hat{\bf s}(t)$ and it returns the
model output $\hat{\bf{y}}(t)$ at time $t$, _i.e._ ,

$\displaystyle\mathcal{M}_{y}(\hat{\bf{x}},{\bf u},\hat{\bf
s}):\hat{\bf{y}}(t)\\!=\hat{C}_{\hat{\bf
s}(t)}\hat{\bf{x}}(t)+\hat{D}_{\hat{\bf s}(t)}{\bf u}(t),$ (3)

where the matrices $\hat{C}_{i}\in\mathbb{R}^{n_{y}\times n_{x}}$ and
$\hat{D}_{i}\in\mathbb{R}^{n_{y}\times n_{u}}$ (for $i=1,\ldots,K$) have to be
estimated from data.

For brevity, we introduce the following notation:
$\Theta_{x}=\\{\hat{A}_{i},\hat{B}_{i}\\}_{i=1}^{K}$,
$\Theta_{y}=\\{\hat{C}_{i},\hat{D}_{i}\\}_{i=1}^{K}$ and the mode sequence
${\bf S}=\left(\hat{\bf s}(0),\ldots,\hat{\bf s}(t_{N-1})\right)$

The resulting continuous-time LSS model is then given by:

$\displaystyle\dot{\hat{{\bf{x}}}}(t)$
$\displaystyle=\mathcal{M}_{x}(\hat{\bf{x}}(t),{\bf u}(t),\hat{\bf
s}(t);\Theta_{x}(\hat{\bf s}(t))),$ (4a) $\displaystyle\hat{\bf{x}}(0)$
$\displaystyle=\hat{\bf{x}}_{0},$ (4b) $\displaystyle{\hat{\bf{y}}}(t)$
$\displaystyle=\mathcal{M}_{y}(\hat{\bf{x}}(t),{\bf u}(t),\hat{\bf
s}(t);\Theta_{y}(\hat{\bf s}(t))),$ (4c)

where $\Theta_{x}(\hat{\bf s}(t))$ and $\Theta_{y}(\hat{\bf s}(t))$ denote the
matrices $(\hat{A}_{\hat{\bf s}(t)},$ $\hat{B}_{\hat{\bf s}(t)})$ and
$(\hat{C}_{\hat{\bf s}(t)},\hat{D}_{\hat{\bf s}(t)})$ respectively,
corresponding to the _active mode_ $\hat{\bf s}(t)$ at time $t$. Note that,
for the $i$-the active mode, _i.e.,_ $\hat{\bf s}(t)=i$ at time $t$, the maps
$\mathcal{M}_{x}(\cdot)$ and $\mathcal{M}_{y}(\cdot)$ are _linear_ functions
of the model matrices $(A_{i},B_{i})$ and $(C_{i},D_{i})$ respectively.

In this paper, we adopt a method originally introduced in Mavkov et al.
(2020), which exploits the _integral form_ of the Cauchy problem (4a)-(4b), by
defining an _integral_ block $\mathcal{M}_{I}$ as:

$\displaystyle\hat{\bf{x}}_{I}(t)=\mathcal{M}_{I}(\hat{\bf{x}}(t),{\bf
u}(t),\hat{\bf s}(t);\Theta_{x}(\hat{\bf s}(t)))$ (5)

with

$\displaystyle\mathcal{M}_{I}(\hat{\bf{x}}(t),{\bf u}(t),\hat{\bf
s}(t);\Theta_{x}(\hat{\bf s}(t)))$
$\displaystyle=\hat{\bf{x}}(0)+\int_{0}^{t}\mathcal{M}_{x}(\hat{\bf{x}}(\tau),{\bf
u}(\tau),\hat{\bf s}(\tau);\Theta_{x}(\hat{\bf s}(\tau)))d\tau.$

The block diagram in Fig. 1 is a representation of (5), along with the output
equation (4c) producing ${\hat{\bf{y}}}(t)$. Given the mode sequence ${\bf
S}$, if the state $\hat{\bf{x}}(t)$ feeding the LSS model block
$\mathcal{M}_{I}(\cdot)$ is actually generated by the model given in (4), then
the state $\hat{\bf{x}}_{I}(t)$ exactly matches $\hat{\bf{x}}(t)$, i.e.,

$\hat{\bf{x}}(t)=\hat{\bf{x}}_{I}(t)\qquad\forall t\in[t_{0}\;t_{N-1}].$ (6)

$\int\limits_{0}^{t}\dot{\hat{{\bf{x}}}}(\tau)d\tau$$\mathcal{M}_{x}\left(\Theta_{x}\right)$$\mathcal{M}_{y}\left(\Theta_{y}\right)$${\bf{x}}_{0}$$\hat{\bf{x}}_{I}(t)$$\hat{\bf{x}}(t)$${\bf
u}(t)$$\hat{\bf
s}(t)$$\dot{\hat{{\bf{x}}}}(t)$$\mathcal{M}_{I}(\cdot,\cdot)$$\hat{\bf{y}}(t)$$+$$+$
Figure 1: Integral architecture for continuous-time LSS model identification.

###  3.2 Fitting criterion

In the proposed scheme, the LSS model matrices $\Theta_{x},\Theta_{y}$, the
mode sequence ${\bf S}$ and the state signal
$\hat{\bf{x}}(t),t\in[t_{0}\;t_{N-1}]$, are _free_ optimization parameters.
They are concurrently optimized according to a cost function constructed with
the following rationale.

First, the estimated model output $\hat{\bf{y}}$ should match the output
measurements in training dataset $\mathcal{D}$. This objective is achieved by
introducing a _fitting term_ $J_{\bf y}$ in the cost function which penalizes
the mismatch between the model outputs $\hat{\bf{y}}(t_{k})$ and the sampled
measured outputs ${\bf y}(t_{k}),\;k=0,1,\dots,N\\!-\\!1$.

Second, the state signal $\hat{\bf{x}}$ should be compatible with the LSS
model dynamics (4). This can be achieved through an additional _regularization
term_ $J_{\bf x}$, which penalizes the distance between $\hat{\bf{x}}_{I}(t)$
and $\hat{\bf{x}}(t)$, where $\hat{\bf{x}}_{I}$ is defined as in (5). The
regularization term enforces the state $\hat{\bf{x}}$ (which is to be
optimized) to follow the CT model dynamics (4).

Furthermore, a _mode loss_ term
$\mathcal{L}:\mathcal{K}^{N}\rightarrow\mathbb{R}$ (where
$\mathcal{K}=\\{1,\ldots,K\\}$) is imposed on the mode sequence ${\bf S}$ to
take into account the temporal order as well as to incorporate the information
of the switching mechanism, for _e.g._ , Markovian switching in jump Markov
linear systems etc.

The following minimization problem is thus formulated:

$\begin{split}\underset{\hat{\bf{x}}(\cdot),{\bf
S},\Theta_{x},\Theta_{y}}{\text{min}}J(\hat{\bf{x}}(\cdot),{\bf
S},\Theta_{x},\Theta_{y}),\end{split}$ (7a) where
$\begin{split}J&=\underset{J_{\bf
y}}{\underbrace{\sum_{k=0}^{N-1}\lVert{\hat{\bf y}}(t_{k})-{\bf
y}(t_{k})\rVert^{2}}}\\\ &+\alpha\underset{J_{\bf
x}}{\underbrace{\int_{t_{0}}^{t_{N-1}}\lVert\hat{\bf{x}}_{I}(\tau)-\hat{\bf{x}}(\tau)\rVert^{2}\;d\tau}}+\beta\
\mathcal{L}({\bf S}),\end{split}$ (7b) with $\displaystyle\hat{\bf{y}}(t_{k})$
$\displaystyle=\mathcal{M}_{y}(\hat{\bf{x}}(t_{k}),{\bf u}(t_{k}),\hat{\bf
s}(t_{k});\;\Theta_{y}(\hat{\bf s}(t_{k}))),$ (7c)
$\displaystyle\hat{\bf{x}}_{I}(t)$
$\displaystyle=\hat{\bf{x}}(0)+\int_{0}^{t}\mathcal{M}_{x}(\hat{\bf{x}}(\tau),{\bf
u}(\tau),\hat{\bf s}(\tau);\;\Theta_{x}(\hat{\bf s}(t_{k})))\;d\tau.$ (7d)

As introduced in (Bemporad et al., 2018), the following mode loss
$\mathcal{L}:\mathcal{K}^{N}\rightarrow\mathbb{R}$ can be considered

$\displaystyle\mathcal{L}({\bf S})$ $\displaystyle=\mathcal{L}^{\rm
init}(\hat{\bf s}(0))+\sum_{k=1}^{N-1}\mathcal{L}^{\rm mode}(\hat{\bf
s}(t_{k}))$ $\displaystyle+\sum_{k=1}^{N-1}\mathcal{L}^{\rm trans}(\hat{\bf
s}(t_{k}),\hat{\bf s}(t_{k-1})),$ (8)

where $\mathcal{K}=\\{1,\ldots,K\\}$, $\mathcal{L}^{\rm
init}:\mathcal{K}\rightarrow\mathbb{R}$ is the initial mode cost,
$\mathcal{L}^{\rm mode}:\mathcal{K}\rightarrow\mathbb{R}$ is the mode cost,
and $\mathcal{L}^{\rm trans}:\mathcal{K}^{2}\rightarrow\mathbb{R}$ is the mode
transition cost.

The hyper-parameters $\alpha,\beta>0$ act as a tuning knob balancing the
relative importance of the fitting cost $J_{y}$, the regularization cost
$J_{x}$ and the _mode loss_ $\mathcal{L}({\bf S})$. Additionally, in order to
enforce smoothness properties for the estimated state variables, an $\ell_{1}$
regularization term $\|\hat{\bf{x}}(t_{k})-\hat{\bf{x}}(t_{k-1})\|_{1}$ can be
also included in the optimization problem (7a).

###  3.3 Integral approximation

Note that the continuous-time state signal
$\hat{\bf{x}}(t)\in\mathbb{R}^{n_{x}},\;t\in[t_{0}\;t_{N-1}]$ is one of the
problem’s decision variables. Indeed, the optimization problem (7) is
_infinite-dimensional_ and thus computationally intractable. Following the
rationale in (Mejari et al., 2022), we employ numerical techniques to
approximate (7) into a finite-dimensional problem amenable for a tractable
implementation. In particular, the state signal $\hat{\bf{x}}(t)$ is
approximated using a finite-dimensional parameterization. For simplicity of
exposition, we represent the state signal with a _piecewise constant_
parameterization, where $\hat{\bf{x}}(t)$ is constant during the sampling
intervals $[t_{k-1}\ t_{k}]$, $k=0,1,\dots,N-1$. In general, more complex
parametrizations for $\hat{\bf{x}}$ such as piecewise linear or polynomial
could be also used. Moreover, the intervals for the piecewise constant
approximation of $\hat{\bf{x}}$ may not necessarily correspond to the
input/output sampling time $\Delta t$ or the _dwell time_ of the switching
signal.

Furthermore, we approximate the integrals in (7b) and (7d) by applying a
numerical integration scheme. For simplicity, in this work we apply the
classical _rectangular approximation_ rule for the numerical integration of
(7b) and (7d). Other quadrature rules such as trapezoidal or Gaussian
quadrature could be alternatively considered.

Overall, the _piecewise constant_ parametrization of the signals
$\hat{\bf{x}}(t)$, ${\bf u}(t)$, $\hat{\bf s}(t)$ with the _rectangular
quadrature_ of the integrals leads to the following approximation:

$\displaystyle\int_{t_{0}}^{t_{N-1}}\\!\\!\lVert\hat{\bf{x}}_{I}(\tau)-\hat{\bf{x}}(\tau)\rVert^{2}d\tau\\!\approx\\!\sum_{k=1}^{N-1}\lVert\hat{\bf{x}}_{I}(t_{k})-\hat{\bf{x}}(t_{k})\rVert^{2}\Delta
t_{k},$

where $\Delta t_{k}=t_{k}-t_{k-1}$, and (7d) can be approximated with the
following Riemann sum:

$\displaystyle\hat{\bf{x}}_{I}(t_{k})$
$\displaystyle\approx\\!\hat{\bf{x}}(0)\\!+\\!\sum_{j=0}^{k-1}\Delta
t_{j+1}\mathcal{M}_{x}(\hat{\bf{x}}(t_{j}),{\bf u}(t_{j}),\hat{\bf
s}(t_{j});\;\Theta_{x}(\hat{\bf s}(t_{j}))).$ (9)
$\displaystyle=\hat{\bf{x}}(0)+\sum_{j=0}^{k-1}\Delta
t_{j+1}\left(\hat{A}_{\hat{\bf s}(t_{j})}\hat{\bf{x}}(t_{j})+\hat{B}_{\hat{\bf
s}(t_{j})}{\bf u}(t_{j})\right)$ (10)

The sum in the equation above can be also constructed recursively as follows:

$\displaystyle\hat{\bf{x}}_{I}(t_{k}\\!+\\!1)\\!=\\!\hat{\bf{x}}_{I}(t_{k})+\overbrace{\Delta
t_{k+1}\left(\hat{A}_{\hat{\bf s}(t_{k})}\hat{\bf{x}}(t_{k})+\hat{B}_{\hat{\bf
s}(t_{k})}{\bf u}(t_{k})\right)}^{\Delta{\bf{x}}_{k}}.$ (11)

###  3.4 Optimization algorithm

In the following, we report a numerical optimization algorithm in order to
minimize the cost function $J$ in (7a) w.r.t. the parameters
$\\{\hat{\bf{x}},{\bf S},\Theta_{x},\Theta_{y}\\}$. To this end, we employ the
_coordinate-descent_ approach as described in Algorithm 1. With a slight abuse
of notation, the optimization variable $\hat{\bf{x}}$ in Algorithm 1 denotes
the finite-dimensional representation of the state signal $\hat{\bf{x}}$,
i.e., $\hat{\bf{x}}=\\{\hat{\bf{x}}(t_{0}),\ldots,\hat{\bf{x}}(t_{N-1})\\}$.

Algorithm 1 Block coordinate descent for the estimation of states
$\hat{\bf{x}}$, mode sequence ${\bf S}$ and model parameter matrices
$\Theta_{x},\Theta_{y}$.

Input: Training dataset $\mathcal{D}=\\{{\bf u}(t_{k}),{\bf
y}(t_{k})\\}_{k=0}^{N-1}$; initial guess $\hat{\bf{x}}^{(0)},{\bf S}^{(0)}$;
tuning parameter $\alpha,\beta$; tolerance $\epsilon$, maximum number of
iterations $n_{\mathrm{max}}$.

  1. 1.

Iterate for $n=1,\ldots$

    1. 3..1.

$\Theta_{x}^{(n)},\Theta_{y}^{(n)}\leftarrow\mathrm{arg}\underset{\Theta_{x},\Theta_{y}}{\text{min}}J(\hat{\bf{x}}^{(n-1)},{\bf
S}^{(n-1)},\Theta_{x},\Theta_{y})$

    2. 3..2.

${\bf S}^{(n)}\leftarrow\mathrm{arg}\underset{{\bf S}}{\text{min}}\
J(\hat{\bf{x}}^{(n-1)},{\bf S},\Theta_{x}^{(n)},\Theta_{y}^{(n)})$

    3. 3..3.

$\hat{\bf{x}}^{(n)}\leftarrow\mathrm{arg}\underset{\hat{\bf{x}}}{\text{min}}\
\ J(\hat{\bf{x}},{\bf S}^{(n)},\Theta_{x}^{(n)},\Theta_{y}^{(n)})$

  2. 2.

Until $\|J^{(n)}-J^{(n-1)}\|\leq\epsilon$ or $n=n_{\mathrm{max}}$

Output: Estimated $\\{\hat{\bf{x}}(t_{k})\\}_{k=0}^{N-1}$, ${\bf S}$ and
$\Theta_{x},\Theta_{y}$.

Given an initial guess $\hat{\bf{x}}^{(0)}$ of the state and ${\bf S}^{(0)}$
of the mode sequence, at each iteration $n\geq 1$, Algorithm 1 alternates
between three steps: Step 3.3..3..0., Step 3.3..3..0. and Step 3.3..3..0.. In
particular, at Step 3.3..3..0., model parameters $\Theta_{x},\Theta_{y}$ are
computed by solving (7a) for a _fixed_ state $\hat{\bf{x}}^{(n-1)}$ and mode
sequence ${\bf S}^{(n-1)}$ obtained at the iteration $(n-1)$. At Step
3.3..3..0., the mode sequence is estimated for _fixed_ states
$\hat{\bf{x}}^{(n-1)}$ and fixed model parameters $\Theta_{x}^{(n)}$ and
$\Theta_{y}^{(n)}$ obtained from Step 3.3..3..0. at the $n$-th iteration.
Subsequently, at Step 3.3..3..0., the state sequence $\hat{\bf{x}}^{(n)}$ is
estimated by minimizing the cost (7a) for _fixed_ model parameters
$\Theta_{x}^{(n)}$ and $\Theta_{y}^{(n)}$ obtained from Step 3.3..3..0. and
mode sequence ${\bf S}^{(n)}$ computed from Step 3.3..3..0.. The procedure
continues until a maximum number of iterations is reached, or a certain
convergence criterion is met (Step 3.).

###### Remark 1

Since the underlying optimization problem is non-convex, convergence of
Algorithm 1 to the global optimal is sensitive to the initial guesses for
states $\hat{\bf{x}}^{(0)}$ and mode sequence ${\bf S}^{(0)}$. A possible
choice to initialize the state sequence is to first identify a continuous-time
LTI state-space model and set $\hat{\bf{x}}^{(0)}$ to the states of the LTI
model with small additive perturbations, _i.e.,_
$\hat{\bf{x}}^{(0)}=\hat{\bf{x}}_{\mathrm{LTI}}+v_{x}$ where,
$v_{x}\sim\mathcal{N}(0,\sigma^{2}_{x})$ with variance $\sigma^{2}_{x}$ chosen
by the user. We remark that, in practice, Algorithm 1 can be run multiple
times with different initial conditions and then choosing the best model
parameters according to a figure of merit.

###### Remark 2

We stress that in Algorithm 1, Steps 3.3..3..0. and 3.3..3..0. can be solved
_analytically_ via ordinary least squares, while Step 3.3..3..0. is solved to
global optimality via dynamic programming. Thus, each sub-problem to be
optimized within an iteration of the block corordinte descent is solved
exactly to its unique optimal solution, which can be utilized to prove the
convergence guarantees, see (Razaviyayn et al., 2013; Tseng, 2001).

In the following section, we detail each step of the coordinate descent
algorithm. Without loss of generality, for brevity, we set
$\hat{C}_{i}=\hat{C}$, $\hat{D}_{i}=0$ for all $i=1,\ldots,K$.

####  3.4.1 Step 3.3..3..0.: Optimization over model parameters
$\Theta_{x},\Theta_{y}$

For a fixed mode sequence ${\bf S}$ and a fixed state sequence $\hat{\bf{x}}$,
the cost function $J(\hat{\bf{x}}^{(n-1)},{\bf
S}^{(n-1)},\Theta_{x},\Theta_{y})$ in (7b) can be optimized over the unknown
model parameters $\Theta_{x},\Theta_{y}$. This leads to a _least-squares_
problem described as follows.

Let $\theta_{i}=\left[\mathrm{vec}(\hat{A}_{i})^{\top}\
\mathrm{vec}(\hat{B}_{i})^{\top}\right]^{\top}\in\mathbb{R}^{n_{\theta}}$,
with $n_{\theta}=n_{x}(n_{x}+n_{u})$ and let us define the matrix
$\Phi(t_{j})\in\mathbb{R}^{n_{\theta}\times n_{x}}$ as follows

$\Phi(t_{j})=\Delta t\left[\begin{smallmatrix}\hat{\bf{x}}(t_{j})\\\ {\bf
u}(t_{j})\end{smallmatrix}\right]\otimes I_{n_{x}}$ (12)

with $\otimes$ denoting the Kronecker product.

The approximated state evolution eq. (9) can be written as,

$\displaystyle\hat{\bf{x}}_{I}(t_{k})$
$\displaystyle\approx\\!\hat{\bf{x}}(0)\\!+\\!\sum_{j=0}^{k-1}\Phi^{\top}(t_{j})\theta_{\hat{\bf
s}(t_{j})}$ (13)

Let us define the matrix $P\in\mathbb{R}^{(N-1)\times K}$ such that, for each
of its row $i=1,\ldots,N-1$, the $j$-th column is set to $1$ if the active
mode at time $t_{i}$ is $\hat{\bf s}(t_{i})=j$, _i.e._ , the $(i,j)$-th entry
$P_{i,j}$ is defined as

$\displaystyle P_{i,j}$ $\displaystyle=1\ \ \mathrm{if}\ \hat{\bf s}(t_{i})=j$
$\displaystyle=0\ \ \mathrm{otherwise}$

and let us define $\bar{P}=P\otimes
I_{n_{\theta}}\in\mathbb{R}^{(N-1)n_{\theta}\times Kn_{\theta}}$.

With the matrices defined above, the relation (13) can be written in the
matrix form as follows:

$\displaystyle\underset{\Delta\hat{\bf{x}}}{\underbrace{\left[\begin{smallmatrix}\hat{\bf{x}}_{I}(t_{1})-\hat{\bf{x}}(0)\\\
\hat{\bf{x}}_{I}(t_{2})-\hat{\bf{x}}(0)\\\ \vdots\\\
\hat{\bf{x}}_{I}(t_{N-1})-\hat{\bf{x}}(0)\end{smallmatrix}\right]}}=$ (14)
$\displaystyle\underset{\Psi}{\underbrace{\left[\begin{smallmatrix}\Phi^{\top}(t_{0})&0&0&\cdots&0\\\
\Phi^{\top}(t_{0})&\Phi^{\top}(t_{1})&0&\cdots&0\\\ \vdots&&&&\\\
\Phi^{\top}(t_{0})&\Phi^{\top}(t_{1})&\Phi^{\top}(t_{2})&\cdots&\Phi^{\top}(t_{N-2})\end{smallmatrix}\right]}}\bar{P}\left[\begin{smallmatrix}\theta_{1}\\\
\theta_{2}\\\ \vdots\\\ \theta_{K}\end{smallmatrix}\right]$ (15)

Based on the above definitions, the cost function (7b) can be re-written as

$\displaystyle J=\underset{J_{\bf
y}}{\underbrace{\left\|\tilde{C}\hat{\bf{x}}-{\bf
y}\right\|^{2}}}+\alpha\underset{J_{\bf
x}}{\underbrace{\left\|\Delta\hat{\bf{x}}-\Psi\bar{P}\Theta\right\|^{2}\;\Delta
t}}$ (16)

where $\tilde{C}=blk(\hat{C})$ is a block-diagonal matrix, $\hat{\bf{x}}$ and
${\bf y}$ are the sequences of estimated states and measured outputs
respectively. Note that, for a given state estimates $\hat{\bf{x}}$ and a
given mode sequence ${\bf S}$, the matrices $\Psi$ and $\bar{P}$ can be pre-
computed and thus, (16) is a _least-squares_ problem in the unknown model
parameters $\\{\theta_{i}\\}_{i=1}^{K},\tilde{C}$ (_i.e._ ,
$\\{\hat{A}_{i},\hat{B}_{i}\\}_{i=1}^{K},\hat{C}$), which can be solved
analytically.

####  3.4.2 Step 3.3..3..0.: Optimization over mode sequence ${\bf S}$

Given the estimates of the model parameters
$\\{\hat{A}_{i},\hat{B}_{i}\\}_{i=1}^{K},\hat{C}$ computed at Step 3.3..3..0.
and given a fixed state sequence $\hat{\bf{x}}$, the cost function
$J(\hat{\bf{x}}^{(n-1)},{\bf S},\Theta_{x}^{(n)},\Theta_{y}^{(n)})$ in (7b)
can be optimized over the unknown mode sequence ${\bf S}$ via discrete
_Dynamic Programming_ (DP) algorithm (Bemporad et al., 2018).

The DP algorithm to estimate the mode sequence ${\bf S}$ is summarized as
follows. Let $\ell(\hat{\bf{x}}(t_{k}),{\bf u}(t_{k}),\theta_{\hat{\bf
s}(t_{k})})$ be the transition cost defined as

$\displaystyle\ell(\hat{\bf{x}}(t_{k}),{\bf u}(t_{k}),\theta_{\hat{\bf
s}(t_{k})})$
$\displaystyle=\left\|\hat{\bf{x}}(t_{k+1})-\hat{\bf{x}}(t_{k})-\Delta
t\hat{A}_{\hat{\bf s}(t_{k})}\hat{\bf{x}}(t_{k})-\Delta t\hat{B}_{\hat{\bf
s}(t_{k})}{\bf u}(t_{k})\right\|^{2}$

We compute a matrix of cost $V\in\mathbb{R}^{K\times(N+1)}$ and a matrix of
indices $U\in\mathbb{R}^{K\times N}$ as follows:

First, the terminal cost $V_{i,N}\in\mathbb{R}$ is computed for all modes

$\displaystyle V_{i,N-1}=\mathcal{L}^{\rm
mode}(i)+\ell(\hat{\bf{x}}(t_{N-1}),{\bf u}(t_{N-1}),\theta_{i})\ \
i=1,\ldots,K.$

Next, the cost $V_{i,k}$ and indexes $U_{i,k}$ at time $t_{k}$ are computed
with following dynamic programming recursions, backwards in time for
$k=N-2,\ldots,1,0$

$\displaystyle U(i,k)$ $\displaystyle=\underset{j=1,\ldots,K}{\mathrm{arg}\
\mathrm{min}}\ \ \\{V_{j,k+1}+\mathcal{L}^{\rm trans}(j,i)\\},\quad
i=1,\ldots,K$ $\displaystyle V_{i,k}$ $\displaystyle=\mathcal{L}^{\rm
mode}(i)+\ell(\hat{\bf{x}}(t_{k}),{\bf u}(t_{k}),\theta_{i})+V_{U(i,k),k+1}$
$\displaystyle+\mathcal{L}^{\rm trans}(U(i,k),i),$ $\displaystyle V_{i,0}$
$\displaystyle=\mathcal{L}^{\rm
init}(i)+\underset{j=1,\ldots,K}{\mathrm{min}}\ \ \\{V_{j,1}+\mathcal{L}^{\rm
trans}(j,i)\\}$

The minimum cost mode sequence ${\bf S}$ is retrieved forward in time by
setting

$\displaystyle\hat{\bf s}_{t_{0}}$
$\displaystyle=\underset{j=1,\ldots,K}{\mathrm{arg}\ \mathrm{min}}\ \ V_{j,0}$
$\displaystyle\hat{\bf s}_{t_{k}}$ $\displaystyle=U(\hat{\bf s}_{t_{k-1}},k)\
\ \ k=1,\ldots,N-1$

####  3.4.3 Step 3.3..3..0.: Optimization over states $\hat{\bf{x}}$

By assuming that the initial conditions $\hat{\bf{x}}_{I}(t_{0})$ and
$\hat{\bf{x}}(t_{0})$ are equal, i.e.,

$\displaystyle\hat{\bf{x}}_{I}(t_{0})=\hat{\bf{x}}(t_{0}),$ (17)

the relation (11) can be written in the matrix form:

$\displaystyle\left[\begin{smallmatrix}\hat{\bf{x}}_{I}(t_{0})\\\
\hat{\bf{x}}_{I}(t_{1})\\\ \hat{\bf{x}}_{I}(t_{2})\\\ \vdots\\\
\hat{\bf{x}}_{I}(t_{N-1})\end{smallmatrix}\right]=$
$\displaystyle\underset{\tilde{A}}{\underbrace{\left[\begin{smallmatrix}I&0&0&\cdots&0&0\\\
I+\Delta t\hat{A}_{\hat{\bf s}(t_{0})}&0&0&\cdots&0&0\\\ I+\Delta
t\hat{A}_{\hat{\bf s}(t_{0})}&\Delta t\hat{A}_{\hat{\bf
s}(t_{1})}&0&\cdots&0&0\\\ \vdots&&&&&\\\ I+\Delta t\hat{A}_{\hat{\bf
s}(t_{0})}&\Delta t\hat{A}_{\hat{\bf s}(t_{1})}&\Delta t\hat{A}_{\hat{\bf
s}(t_{2})}&\cdots&\Delta t\hat{A}_{\hat{\bf
s}(t_{N-2})}&0\end{smallmatrix}\right]}}\left[\begin{smallmatrix}\hat{\bf{x}}(t_{0})\\\
\hat{\bf{x}}(t_{1})\\\ \hat{\bf{x}}(t_{2})\\\ \vdots\\\
\hat{\bf{x}}(t_{N-1})\end{smallmatrix}\right]$ $\displaystyle+$
$\displaystyle\underset{\tilde{B}}{\underbrace{\left[\begin{smallmatrix}0\\\
\Delta t\hat{B}_{\hat{\bf s}(t_{0})}{\bf u}(t_{0})\\\ \Delta
t\hat{B}_{\hat{\bf s}(t_{0})}{\bf u}(t_{0})+\Delta t\hat{B}_{\hat{\bf
s}(t_{1})}{\bf u}(t_{1})\\\ \vdots\\\ \Delta t\hat{B}_{\hat{\bf s}(t_{0})}{\bf
u}(t_{0})+\Delta t\hat{B}_{\hat{\bf s}(t_{1})}{\bf u}(t_{1})+\ldots+\Delta
t\hat{B}_{\hat{\bf s}(t_{N-2})}{\bf u}(t_{N-2})\end{smallmatrix}\right]}}$
(18)

Thus, based on the above approximation, the cost function (7b) can be also re-
written as

$\displaystyle J=\underset{J_{\bf
y}}{\underbrace{\left\|\tilde{C}\hat{\bf{x}}-{\bf
y}\right\|^{2}}}+\alpha\underset{J_{\bf
x}}{\underbrace{\left\|(\tilde{A}-I)\hat{\bf{x}}+\tilde{B}\right\|^{2}\;\Delta
t}}+\mathcal{L}({\bf S}),$ (19)

where $\tilde{C}=blk(C)$, which can be solved for $\hat{\bf{x}}$ via ordinary
least-squares.

##  4 SIMULATION EXAMPLE

The performance of the proposed algorithm is assessed via a simulation case
study. All computations are carried out on an i7 1.9-GHz Intel core processor
with 32 GB of RAM running MATLAB R2019a.

We consider a continuous-time linear switched system governed by (1) having
$K=2$ modes with subsystem matrices given as follows (Goudjil et al., 2020):

$\displaystyle\left[\begin{array}[]{c|c}A_{1}&B_{1}\\\ \hline\cr C_{1}&D_{1}\\\ \end{array}\right]=\left[\begin{array}[]{c c | c }0&-120&120\\\ 1&-4&-12\\\ \hline\cr 0&1&0\end{array}\right],$ $\displaystyle\left[\begin{array}[]{c|c}A_{2}&B_{2}\\\ \hline\cr C_{2}&D_{2}\\\ \end{array}\right]=\left[\begin{array}[]{c c | c }0&-50&53\\\ 1&-1.8&25\\\ \hline\cr 0&1&0\end{array}\right]$

The system is excited with a zero-mean Gaussian input signal having unit
variance, ${\bf u}(t)\sim\mathcal{N}(0,1)$. The dynamics switches between the
two subsystems with a Markov switching signal, such that the true mode ${\bf
s}_{t_{k}}$ has $\pi=10\%$ probability of being different from ${\bf
s}_{t_{k-1}}$, starting from ${\bf s}_{t_{0}}=1$. The system belongs to a
class of switched models termed as continuous-time jump Markov linear systems
(Costa et al., 2013). Training dataset of $N=400$ samples is gathered,
sampling the output and input trajectories with a sampling time of $\Delta
t=0.01$ s. The output is corrupted by an additive white Gaussian noise ${\bf
y}(t_{k})={\bf y}^{\rm o}(t_{k})+\eta(t_{k})$ where
$\eta(t_{k})\sim\mathcal{N}(0,\sigma_{\eta}^{2})$ with $\sigma_{\eta}=0.025$,
which corresponds to signal-to-noise ratio of $30$ dB.

For identification, we consider an LSS model structure (4) with state
dimension set to the true system dimension $n_{x}=2$ and number of modes set
to $K=2$. The model matrices $\Theta_{x}$, $\Theta_{y}$, the mode sequence
${\bf S}$ and the state sequence $\hat{\bf{x}}$ are estimated by running the
coordinate descent Algorithm 1 for $n_{\max}=1000$ iterations. The average
computational time for each iteration of the algorithm is $0.1$ s, which
includes the time to compute the model matrices and states via ordinary least-
squares and estimation of mode sequence via dynamic programming recursions. In
total, the entire identification problem is completed in about 500 s, with 5
different initial guesses.

To asses the convergence properties, cost $J$ is plotted in Fig. 2 against the
iterations of the coordinate descent algorithm.

Figure 2: Training: Cost function _vs_ number of iterations.

As mentioned in Remark 1, the initial guess $\hat{\bf{x}}^{(0)}$ for the state
sequence is set to the states of an identified CT LTI state-space model111The
CT LTI model is identified using MATLAB’s system identification toolbox with
command $\mathrm{ssest}$ which employs $\mathrm{N4SID}$ subspace algorithm.
with small additive perturbations, _i.e.,_
$\hat{\bf{x}}^{(0)}=\hat{\bf{x}}_{\mathrm{LTI}}+v_{x}$ where,
$v_{x}\sim\mathcal{N}(0,\sigma^{2}_{x}I)$ with $\sigma_{x}=0.01$. The initial
guess for the mode sequence ${\bf S}^{(0)}$ is chosen randomly.

For the mode loss $\mathcal{L}({\bf S})$ (see eq. (3.2)), we set initial mode
cost $\mathcal{L}^{\rm init}({\bf s}_{t_{0}})=0$, mode cost $\mathcal{L}^{\rm
mode}({\bf s}_{t_{k}})=0$. The transition mode cost $\mathcal{L}^{\rm trans}$
is chosen as follows:

$\mathcal{L}^{\rm trans}({\bf s}_{t_{k}},{\bf
s}_{t_{k-1}})=\left\\{\begin{matrix}-\tau\log(1-\pi)\quad\mathrm{if}\ {\bf
s}_{t_{k}}={\bf s}_{t_{k-1}}\\\ \ \tau\log(\pi)\quad\quad\quad\ \mathrm{if}\
{\bf s}_{t_{k}}\neq{\bf s}_{t_{k-1}}\end{matrix}\right.$ (20)

with $\pi=0.1$ and $\tau=10^{-6}$. The regularization hyper-parameter $\alpha$
is set to $0.01$. The hyper-parameters $\tau$ and $\alpha$ are chosen via a
grid search.

Table 1: True _vs_ estimated transfer functions for mode $1$ and mode $2$. Mode | True | Estimated
---|---|---
$G_{1}(s)$ | $\frac{-12s+120}{s^{2}+4s+120}$ | $\frac{-12.01s+122.3}{s^{2}+3.98s+119}$
$G_{2}(s)$ | $\frac{25s+53}{s^{2}+1.8s+50}$ | $\frac{24.97s+52.9}{s^{2}+1.89s+49.3}$

The true and estimated transfer functions $G(s)$ of the two linear subsystems
are reported in Table 1. The corresponding Bode plots of the subsystems are
depicted in Fig. 3. The obtained results show that the model parameters of
each subsystem have been identified with high accuracy and input-output
behavior of the estimated linear subsystems matches closely to that of the
true subsystems.

Figure 3: Bode plot: true (red) _vs_ estimated (dashed blue) model for mode
$1$ (top panel) and mode $2$ (bottom panel).

The performance of the proposed identification algorithm is further assessed
in terms of mode sequence estimation, quantified via a _mode fit_ (MF) index
$\frac{100}{N}\sum_{k=0}^{N-1}\delta(\hat{\bf s}_{t_{k}},{\bf s}_{t_{k}})$
where $\delta(i,j)$ is the Kronecker delta function. Fig. 4 shows the true
_vs_ the estimated mode sequence. Only $7$ out of $400$ modes have been
incorrectly classified which leads to a mode fit of $98.25\%$. It is clear
from Fig. 4 that, starting from a random initial mode sequence, the proposed
algorithm is able to reconstruct the true mode sequence accurately.

Figure 4: Estimation of mode sequence: Top panel: Random intial guess ${\bf
S}^{(0)}$, Bottom panel: True (red) _vs_ estimated (blue) mode sequence with
mode fit percent = $98.25\%$.

Furthermore, we asses the effect of varying noise levels $\sigma_{\eta}$ as
well as regularization hyper-parameters $\tau$ on the performance of the
proposed algorithm quantified via mode-fit index. The results are summarized
in Fig. 5, which shows the percentage of correctly classified modes against
the values of $\tau$ for different noise standard deviations
$\sigma_{\eta}=\\{0.02,0.03,0.05,0.08\\}$ corresponding to signal-to-noise
ratios $\\{32,28,24,20\\}$ dB, respectively. From the choice of the mode
transition loss $\mathcal{L}^{\rm trans}$ in (20), higher values of $\tau$
implies more penalty on the change of mode. In other words, for large values
of $\tau$, mode change is discouraged and only single mode is recognized
(typically, the value of the initial mode $\hat{\bf s}_{t_{0}}$ is retained),
leading to a lower mode fit percent as seen in Fig. 5. The hyperparameter
$\tau$ and the choice of $\mathcal{L}^{\rm trans}$ thus act as a tuning knob,
which can be chosen via cross-validation, depending upon either fast or slow
switching dynamics.

Figure 5: Mode fit percent _vs_ $\tau$ for different noise levels. Optimal
value of $\tau$ shown with dashed lines.

Finally, in order to analyze the statistical properties and robustness of the
proposed algorithm, we perform a Monte-Carlo (MC) analyses with $50$ MC runs.
At each MC run, data is gathered by exciting the system with a new realization
of the input, switching signal and noise. The variance of the noise
distribution is set such that the average SNR for each run is $30$ dB.
Algorithm 1 is run for $5$ different initial guesses for the states
$\hat{\bf{x}}^{(0)}$ and mode sequence ${\bf S}^{(0)}$, with
$n_{\mathrm{max}}=1000$ iterations setting $(\alpha,\tau)=(0.01,3\cdot
10^{-7})$. Among $5$ different initializations, model parameters obtained from
the run having maximum _best fit rate_ :
$\mathrm{BFR}=100\left(1-\sqrt{\frac{\sum_{k=0}^{N-1}(\hat{\bf{y}}(t_{k})-{\bf
y}(t_{k}))^{2}}{\sum_{k=0}^{N-1}({\bf y}(t_{k})-\bar{{\bf y}})^{2}}}\right)\%$
are chosen.

Figure 6: Monte-Carlo analysis: Best fit rate and mode fit percentage over
$50$ MC runs.

The box-plots of the mode fit (MF) and BFR indexes over $50$ Monte-Carlo runs
are shown in Fig. 6. We observe that satisfactory performance is obtained in
terms of reconstruction of the output as well as the mode sequence. We remark
that, although convergence is not guaranteed for every run of the algorithm,
in practice, running the algorithm with $5$ different initializations was
sufficient to achieve accurate model parameter estimates.

##  5 CONCLUSIONS

In this work, we have presented an integral architecture for continuous-time
identification of switched state-space models. The proposed approach can be
seen as the first step towards developing a generic framework for direct
identification of continuous-time state-space hybrid dynamical systems. The
presented analysis has shown that satisfactory results have been achieved for
identifying a Markov jump linear system, in terms of reconstruction of the
mode sequence as well as capturing the input-output behaviours of the linear
submodels. Future works will focus on developing refinement strategies in
order to improve the estimation of mode sequence and robustness w.r.t. to
initial conditions.

## References

  * Abdollahi and Khorasani (2011) Abdollahi, F. and Khorasani, K. (2011).  A decentralized markovian jump $\mathcal{H}_{\infty}$ control routing strategy for mobile multi-agent networked systems.  _IEEE Transactions on Control Systems Technology_ , 19(2), 269–283.
  * Bako (2011) Bako, L. (2011).  Identification of switched linear systems via sparse optimization.  _Automatica_ , 47(4), 668 – 677.
  * Bako et al. (2009) Bako, L., Mercère, G., Vidal, R., and Lecoeuche, S. (2009).  Identification of switched linear state space models without minimum dwell time.  In _Proc. 15th IFAC Symposium on System Identification_ , 569–574. Saint-Malo, France.
  * Bemporad et al. (2018) Bemporad, A., Breschi, V., Piga, D., and Boyd, S. (2018).  Fitting jump models.  _Automatica_ , 96, 11–21.
  * Breschi et al. (2016) Breschi, V., Piga, D., and Bemporad, A. (2016).  Piecewise affine regression via recursive multiple least squares and multicategory discrimination.  _Automatica_ , 73, 155–162.
  * Costa et al. (2013) Costa, O.L.V., Fragoso, M.D., and Todorov, M.G. (2013).  _Continuous-Time Markov Jump Linear Systems_.  Springer.
  * Doucet et al. (2001) Doucet, A., Gordon, N., and Krishnamurthy, V. (2001).  Particle filters for state estimation of jump markov linear systems.  _IEEE Transactions on Signal Processing_ , 49(3), 613–624.
  * Du et al. (2021) Du, Y., Liu, F., Qiu, J., and Buss, M. (2021).  Online identification of piecewise affine systems using integral concurrent learning.  _IEEE Transactions on Circuits and Systems I: Regular Papers_ , 68(10), 4324–4336.
  * Garnier (2015) Garnier, H. (2015).  Direct continuous-time approaches to system identification. overview and benefits for practical applications.  _European Journal of control_ , 24, 50–62.
  * Garnier and Wang (2008) Garnier, H. and Wang, L. (2008).  _Identification of Continuous-time Models from Sampled Data_.  Springer Publishing Company.
  * Goudjil et al. (2020) Goudjil, A., Pouliquen, M., Pigeon, E., Gehan, O., and Bonargent, T. (2020).  Continuous-time identification for a class of switched linear systems.  In _Proc. 2020 European Control Conference (ECC)_ , 521–526. Saint Petersburg, Russia.
  * Kersting and Buss (2019) Kersting, S. and Buss, M. (2019).  Recursive estimation in piecewise affine systems using parameter identifiers and concurrent learning.  _International Journal of Control_ , 92(6), 1264–1281.
  * Mavkov et al. (2020) Mavkov, B., Forgione, M., and Piga, D. (2020).  Integrated neural networks for nonlinear continuous-time system identification.  _IEEE Control Systems Letters_ , 4(4), 851–856.
  * Mejari et al. (2020a) Mejari, M., Breschi, V., and Piga, D. (2020a).  Recursive bias-correction method for identification of piecewise affine output-error models.  _IEEE Control Systems Letters_ , 4(4), 970–975.
  * Mejari et al. (2022) Mejari, M., Mavkov, B., Forgione, M., and Piga, D. (2022).  Direct identification of continuous-time LPV state-space models via an integral architecture.  _Automatica_ , 142, 110407.
  * Mejari et al. (2018) Mejari, M., Naik, V.V., Piga, D., and Bemporad, A. (2018).  Energy disaggregation using piecewise affine regression and binary quadratic programming.  In _Proc. 57th IEEE Conference on Decision and Control (CDC)_ , 3116–3121. Miami Beach, FL, USA.
  * Mejari et al. (2020b) Mejari, M., Naik, V.V., Piga, D., and Bemporad, A. (2020b).  Identification of hybrid and linear parameter-varying models via piecewise affine regression using mixed integer programming.  _International Journal of Robust and Nonlinear Control_ , 30(15), 5802–5819.
  * Oh et al. (2008) Oh, S.M., Rehg, J.M., Balch, T., and Dellaert, F. (2008).  Learning and inferring motion patterns using parametric segmental switching linear dynamic systems.  _International Journal of Computer Vision_ , 77, 103–124.
  * Ohlsson and Ljung (2013) Ohlsson, H. and Ljung, L. (2013).  Identification of switched linear regression models using sum-of-norms regularization.  _Automatica_ , 49(4), 1045–1050.
  * Pavlovic et al. (2000) Pavlovic, V., Rehg, J.M., and MacCormick, J. (2000).  Learning switching linear models of human motion.  In _Proc. of the 13th International Conference on Neural Information Processing Systems_ , 942–948. Denver CO, United States.
  * Petreczky et al. (2013) Petreczky, M., Bako, L., and van Schuppen, J.H. (2013).  Realization theory of discrete-time linear switched systems.  _Automatica_ , 49(11), 3337–3344.
  * Piga (2018) Piga, D. (2018).  Finite-horizon integration for continuous-time identification: bias analysis and application to variable stiffness actuators.  _International Journal of Control_ , 93.
  * Piga et al. (2020) Piga, D., Bemporad, A., and Benavoli, A. (2020).  Rao-Blackwellized sampling for batch and recursive Bayesian inference of Piecewise Affine models.  _Automatica_ , 117, 109002.
  * Razaviyayn et al. (2013) Razaviyayn, M., Hong, M., and Luo, Z. (2013).  A unified convergence analysis of block successive minimization methods for nonsmooth optimization.  _SIAM Journal on Optimization_ , 23(2), 1126–1153.
  * Tseng (2001) Tseng, P. (2001).  Convergence of a block coordinate descent method for nondifferentiable minimization.  _Journal of Optimization Theory and Applications_ , 109, 475–494.
  * Verdult and Verhaegen (2004) Verdult, V. and Verhaegen, M. (2004).  Subspace identification of piecewise linear systems.  In _Proc. 43rd IEEE Conference on Decision and Control (CDC)_ , volume 4, 3838–3843. Atlantis, Bahamas.
  * Vidal (2008) Vidal, R. (2008).  Recursive identification of switched ARX systems.  _Automatica_ , 44(9), 2274–2287.