Offline state estimation for hybrid systems via nonsmooth variable projection

TL;DR

This paper introduces an offline method for estimating both discrete and continuous states of hybrid systems by formulating a continuous optimization problem and solving it with a novel nonsmooth variable projection algorithm, demonstrating effectiveness on impact systems.

Contribution

It presents a new offline estimation algorithm for hybrid systems that relaxes discrete states and employs a nonsmooth optimization approach with convergence guarantees.

Findings

Effective estimation on impact systems

Algorithm converges to stationary points

Handles large state changes during switching

Abstract

A hybrid dynamical system switches between dynamic regimes at time- or state-triggered events. We propose an offline algorithm that simultaneously estimates discrete and continuous components of a hybrid system's state. We formulate state estimation as a continuous optimization problem by relaxing the discrete component and use a robust loss function to accommodate large changes in the continuous component during switching events. Subsequently, we develop a novel nonsmooth variable projection algorithm with Gauss-Newton updates to solve the state estimation problem and prove the algorithm's global convergence to stationary points. We demonstrate the effectiveness of our approach on simple piecewise-linear and -nonlinear mechanical systems undergoing intermittent impact.

Tables

Table 1. Table 1 : Discrete states and continuous dynamics for impact oscillator hybrid system models (Sec.  4.1 ). Note that the continuous dynamics q ¨ ¨ 𝑞 \ddot{q} have the same general form for both the piecewise-linear and -nonlinear models, with the spring law k 𝑘 k being a linear or nonlinear function of the continuous state x = ( q , q ˙ ) 𝑥 𝑞 ˙ 𝑞 x=(q,\dot{q}) depending on which model is considered.

Discrete state $w$	Icon	${\ddot{q}}_w\left(x\right)$
$w=\text{A}\downarrow$		$\left[\begin{array}{c}\frac{1}{m_h}\left(-k_1(q,\dot{q})\right)-g\\ \frac{1}{m_t}\left(k_1(q,\dot{q})\right)-g\end{array}\right]$
$w=\text{G}\downarrow$		$\left[\begin{array}{c}\frac{1}{m_h}\left(-k_1(q,\dot{q})\right)-g\\ 0\end{array}\right]$
$w=\text{G}\uparrow$		$\left[\begin{array}{c}\frac{1}{m_h}\left(-k_2(q,\dot{q})\right)-g\\ 0\end{array}\right]$
$w=\text{A}\uparrow$		$\left[\begin{array}{c}\frac{1}{m_h}\left(-k_2(q,\dot{q})\right)-g\\ \frac{1}{m_t}\left(k_2(q,\dot{q})\right)-g\end{array}\right]$

Equations

$$x_{t+1}$$

$$y_t$$

$$r\log \left(r+{\Vert Q^{-1/2}{\sigma }_t\Vert }^2\right)-r\log (r),$$

$$\frac{1}{2}{\Vert R^{-1/2}{\delta }_t\Vert }^2,$$

$$\underset{x_t,w_t\in {\mathcal{D}}^m}{\min }\sum _{t=0}^{T-1}\frac{1}{2}{\Vert R^{-1/2}\left(y_t-{\mathcal{H}}_t(x_t)\right)\Vert }^2+$$

$$r\log \left(r+{\Vert Q^{-1/2}\left(x_{t+1}-\sum _{m=1}^M{\mathcal{F}}_m(x_t)w_t[m]\right)\Vert }^2\right).$$

$$\underset{w\in {\mathcal{D}}^M}{\min }g\left(x_2-\sum _{m=1}^{M}w[m]{\mathcal{F}}_m(x_1)\right)$$

$$=$$

$$\underset{w\in {\mathcal{D}}^M}{\mathrm{argmin}}g\left(x_2-\sum _{m=1}^{M}w[m]{\mathcal{F}}_m(x_1)\right)$$

$$=$$

$$g(x_2-{\mathcal{F}}_1(x_1)),g(x_2-{\mathcal{F}}_2(x_1)),\dots ,g(x_2-{\mathcal{F}}_m(x_1)).$$

$$\displaystyle\min_{x_{t},w_{t}\in\mathcal{D}^{M}}\sum_{t=0}^{T-1}\bigg{(}\frac{1}{2}\left\|R^{-1/2}\left(y_{t}-\mathcal{H}_{t}(x_{t})\right)\right\|^{2}+$$

$$\displaystyle\sum_{m=1}^{M}w_{t}[m]r\log\left(r+\left\|Q^{-1/2}\left(x_{t+1}-\mathcal{F}_{m}(x_{t})\right)\right\|^{2}\right)\bigg{)}.$$

$$\nu \Vert w_{t+1}-w_t{\Vert }_2^2,$$

$$\displaystyle\min_{x_{t},w_{t}}f(x,w):=\sum_{t=0}^{T-1}\bigg{(}\frac{1}{2}\left\|R^{-1/2}\left(y_{t}-\mathcal{H}_{t}(x_{t})\right)\right\|^{2}$$

$$+\sum _{m=1}^M{w}_t[m]r\log \left(r+{\Vert Q^{-1/2}\left(x_{t+1}-{\mathcal{F}}_m(x_t)\right)\Vert }^2\right)$$

$$\displaystyle+\nu\|w_{t+1}-w_{t}\|_{2}^{2}+I(w_{t}|\Delta^{M})\bigg{)},$$

$$I(w_t|{\Delta }^M):=\{\begin{array}{ll}0,& w_t\in {\Delta }^M;\\ \infty ,& w_t\notin {\Delta }^M.\end{array}$$

$$v(x):=\underset{w}{\min }f(x,w)$$

$$v_{\beta }(x):=\underset{w}{\min }f(x,w)+\frac{\beta }{2}\Vert w{\Vert }^2.$$

$$\nabla v_{\beta }(x)={\partial }_{x}f(x,w){|}_{w=w(x)}.$$

$$\underset{x}{\min }{v}_{\beta }(x)=$$

$$+$$

$$F(x)=\left(\begin{array}{c}{f}_1(x)\\ f_2(x)\end{array}\right)$$

$$f_1(x)=$$

$$+\nu \Vert w_{t+1}(x)-w_t(x){\Vert }_2^2+\frac{\beta }{2}\Vert w(x){\Vert }^2$$

$$f_2(x)=$$

$$\rho \left(\begin{array}{c}c\\ u\end{array}\right)=c+\frac{1}{2}\Vert u{\Vert }_{R^{-1}}^2+{\delta }_{[0,+\infty ]}(c).$$

$$d^{\ast }\in {\text{argmin}}_d$$

$$\in {\text{argmin}}_d$$

$$+\frac{1}{2}{d}^{T}U(x)d$$

$$\in {\text{argmin}}_d$$

$$+\nabla v_{\beta }(x{)}^{T}d$$

$$U=\left[\begin{array}{cccc}{U}_1& A_2^T& 0& \\ A_2& U_2& A_3^T& 0\\ 0& \ddots & \ddots & \ddots \\ & 0& A_T& U_T\end{array}\right]$$

$$A_t=-r\sum _{m=1}^M{w}_{t-1,m}(x)\frac{Q^{-1}\nabla {\mathcal{F}}_m(x_{t-1})}{r+\Vert x_t-{\mathcal{F}}_m(x_{t-1}){\Vert }^2},$$

$$U_t=$$

$$+$$

$$U_T=\frac{r{w}_{T-1,m}(x)Q^{-1}}{r+\Vert x_T-{\mathcal{F}}_m(x_{T-1}){\Vert }^2}.$$

Full text

††thanks: This material is based on work supported in part by the Washington Research Foundation Data Science Professorship, in part by the U. S. Army Research Laboratory and the U. S. Army Research Office under grant #W911NF-16-1-0158, and in part by the National Science Foundation Cyber-Physical Systems program award #1565529.

, , ,

1 Introduction

This paper considers the problem of using noisy measurements from a piecewise-continuous trajectory to estimate a hybrid system’s state. The state estimation problem has been extensively studied in classical dynamical systems whose states evolve according to one (possibly time–varying) smooth model. This problem is fundamentally more challenging for hybrid systems since the set of discrete state111We refer to the discrete component of the hybrid system state as the discrete state, and similarly refer to the continuous state, although the state of the hybrid system is specified by both the discrete and continuous components.

sequences generally grows combinatorially in time.

When the discrete state sequence and switching times are known a priori or directly measured, only the continuous state needs to be estimated, yielding a classical state estimation problem; this approach has been applied to piecewise-linear systems [30, Chap. 4.5] and to nonlinear mechanical systems undergoing impacts [25]. When the discrete state is not known or measured, estimating both the discrete and continuous states simultaneously improves estimation performance. One approach uses a bank of filters, each tuned to one discrete state, and selects the discrete states as the filter with the lowest residual [10, §4.1]. This filter bank method has been applied to hybrid systems with linear dynamics [8, §4.1] [20], nonlinear dynamics [11], and jumps in the continuous state when the discrete state changes [9]. Likewise, particle filter methods for hybrid systems [14, 17, 29] use a collection of filters, identified as particles, and are applicable to more general nonlinear process dynamics. Particle filters and filter banks are effective when the number of discrete states and dimension of continuous state spaces are small.

Another approach formulates a moving-horizon estimator over both the continuous and discrete states, resulting in a mixed-integer optimization problem [13]. The inherently discrete nature of the problem formulation enables estimation of the exact sample when the discrete state switches, at the expensive of combinatorial growth of the set of discrete decision variables as the horizon increases. Multiple methods have been developed to mitigate the challenge posed by this combinatorial complexity. One approach entails summarizing past measurements and state estimates with a penalty term in the the objective function [19]. Another approach, applicable to systems with bounded noise, entails restricting the set of possible discrete state sequences using a priori knowledge of the system [1, 2].

An alternative approach to circumventing the combinatorial challenge entailed by exactly estimating the discrete state sequence involves relaxing the discrete state estimate to take on continuous values as in [7, 22]. The latter reference uses a sparsity-promoting convex program whose objective incorporates a nonsmooth penalty across all possible discrete state sequences, and guarantees the estimate converges to the true continuous and discrete states. Both approaches are formulated for piecewise-linear systems whose continuous states do not jump when switching between subsystems; in the language of hybrid systems, the continuous states are reset using the identity function.

Our approach and contributions

We propose an offline algorithm for estimating the state of hybrid systems with nonlinear dynamics, non–identity resets, and noisy process and observation models. Our starting point is the optimization perspective on generalized and robust state estimation [3, 4]. To formulate state estimation as a continuous optimization problem, we relax the discrete state to take on continuous values as in prior work. Unlike prior work, we model process noise using the Student’s $t$ distribution, which allows large innovations and makes the method applicable to systems with non-identity resets.

In combination, these elements yield a nonsmooth nonconvex continuous optimization formulation (Sec. 2). We develop a Gauss-Newton type algorithm to solve this problem and prove the algorithm globally converges to stationary points (Sec. 3). The problem formulation and algorithm is evaluated on piecewise-linear and -nonlinear hybrid system models (Sec. 4).

2 Problem formulation

2.1 Process and observation models

We consider a class of discrete-time switched systems

[TABLE]

where $m\in M$ indexes the continuously-differentiable process model $\mathcal{F}_{m}\colon\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$, $\mathcal{H}_{t}\colon\mathbb{R}^{n}\rightarrow\mathbb{R}^{d}$ is the continuously-differentiable measurement model that generates observations $y_{t}\in\mathbb{R}^{d}$ of the hidden continuous state $x_{t}\in\mathbb{R}^{n}$, $\sigma_{t},\delta_{t}$ are process and measurement noises, and $w_{t}\in\mathcal{D}^{M}$ is a one-hot vector222 $w\in\mathbb{R}^{M}$ is one-hot if $w[i]\in\{0,1\}$ for all $i\in\left\{1,\dots,M\right\}$ and $1^{T}w=1$; $\mathcal{D}^{M}\subset\mathbb{R}^{M}$ denotes the set of one-hot vectors.that indicates which process model is active at time $t$. Note that the observations do not depend explicitly on the active model $m$, which must be inferred from measurements of the continuous state $x_{t}$.

The model $m$ that is active during each time step may be determined by an exogenous signal, prescribed as a function of time or state, or some combination thereof. Thus, the equation in (1) can represent the process and observation models of a wide variety of hybrid systems. We are motivated theoretically and experimentally to focus on cases where the active model is constant for many time steps, only occasionally switching to a new model. Such cases arise, for instance, when sampling trajectories of a hybrid dynamical system with a fixed timestep (so long as the timestep is much smaller then the average dwell time [21]).

When the discrete state changes in a hybrid system, the continuous state may change abruptly according to a reset map. As an example, the velocity of a rigid mass changes abruptly when it impacts a rigid surface [24]. Empirically, these discrete reset dynamics are much more poorly characterized than their continuous counterparts. For instance, whereas the ballistic trajectory of a rigid mass is well-approximated by Newton’s laws, the abrupt change in velocity that occurs at impact is not consistent with any established impact law [18]. Including such a reset in the system model (1) will introduce bias into the state estimate because the model will generate erroneous predictions at resets, diminishing the accuracy of estimated states at nearby times. This observation motivates us in the next section to account for the effect of unknown resets as part of the process noise.

2.2 Process noise and observation noise models

Instead of incorporating continuous state resets explicitly into the model (1), we introduce a distributional assumption on the process noise $\sigma_{t}$ that accepts large instantaneous changes in the continuous state estimate. Specifically, we assume that process noise $\sigma_{t}$ follows a Student’s $t$ distribution. Compared with the commonly-used Gaussian distribution, the heavy-tailed Student’s $t$ is tolerant to large deviations in the estimate of the hidden continuous state $x_{t}$ [6]. Hence, the Student’s $t$ error model allows an instantaneous change in the state that is consistent with (1) before and after the change. The negative log-likelihood of the Student’s $t$ (as a function of $\sigma_{t}$) is given by

[TABLE]

where $r$ is the degrees of freedom parameter of the Student’s $t$, and $Q$ is the covariance matrix.

If the continuous state $x_{t}$ was known, then any residual between the predicted observations $\mathcal{H}_{t}(x_{t})$ and actual measurements $y_{t}$ at time $t$ is due to measurement noise; in particular, the residual does not exhibit large deviations due to continuous state resets at switching times. Thus, we assume the measurement noise $\delta_{t}$ follows the usual Gaussian distribution, with negative log-likelihood

[TABLE]

where $R$ is the covariance matrix. The plots below provide a comparison between the probability density (left) and the negative log-likelihood (right) for the scalar Gaussian (solid blue) and Student’s $t$ distributions (dashed red; degree-of-freedom $r=1$).

probability density negative log-likelihood

2.3 State estimation problem formulation

The maximum a posteriori (MAP) likelihood for estimating states of (1) is given by

[TABLE]

Problem (2) is a nonlinear mixed-integer program with respect to both the continuous ($x_{t}$) and discrete ($w_{t}$) decision variables. We can significantly simplify the structure by establishing the following lemma.

Lemma 1** (Formulation Equivalence).**

Given $w\in\mathcal{D}^{M}$, any vectors $x_{1},x_{2}$, models $\mathcal{F}_{i}$, and any penalty functional $g$, we have

[TABLE]

Proof:

Since $w\in\mathcal{D}^{M}$ for both problems, there are only $m$ possible values for both objective functions, i.e.

[TABLE]

Hence, the minimum objective value for both problems will be $\min_{i}g(x_{2}-\mathcal{F}_{i}(x_{1}))$ and every minimizer is a one-hot vector that selects a minimum value. $\square$

Therefore, an equivalent formulation to (2) is given by

[TABLE]

Although still a mixed-integer program, this reformulation exhibits linear coupling between the discrete variables $w_{t}$ and continuous variables $x_{t}$. We will leverage this linear coupling when we develop our estimation algorithm based on the relaxed problem formulation introduced in the next section.

2.4 Relaxed state estimation problem formulation

Ultimately, the discrete state estimate will be specified as a one-hot vector, $w_{t}\in\mathcal{D}^{M}\subset\mathbb{R}^{M}$. To formulate a continuous optimization problem that approximates the mixed-integer problem formulated in the previous section, we relax the decision variable $w_{t}$ to take values in the convex hull $\Delta^{M}$ of $\mathcal{D}^{M}$.333 We use ${\Delta^{M}:=\{w\in[0,1]^{M}:1^{T}w=1\}}$ to denote the simplex in $\mathbb{R}^{M}$.

The optimal relaxed $w_{t}$ will generally lie on the interior of the simplex, so we project the result from our relaxed optimization problem to return the one-hot discrete state estimate. Since this relaxation-optimization-projection process tends to induce frequent changes in the discrete state estimate, we introduce a smoothing term on $w_{t}$,

[TABLE]

yielding the continuous relaxation of (3) given by

[TABLE]

where $x$ is the concatenated variable containing all $x_{t}$, $w$ is the concatenated variable containing all $w_{t}$, and convex indicator function $I$ is defined by

[TABLE]

The optimal relaxed discrete state estimate $w_{t}\in\Delta^{M}$ is projected onto $\mathcal{D}^{M}$ by choosing the (unique) one-hot vector whose $\operatornamewithlimits{argmax}_{i}{w}_{t}[i]$ component is equal to 1.

3 State estimation algorithm

We develop an algorithm to solve the relaxed state estimation problem formulated in (4) using two key ideas:

nonsmooth variable projection; 2. 2.

Gauss-Newton descent with Student’s $t$ penalties.

These two ideas are explained in the next two subsections, followed by a convergence analysis in the third subsection.

3.1 Nonsmooth variable projection

The first idea is to pass to the value function, projecting out (partially minimizing over) the $w$ variables. Define

[TABLE]

with $f(x,w)$ as in (4). The objective $f(x,w)$ is convex in $w$, but not strictly convex. To guarantee differentiability of $v(x)$, we add a smoothing term and consider

[TABLE]

where $\beta$ is usually taken to be a very small number (e.g. $10^{-4}$) so that the added term has minimal effect on the original value function. (The minimizer of $v_{\beta}$ is different from that of $v$.) The function $v_{\beta}(x)$ is a Moreau envelope [27, Def 1.22] of the true value function $v$; we refer the interested reader to [5] for details and examples concerning the Moreau envelope specifically (and nonsmooth variable projection more broadly). The unique minimizer $w(x)$ can be found quickly and accurately since the minimization problem with respect to $w$ is strongly convex: projected gradient descent converges linearly and can be accelerated using the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) [12] approach. With the minimizer $w(x)$, the gradient of $v_{\beta}$ is readily computed as

[TABLE]

Plugging $w(x)$ back into (4) we obtain the problem

[TABLE]

where $w_{t,m}(x)\equiv w_{t}[m](x)$.

3.2 Gauss-Newton descent with Student’s $t$ penalties

In this section we write the objective function (8) as a convex composite function and specify the line search we use to update $x$. The particular line search method was first proposed in [15] for a general class of algorithms for convex composite problems; the Gauss-Newton update derived here is a special case of this general line search scheme.

To cast the objective into a convex composite function, let $v_{\beta}=\rho\circ F$, where

[TABLE]

with

[TABLE]

and

[TABLE]

At each iteration, we choose a search direction $d^{*}(x)$ that

[TABLE]

where the equivalence is obtained by dropping terms independent of $d$. In general $U(x)$ can be any positive semidefinite matrix that varies continuously with respect to $x$, but for our particular objective function involving Student’s t penalty, $U(x)$ is chosen to be a Hessian approximation of the Student’s $t$ term in $f_{1}(x)$. Therefore the update can be interpreted as a Gauss-Newton style update. This approximation, proposed in [6, (5.5), (5.6)], is employed here because of its significant computational advantage; it is of the form

[TABLE]

with

[TABLE]

[TABLE]

for $1\leq t\leq T-1$, and

[TABLE]

We can rewrite $U(x)$ as

[TABLE]

where

[TABLE]

and

[TABLE]

Clearly $U(x)$ is positive semidefinite; we show in Lemma 3 that $U(x)$ is actually positive definite, so problem (9) reduces to the block triadiagonal linear system

[TABLE]

Given $d^{*}(x)$, the new $x^{+}$ is of the form

[TABLE]

where $\delta$ is a step size selected using the Armijo-type [26, Sec. 3.1] line search criterion.

[TABLE]

with

[TABLE]

When $d=0$, we have $\Delta(x;0)=0$444We overload $\Delta$ here to match the notation in [15, 6]; $\Delta(x;d^{*})$ should not be confused with $\Delta^{M}$, which is used to denote the simplex containing relaxed state estimates., and since we choose the minimizing

[TABLE]

we have $\Delta(x;d^{*})\leq 0$. Further,

[TABLE]

by [15, Thm. 3.6]. In other words, stationarity is achieved when $\Delta(x;d^{*})=0$. When $\Delta(x;d)<0$, we are guaranteed to have descent

[TABLE]

since $U(x)$ is positive semidefinite. This condition ensures that the line search step (11) is well-defined [15, Lemma 2.3].

Our approach is summarized in Algorithm 1. The positive parameter $\epsilon$ in the algorithm specifies the stopping condition. Finally, we project the relaxed discrete state estimate ${w}_{t}\in\Delta^{M}$ to obtain a discrete state estimate in $\mathcal{D}^{M}$ as described in Section 2.4.

3.3 Convergence of state estimation algorithm

The convergence of Algorithm 1 to a stationary point for a general class of convex composite objective functions is established in [15] and [6]. In particular  [6, Theorem 5.1] establishes the possible outcomes when applying this type of algorithm; informally, either the algorithm converges or the search direction $d_{k}$ diverges. In the remainder of this section we provide two technical results needed to formalize this intuition:

• •

Lemma 2 establishes a set of sufficient conditions that prevent divergence ($\|d^{(k)}\|\to\infty$);

• •

Lemma 3 proves that the sufficient conditions are satisfied.

Lemma 2**.**

Let $\Lambda=\{y|\rho(y)\leq v_{\beta}(x^{(0)})\}$. If $F^{-1}(\Lambda)=\{x|F(x)\in\Lambda\}$ is bounded and $U(x)$ is positive definite for all $x\in F^{-1}(\Lambda)$, then the hypotheses in [6, Theorem 5.1] are satisfied and the sequence of search directions $\{d^{(k)}\}$ is bounded.

Proof: The hypotheses in [6, Theorem 5.1] require that $F^{(1)}$ to be bounded and uniformly continuous on the set $S=\bar{co}(F^{(-1)}(\Lambda))$ where $\bar{co}$ stands for the closed convex hull. $F^{(1)}$ is continuous on $S$ since $f_{1}^{(1)}$ exists and is continuous by property of Moreau envelope and proximal operator, and $f_{2}^{(1)}$ is continuous trivially. Further, given that $S$ is closed by definition and bounded by assumption, it is compact. Hence $F^{(1)}$ is bounded and uniformly continuous on $S$.

Now we need to show that the sequence of search direction is bounded. At any iteration, the search direction $d$ we choose satisfies

[TABLE]

where the first inequality relies on $\rho\geq 0$ and on the positive semidefinite property of $U(x)$; the second inequality comes from $\Delta(x;d)\leq 0$; the third inequality results from the line search condition that creates a decreasing sequence $\{\rho(F(x^{(k)})\}$.

Since $\rho(F(x^{0}))$ is finite, $d^{T}U(x)d<\infty$ for all iterations. Because $\Lambda$ is closed by closedness of $\rho$ and $F$ is continuous, $F^{-1}(\Lambda)$ is also closed. Along with its boundedness by assumption, $F^{-1}(\Lambda)$ is compact. Since $x\in F^{-1}(\Lambda)\mapsto\lambda_{\min}(U(x))$ is continuous, its image is bounded, hence given that $U(x)$ is positive definite there exists some $\lambda_{\min}>0$ for all $x\in F^{-1}(\Lambda)$. Therefore $0<\lambda_{\min}\|d\|^{2}\leq d^{T}U(x)d<\infty$, which implies that $d^{(k)}$ cannot be unbounded. $\square$

Lemma 3**.**

$F^{-1}(\Lambda)$* is bounded for problem (8) and $U(x)$ is positive definite for all $x\in F^{-1}(\Lambda)$.*

Proof: First note that $\Lambda$ is bounded by the coercivity of $\rho$. This implies that for an unbounded sequence $\|x^{(k)}\|\to\infty$, we still have $f_{1}(x^{(k)})<\infty$ and $\|f_{2}(x^{(k)})\|<\infty$.

If $\|x^{(k)}\|\to\infty$, then we can find some $t+1$ and a subsequence $J$ such that $\lim_{k\in J}\|x_{t+1}^{(k)}\|=\infty$. By the definition of $f_{1}$ and $f_{1}(x^{(k)})<\infty$, $\lim_{k\in J}\|\mathcal{F}_{i}(x_{t}^{(k)})\|=\infty$, which further implies that $\lim_{k\in J}\|x_{t}^{(k)}\|=\infty$. Iteratively this means that $\lim_{k\in J}\|x_{t}^{(k)}\|=\infty$ for all $t$, in particular for the given starting point $x_{0}$, but that is not possible.

To show that $U(x)$ in (10) is positive definite, recall that we can rewrite $U(x)$ as

[TABLE]

If there exists some $d$ such that $d^{T}U(x)d=0$, then

[TABLE]

since $\tilde{Q}_{m}(w(x))^{-1}=\text{diag}(\tilde{Q}_{m,t}(w(x))^{-1})$, and

[TABLE]

are positive semidefinite. However because each $w_{t}\in\Delta$, there has to be some $\tilde{Q}_{m,t}^{-1}\succ 0$ for each $t$. Therefore $U(x)$ must be positive definite for all $x\in F^{-1}(\Lambda)$. $\square$

4 Experiments

To evaluate the proposed approach to state estimation for hybrid systems, we apply our algorithm to linear and nonlinear impact oscillators. In addition to being well-studied ([16, §1.2], [28]), these mechanical systems were chosen since they are among the simplest physically-relevant models that satisfy our assumptions, including non–identity reset maps. The parameter and trajectory regime considered in what follows is representative of a jumping robot constructed from one limb of a commercially-available quadrupedal robot [23] and controlled with an event-triggered stiffness adjustment; Figure 1(a) contains a photograph of the limb. The jumping robot’s hip and foot are constrained to move vertically in a gravitational field, so the rigid pantograph mechanism depicted in Figure 1(b) has two mechanical degrees-of-freedom (DOF) coupled through nonlinear pin-joint constraints. These two DOF are preserved, but their nonlinear coupling is neglected, in the piecewise-linear model illustrated in Figure 1(c). The hybrid dynamics of these linear and nonlinear impact oscillators are specified in Section 4.1

We perform two sets of experiments. The first set of experiments in Sec. 4.2 concern the piecewise-linear model depicted in Figure 1(c) and explore the consequences of our modelling assumptions and the efficacy of our proposed algorithm:

• •

Sec. 4.2.1 demonstrates the advantage of employing a Student’s $t$ distribution for process noise as compared to a Gaussian distribution;

• •

Sec. 4.2.2 demonstrates the superior convergence rate yielded by Gauss-Newton descent directions as compared to gradient (steepest) descent;

• •

Sec. 4.2.3 demonstrates the advantage of smoothing the relaxed discrete state estimate; and

• •

Sec. 4.2.4 demonstrates the algorithm’s performance when onboard measurements are used instead of offboard measurements.

The second set of experiments in Sec. 4.3 evaluate our proposed approach using the nonlinear model depicted in Figure 1(b).

Since this section is devoted to comparing estimated states to ground truth simulation results, and since our approach entails the determination of a relaxed discrete state estimate en route to obtaining the discrete state estimate, we now introduce notation that distinguishes these quantities:

• •

$w_{t}$ $\in\mathcal{D}^{M}$ denotes the ground truth discrete state;

• •

$\widetilde{w}_{t}$ $\in\Delta^{M}$ denotes the relaxed discrete state estimate;

• •

$\widehat{w}_{t}$ $\in\mathcal{D}^{M}$ denotes the discrete state estimate.

This notational distinction was not introduced previously in the interest of readability since there was no ambiguity entailed by overloading notation in the problem formulation and algorithm specification.

4.1 Impact oscillator hybrid system models

The continuous state $x=(q,\dot{q})\in\mathbb{R}^{4}$ for the jumping robot hybrid system model consists of the two-dimensional configuration vector $q\in\mathbb{R}^{2}$ and corresponding velocity $\dot{q}\in\mathbb{R}^{2}$, where $q[1]$ and $q[2]$ denote the vertical height of the hip and foot, respectively. The foot is not permitted to penetrate the ground, $q[2]\geq 0$, so the first part of the discrete state indicates whether this constraint is active: A (air) if $q[2]>0$, G (ground) if $q[2]=0$. To compensate for energy losses at impact, an event-triggered controller stiffens or softens a spring based on which direction the hip is traveling, so the second part of the discrete state indicates the direction of travel for $q[1]$: $\uparrow$ if up, $\downarrow$ if down. With $\ddot{q}_{m}(q,\dot{q})\in\mathbb{R}^{2}$ denoting the acceleration of the hip and foot in discrete state $m\in\left\{\text{A}\!\downarrow,\text{G}\!\downarrow,\text{G}\!\uparrow,\text{A}\!\uparrow\right\}$, formulae for this acceleration are given in Table 1. At the moment of impact (when the discrete state changes from $w_{t}\in\{\text{A}\!\downarrow,\text{A}\!\uparrow\}$ to $w_{t+1}\in\{\text{G}\!\downarrow,\text{G}\!\uparrow\}$) the foot velocity $\dot{q}[2]$ is instantaneously reset to [math], corresponding to perfectly plastic impact. An example of the jump in continuous state when transitioning from $\text{A}\!\downarrow$ to $\text{G}\!\downarrow$ on the foot velocity $\dot{q}[2]$ is shown in Figure 2 near time 17.5s.

4.2 Piecewise-linear impact oscillator experiment

In this subsection, we employ the linear spring laws

[TABLE]

[TABLE]

with parameter values $m_{h}=3,m_{t}=1,g=2$.

In our first demonstration the observed states are $q[1]$ and $q[2]$, position of the hip and foot, leaving the velocities unobserved:

[TABLE]

State estimation results for this system are shown in Figure 5.

In the remainder of this subsection, we demonstrate the effects of the choices we made in our problem formulation (Sec. 2) and algorithm derivation (Sec. 3) using the piecewise-linear model as a running example. We also consider a variation where the measurements correspond to the leg length and velocity, which are more representative of the onboard measurements available to an autonomous robot operating outside of the laboratory.

4.2.1 Student’s $t$ versus Gaussian process noise

Figure 2 compares the estimation of foot velocity using Student’s $t$ with $r=0.01$ versus using Gaussian for the process noise distribution; in both cases the true discrete state is given. The estimated trajectory for both distributions match the true simulated trajectory away from jumps, while near jumps, such as around times $16.6$s and $17.5$s, using the Student’s $t$ distribution enables closer tracking of the instantaneous change in the true foot velocity $\dot{q}[2]$ than when using a Gaussian distribution.

4.2.2 Gauss-Newton versus gradient (steepest) descent

We empirically compared convergence rates for continuous state $x_{t}$ updates obtained using Gauss-Newton and gradient (steepest) descent directions (Algorithm 1, line 3). Figure 3 shows the log loss versus algorithm iteration for the two methods; the actual discrete state $w_{t}$ was taken as given to perform this comparison. As expected, the objective value decreases significantly faster when the search direction is determined by the Gauss-Newton scheme as compared to the direction of steepest descent, reaching the stopping criterion in ten times fewer iterations in our tests.

4.2.3 Smoothing the relaxed discrete state versus not

If the continuous states are given, the discrete state estimate returned by our algorithm (skipping lines 2 and 3 of Algorithm 1) is very close to the true discrete state regardless of whether a smoothing term is included in the relaxed problem formulation. When simultaneously estimating both the continuous and discrete states, the smoothing term becomes crucial, as illustrated by comparing the discrete state estimates ($\widehat{w}_{t}$) in Figure 4 (without smoothing) and Figure 5 (with smoothing). In particular, the estimated discrete state switches rapidly without smoothing, whereas with smoothing the discrete state tends to remain constant for many samples and change mostly near ground-truth switching times.

4.2.4 Onboard versus offboard measurements

In the laboratory, the positions of the robot hip and foot can be directly measured offboard, e.g. with an external camera system. Outside of the laboratory, only the relative position of the hip and foot can be directly measured onboard our robot. Thus, we are motivated by this practical consideration to evaluate our algorithm’s performance in the case where only the relative position and velocity of the hip and foot are measured,

[TABLE]

Although the full hybrid system state is formally unobservable with these relative measurements, our algorithm nevertheless yields good estimates of the discrete state as shown in Figure 6; due to large errors in the estimate of (unobservable) continuous states, we omit those results from the figure.

4.3 Piecewise-nonlinear impact oscillator experiment

To test Algorithm 1 on a nonlinear model, we included the kinematic constraints depicted in Figure 1(b), resulting in a nonlinear spring force. In this model we set the two spring laws to be the same $k_{1}=k_{2}$, decreasing the number of discrete states from four to two: $w=\text{A}$ when $q[2]>0$ and $w=\text{G}$ when $q[2]=0$. State estimation results compare favorably with the analogous results from the piecewise-linear system when using either absolute position measurements $\mathcal{H}_{\text{pos}}$ (12) (compare Figure 7 with Figure 5) or relative measurements $\mathcal{H}_{\text{relative}}$ (13) (compare Figure 8 with Figure 6).

In Figure 7 we see that the model can estimate continuous and discrete states in the nonlinear setting. However, we do notice that the estimated trajectories are not as close to ground truth as in the linear case. In particular, when $q[2]$ has a value only slightly greater than 0 (e.g. between times 3s and 4s), the algorithm fails to detect the transition between $w=\text{A}$ and $w=\text{G}$.

5 Conclusion

We proposed a new state estimation algorithm for hybrid systems, analyzed its convergence properties, and evaluated its performance on piecewise-linear and -nonlinear hybrid systems with non-identity resets. The algorithm leverages a relaxed state estimation problem formulation where the decision variables corresponding to the discrete state are allowed to take on continuous values. This relaxation yields a continuous optimization problem that can be solved using recently-developed nonsmooth variable projection techniques. The effectiveness of the approach was demonstrated on hybrid system models of mechanical systems undergoing impact.

{ack}