1

arXiv:1604.00217v2 [cs.SY] 4 Apr 2018

Moving horizon estimation for discrete-time linear
systems with binary sensors: algorithms and
stability results
G. Battistellia , L. Chiscia , S. Gherardinia,b
a Università di Firenze, Dipartimento di Ingegneria dell’Informazione (DINFO), Via di Santa Marta 3,
50139 Firenze, Italy.
b CSDC, Università di Firenze, INFN, and LENS, Via G. Sansone 1, I-50019 Sesto Fiorentino, Italy,
and QSTAR, Largo E. Fermi 2, I-50125 Firenze, Italy.
{giorgio.battistelli, luigi.chisci, stefano.gherardini}@unifi.it

Abstract
The paper addresses state estimation for linear discrete-time systems with binary (threshold) measurements. A
Moving Horizon Estimation (MHE) approach is followed and different estimators, characterized by two different
choices of the cost function to be minimized and/or by the possible inclusion of constraints, are proposed. Specifically,
the cost function is either quadratic, when only the information pertaining to the threshold-crossing instants is
exploited, or piece-wise quadratic, when all the available binary measurements are taken into account. Stability
results are provided for the proposed MHE algorithms in the presence of unknown but bounded disturbances and
measurement noise. Performance of the proposed techniques is also assessed by means of simulation examples.

Keywords: State estimation; moving-horizon estimation; binary measurements; stability analysis.
I. I NTRODUCTION
Binary (threshold) sensors whose output can take two possible values according to whether the sensed variable
exceed or not a given threshold, are nowadays commonly exploited for monitoring/control aims in a wide range of
application domains. A non-exhaustive list of existing binary sensors includes: industrial sensors for brushless
dc motors, liquid levels, pressure switches; chemical process sensors for vacuum, pressure, gas concentration
and power levels; switching sensors for exhaust gas oxygen (EGO or lambda sensors), ABS, shift-by-wire in
automotive applications; gas content sensors (CO, CO2 , H2 , etc.) for gas & oil industry; traffic condition indicators
for asynchronous transmission mode (ATM) networks; medical sensors/analyses with dichotomous outcomes. In
some applications, binary sensors represent the only viable solution for real-time monitoring. In any case, they
provide a remarkably more cost-effective alternative to traditional (continuous-valued) sensors at the price of an
accuracy deterioration which can, however, be compensated by using many binary sensors (for different variables
and/or thresholds) in place of a single one or few traditional sensors. Moreover, binary (threshold) measurements
arise naturally in the context of networked state estimation when, in order to save bandwidth and reduce the energy
consumption due to data transmission, the measurements collected by each remote sensor are compared locally with
a (possibly time-varying) threshold and only information pertaining to the threshold-crossing instants is transmitted
to the fusion center. This latter setting falls within the framework of event-based or event-triggered state estimation
[1], [2], [3], and is more challenging as compared to the usually addressed settings due to the minimal information
exchange.
The above arguments, as well as the difficulties due to the very limited information provided by binary measurements, has motivated the work on the exploitation of binary measurements for estimation purposes. In particular,
[4], [5] investigated observability and observer design for linear time-invariant (LTI) continuous-time systems under
binary-valued output observations. The work in [6], [7] addressed system identification using binary sensors. Specific
attention was also devoted to state estimation of hybrid nonlinear systems with binary/quantized sensors [8] and to
target tracking with binary sensor networks [9]. A possible solution for coping with the high nonlinearity associated
with binary measurements within a stochastic framework is particle filtering [10], [11]. However such techniques,
while effective in many contexts, suffer from the so-called curse of dimensionality (i.e., the exponential growth

2

of the computational complexity as the state dimension increases) and from the lack of guaranteed stability and
performance (being based on Monte Carlo integration).
The present paper addresses state estimation for linear discrete-time systems with binary (threshold) output measurements by following a moving horizon estimation (MHE) approach. MHE techniques were originally introduced
to deal with uncertainties in the system knowledge [12] and, in recent years, have gathered an increasing interest
thanks to their capability of taking explicitly into account constraints on state and disturbances in the filter design
[13], and on the possibility of having guaranteed stability and performance even in the nonlinear case [14], [15],
[16]. In fact, MHE has been successfully applied in many different contexts, ranging from switching and large-scale
systems [17], [18], [19], [20], [21] to networked systems [22], [23], [24].
In this paper, the state estimation problem with binary measurements is cast in a deterministic framework, in the
sense that no probabilistic description of the plant disturbance and noises is supposed to be available. The estimates
are computed by minimizing suitable cost functions defined over a given time-horizon (advancing in time) of finite
length, possibly subject to linear inequality constraints accounting for the threshold measurements. Specifically,
two different approaches are proposed and analyzed. In the first approach, only the threshold-crossing instants are
taken into account in the definition of the cost function, by penalizing the distance of the expected continuous
outputs (based on the state estimates) from the threshold at those instants. The main advantage of this solution is
that the resulting cost function is quadratic. The second approach, instead, exploits all the available information by
defining a piece-wise quadratic cost function which accounts for all the available binary measurements, but requires
the solution of a convex optimization problem at each time instant. Both unconstrained and constrained MH state
estimators will be presented for the two different choices of the cost function and stability results will be proved,
assuming unknown but bounded disturbances.
Summarizing, the paper provides the following contributions.
• Design of novel receding-horizon state estimators for linear discrete-time systems subject to binary (threshold)
measurements using either a quadratic or a piecewise quadratic cost function to be minimized and, independently, either including or not constraints.
• Stability analysis showing that all proposed estimators, irrespectively of the cost being used and of the inclusion
of constraints, guarantee an asymptotically bounded estimation error under bounded disturbances and suitable
observability assumptions.
• Performance comparison demonstrating the effectiveness, in terms of both estimation accuracy and computational cost, of our approach.
Some of the results of this paper have been preliminarily presented, without proof, in [25].
The rest of the paper is structured as follows. Section 2 formulates the estimation problem of interest. Section
3 discusses how to solve the problem by means of the MHE approach, with different variants depending on the
choice of the cost function as well as on the inclusion or not of constraints. Section 4 deals with the stability
analysis of the proposed MH estimators. In section 5, some numerical examples are presented in order to evaluate
and compare the proposed estimators. Finally, section 6 ends the paper with concluding remarks and perspectives
for future work.
II. P ROBLEM FORMULATION AND PRELIMINARY CONSIDERATIONS
The following notation will be used throughout the paper: col(·) denotes the matrix obtained by stacking its
arguments one on top of the other; diag(m1 , . . . , mq ) denotes the diagonal matrix whose diagonal elements are the
scalars m1 , . . . , mq ; further, given a matrix M, vec(M) denotes the linear transformation which converts the matrix
4

M into a column vector and kvkM = v0 Mv. Finally, ⊗ denotes the Kronecker product.
Let us consider the problem of recursively estimating the state of the discrete-time linear dynamical system
xt+1 = Axt + But + wt
zti = Ci xt + vti , i = 1, . . . , p

(1)

from binary (threshold) measurements
yti

=

hi (zti )


=

+1, if zti ≥ τ i
−1, if zti < τ i

(2)

3

p
In (1)-(2): xt ∈ Rn is the state to be estimated;
ut ∈ Rm is a known input; zt = col zti i=1 ∈ R p ; τ i is the threshold
p

p
of the i−th binary sensor; A, B,C = col Ci i=1 are matrices of compatible dimensions; wt and vt = col vti i=1 are
the process and, respectively, measurement noises assumed unknown but bounded. Notice from (1)-(2) that sensor
i provides a binary measurement yti ∈ {−1, +1} (two-level measurement quantization) according to whether the
noisy linear function of the state zti = Ci xt + vti falls below or above the threshold τ i . The problem (1)-(2) clearly
includes, as a special instance, the case of quantized sensors with an arbitrary number of levels. In fact, a d-level,
for generic d ≥ 2, quantizer can be easily realized by using d − 1 binary (threshold) sensors for the same physical
variable but with appropriate different thresholds. The considered setting with multiple binary sensors (which can
measure the same physical variable with different thresholds but also different physical variables) is clearly more
general.
It is worth to point out that the system (1)-(2) represents a very special instance of a linear system with output
nonlinearity, i.e. a Wiener system [26]. However, due to the discontinuous nature of the measurement function (2), all
those state estimation techniques for Wiener systems that require a certain smoothness of the output nonlinearity (see
for example [27] and the references therein) cannot be applied. In fact, while general-purpose nonlinear estimators
accounting for such a discontinuity (e.g., the particle filter) could be used, the peculiar nature of the considered
output nonlinearity deserves special attention and, for optimal exploitation of the poor available information, the
development of ad-hoc receding-horizon estimators that will be presented in the sequel.
Before addressing the estimation problem, some preliminary considerations on the information provided by binary
multisensor observations are useful. With this respect, it has been pointed out in [5] that, in the continuous-time case,
the information provided by a binary sensor of the form (2) is strictly related to the threshold-crossing instants. In
fact, in this case, at every instant corresponding to a discontinuity of the binary signal yi , it is known that the signal
zi is equal to the threshold value τ i , implying that the linear measurement zi = τ i is available. Hence, observability
with binary sensors for continuous-time linear systems can be analyzed within the more general framework of
observability for irregularly sampled systems [5]. In particular, observability can be ensured when the number of
threshold-crossing instants (which corresponds to the number of available irregularly sampled linear measurements)
is sufficiently large.
The situation is, however, different for discrete-time systems. To see this, consider a generic time instant k in
which the binary signal yik changes sign, i.e., yik yik+1 < 0. Then, it is not possible to state, as in the continuous-time
case, that zik coincides with the threshold τ i . Conversely, it can be simply concluded that there exists α ∈ [0, 1] such
that
α zik + (1 − α) zik+1 = τ i ,
(3)
the exact value of α being clearly unknown and unobservable from the binary measurements. Notice that (3)
simply states that if the binary output yik switches from discrete time k to k + 1, then the threshold τ i must lie in the
interval between zik and zik+1 . In view of (3), such discrete time instants k at which the output of some binary sensor
changes value will be more appropriately referred to as output switching or simply switching instants, instead of
threshold-crossing instants like in the continuous-time case considered in [5]. It is easy to see that (3) corresponds
to an uncertain linear measurement
α zik + (1 − α) zik+1 = Ci xk + δki + ηki ,

(4)

where δki is the uncertainty and ηki the measurement noise given by
δki = (1 − α)Ci (A − I)xk + (1 − α)Ci Buk ,
ηki = α vik + (1 − α) vik+1 + (1 − α)Ci wk .
As a consequence, even in presence of bounded disturbances, the uncertainty associated with the measurement (3)
depends on xk and uk . Recalling that, in general in the context of state estimation for uncertain systems, boundedness
of the state trajectories is a prerequisite for the boundedness of the estimation error - see, for instance, the discussion
in Section 2.1 of [28] - our attention will be restricted to the case of bounded state and input trajectories by making
the following assumption.
A1 At any time t, the vectors xt , ut , wt , vti , i = 1, . . . , p, belong to the compact sets X, U, W , and V i , i = 1, . . . , p,
respectively.

4

In practice, the compact sets X, U, W , V i need not be known by the estimator; they will only be used for stability
analysis purposes.
Remark 1: When the discrete-time system is obtained by sampling a continuous-time one with system matrices
(Ac , Bc ,Ci ), then the amplitude of the uncertainty δki can be related to the sampling interval Ts . In fact, it turns out
R
that, since in this case A = eAc Ts and B = 0Ts eAct B dt, δki vanishes as Ts goes to zero and, in addition, when Ts is
small δki ≈ Ts (1 − α)Ci Ac xk + (1 − α)Ci Bc uk .
III. M OVING HORIZON ESTIMATION FOR BINARY SENSORS
In order to estimate the state xt of the linear system (1) given the binary measurements (1)-(2), a MHE approach
is adopted. Then, by considering a sliding window Wt = {t − N,t − N + 1, . . . ,t}, the goal is to find estimates of
the state vectors xt−N , . . . , xt on the basis of the information available in Wt and of the state prediction xt−N at the
beginning of Wt . Let us denote by x̂t−N|t , . . . , x̂t|t the estimates of xt−N , . . . , xt , respectively, to be obtained at any
stage t.
Following the discussion at the end of the previous section, a first natural approach for constructing a MH estimator
would amount to considering the information provided by the switching instants inside the sliding window Wt , in
order to define the cost-function to be minimized. Accordingly, for any time instant t ≥ N and for any sensor index
i, let us define the set Iti of switching instants as
Iti = {k ∈ Wt : k + 1 ∈ Wt and yik yik+1 < 0}.

(5)

Then, the following least-squares cost function can be defined
t−1

JtA = kx̂t−N|t − xt−N k2P +

∑
k=t−N

p

kx̂k+1|t − Ax̂k|t − Buk k2Q + ∑

∑i kCi x̂k|t − τ i k2Ri ,

(6)

i=1 k∈It

where the positive definite matrices P ∈ Rn×n , Q ∈ Rn×n and the positive scalars Ri , i = 1, . . . p, are design parameters
to be suitably chosen. The first term, weighted by the matrix P, penalizes the distance of the state estimate at the
beginning of the sliding window from the prediction xt−N . The second contribution, weighted by the matrix Q, takes
into account the evolution of the state in terms of the state equation (1). Finally, for each sensor i the third term
weighted by the scalar Ri penalizes the distances of the expected output (based on the state estimates) Ci x̂k|t from
the threshold τ i at the switching instants. Notice that considering the distance from the threshold at the switching
instant is equivalent, for sampled-data systems, to considering the beginning of the time interval [kTs , (k + 1)Ts ]
in which the threshold crossing happens. As a matter of fact, since for a sampled-data system a binary sensor
does not provide a precise information on the threshold crossing instant in the interval [kTs , (k + 1)Ts ], considering
the distance from the threshold at the beginning of the time interval is just a choice, not necessarily optimal. As
an alternative, with little modifications, one could consider for instance the middle point of the interval. Such
modifications would not affect the properties (e.g. stability) of the estimator.
Thus, at each time t ≥ N, the estimates in the window Wt can be obtained by solving the following optimization
problem.
Problem EtA : Given the prediction xt−N , the input sequence {ut−N , . . . , ut−1 }, and the sets Iti , i = 1, . . . , p, find
◦
◦ that minimize the cost function (6).
the optimal estimates x̂t−N|t
, . . . , x̂t|t
A , different prediction
Concerning the propagation of the estimation procedure from Problem EtA to Problem Et+1
strategies may be adopted. For instance, a first possibility consists of assigning to xt−N+1 the value of the estimate
◦
of xt−N+1 made at time instant t, i.e., x̄t−N+1 = x̂t−N+1|t
. As an alternative, following [15], the state equation of the
◦
noise-free system can be applied to the estimate x̂t−N|t
. In this case, the predictions are recursively obtained by
◦
xt−N+1 = Ax̂t−N|t
+ But−N , t = N, N + 1, . . . .

(7)

Such a recursion is initialized with some a priori prediction x0 of the initial state vector. Hereby, this latter possibility
will be adopted as it will facilitate the derivation of the stability results (see Section 4).
The main positive feature of Problem EtA is that it admits a closed-form solution since the cost function (6)
depends quadratically on the estimates x̂t−N|t , . . . , x̂t|t (for the readers’ convenience an explicit expression for the

5

solution is reported in the Appendix). On the other hand, such a cost takes into account only the information
pertaining to the switching instants, which, however, is intrinsically uncertain as discussed in the previous section.
In order to overcome such a limitation, a different cost function can be considered by taking into account all the
time instants in the sliding window Wt . To this end, for any sensor i = 1, . . . , p, let us define the functions


1, if zi − τ i yi < 0
i i i
ω (z , y ) =
(8)
0, otherwise
Suppose now that at time k the sensor i provides a measurement yik = 1. Then, the information provided by such a
measurement is that the linear measurement zik is above the threshold τ i , i.e., belongs to the semi-interval [τ i , +∞).
Such information can be included in the cost function by means of a term of the form ω(Ci x̂k|t , 1) kCi x̂k|t − τ i k2Ri
which penalizes the distance of the expected output Ci x̂k|t from [τ i , +∞). Similarly, in the case yik = −1, a term
of the form ω(Ci x̂k|t , −1) kCi x̂k|t − τ i k2Ri can be used to penalize the distance of the expected output Ci x̂k|t from
(−∞, τ i ]. Summing up, the inclusion of such terms gives rise to a cost function of the following form
p

t−1

JtB = kx̂t−N|t − xt−N k2P +

∑

kx̂k+1|t − Ax̂k|t − Buk k2Q + ∑

k=t−N

t

∑

ω i (Ci x̂k|t , yik )kCi x̂k|t − τ i k2Ri .

(9)

i=1 k=t−N

While a closed-form expression for the global minimum of (9) does not exist, since JtB is piece-wise
quadratic, it is

B
i
i
i
easy to see that the cost Jt enjoys

 some nice properties. In fact, while each function ω C x̂k|t , yk per se is discontinuous, the product ω i Ci x̂k|t , yik kCi x̂k|t − τ i k2Ri is continuous since at the points of discontinuity of ω i Ci x̂k|t , yik , i.e.,
for Ci x̂k|t = τ i , the product vanishes. Further, for similar reasons, also the derivative 2ω i Ci x̂k|t , yik Ri (Ci )0 (Ci x̂k|t −τ i )
of the product turns out to be continuous even at Ci x̂k|t = τ i . Thus the product ω i Ci x̂k|t , yik kCi x̂k|t − τ i k2Ri is
continuously differentiable on Rn . Hence, the overall cost function JtB is continuously differentiable with respect to
the estimates x̂t−N|t , . . . , x̂t|t and also strictly convex (since P > 0 and Q > 0). Hence, standard optimization routines
can be used in order to find its global minimum. Clearly, since an optimization has to be performed, it is also
reasonable to include constraints accounting for the available information on the state trajectory so that the solver
can work on a bounded solution set. In particular, in order to preserve convexity, it is advisable to consider a
convex set X containing X (if X is convex, one can simply set X = X; in general, choosing X as a convex
polyhedron is preferable so that only linear constraints come into play). Then, at any stage t = N, N + 1, . . ., the
following optimization problem has to be solved.
Problem EtB : Given the prediction xt−N ,
i
{yt−N , . . . , yti , i = 1, . . . , p}, find the optimal
◦ ∈ X for k = t − N, . . . ,t.
constraints x̂k|t

the input sequence {ut−N , . . . , ut−1 }, the measurement sequences
◦
◦ that minimize the cost function (9) under the
estimates x̂t−N|t
, . . . , x̂t|t

Also in this case, the predictions xt−N are supposed to be recursively obtained via equation (7) starting from
a prior prediction x0 . Of course, if no information on the set X is available or if it is preferable to resort to an
unconstrained optimization routine, one can simply let X = Rn .
As a final remark, it is worth pointing out that for the two previously presented optimization problems there is
a trade-off between estimation accuracy and computational cost. In fact, the cost in Problem EtA is quadratic but
accounts only for part of the information provided by the sensors, while Problem EtB accounts for all the available
information but requires a convex optimization program to be solved.
Some considerations on the computational complexity of the proposed approaches are in order. The solution of
Problem EtA requires simply the minimization of a strictly convex quadratic form in (n+1)N variables,where n is the
plant order. Standard techniques like Gaussian elimination can solve this kind of problems with complexity O(n3 N 3 )
but faster algorithms are available. This means that this approach is much computationally cheaper as compared to
particle filtering algorithm which usually require in the order of O(10n ) particles to provide satisfactory performance.
As for the solution of Problem EtB , it entails the minimization of a convex and continuously differentiable piecewise
quadratic cost function. It is known that this kind of problems can be solved in finite time by means of sequential
quadratic programming [29]. Further, many computationally efficient algorithms are available which are able to
handle problems with hundreds of optimization variables [30], [31] and enjoys super-linear convergence [32].
Nevertheless, application of Problem EtB is possible only when the number n of state variables is not too large and
the sampling interval is sufficiently long so as to allow the optimization to terminate. In the other cases, one must
resort to Problem EtA .

6

A. Accounting for additional constraints
Provided that some information on the bounds of the process disturbance wt and measurement noises vti is available, additional constraints can be considered in the determination of the state estimates. For instance, considering
a convex (usually polyhedral) set W containing W , one can impose the constraints
x̂k+1|t − Ax̂k|t − Buk ∈ W ,

k = t − N, . . . ,t − 1

(10)

in the solution of the optimization problem. Moreover, assuming the knowledge of upper bounds ρVi on the
amplitudes |vti |, i = 1, . . . , p, of the measurement noises, for each k and each i, the constraints
(
Ci x̂k|t < τ i + ρVi , if yik = −1
(11)
Ci x̂k|t > τ i − ρVi , if yik = 1
can be imposed. With this respect, it is an easy matter to see that the constraints in (11) define a polyhedron in the
state space as summarized in the following proposition (the proof is reported in the Appendix).

Proposition 1: Given the vector χ̂t = vec [x̂t−N|t · · · x̂t|t ]0 of the estimates in the observation window, the constraints in (11), for k = 0, . . . , N and i = . . . , p, can be written in compact form as
Γt χ̂t < γt ,

(12)

where
Γt = [Φt (C ⊗ IN )] ∈ R pN×nN ,


γt = Φt vec(T 0 ) + vec(V ) ∈ R pN ,
p
1
2
Φt = −diag(yt−N
, . . . , yt1 , yt−N
, . . . , yt2 , . . . , yt−N
, . . . , ytp ) ∈ R pN×pN ,


 1

ρV1 · · · ρV1
τ · · · τ1

..
..  ∈ R p×N , V =  ..
..
..  ∈ R p×N .
T =  ...
 .
.
.
.
. 
p
p
p
τ ··· τ
ρV · · · ρVp

(13)

While the inclusion of the constraints (10) and (12) in the convex optimization problem EtB is natural, in some
circumstances it may be interesting to combine them also with the quadratic cost JtA . For example, minimizing
JtA under the linear constraints (12) can be a way to account for the information concerning the non switching
instants without the necessity of considering the piece-wise quadratic cost. In fact, this would result in a quadratic
programming problem (being the cost quadratic and the constraints linear) for which many efficient solvers are
available. It is worth to point out that what is, among the above mentioned options, the best choice clearly depends
on the situation under consideration and, in particular, on the available computational resources, on the available
information (the bounds ρVi may be unknown), and on the necessity (or not) of having estimates satisfying the
constraints (since clearly this property is guaranteed only if the constraints are taken into account in the estimator
design). Nevertheless, in the next Section it will be shown that both costs JtA and JtB imply some nice stability
properties of the resulting MH estimator.
Remark 2: Notice that while the considered system dynamics is linear, we do not have access to the linear
measurements zt = Cxt + vt but rather to the nonlinear (binary) measurements yti = hi (zti ) , for which we cannot
apply neither the Kalman filter due to nonlinearity of hi (·) nor the extended Kalman filter due to the discontinuous
nature of hi (·). It is however worth noting that the simplified quadratic cost JtA amounts to considering a fictitious
linear measurement of the form Ci xk = τ i + ηki for each switching instant k in the observation window. In this case
and supposing that no constraints are imposed, the estimates could be computed also via a Kalman-like filter. In
all the other cases, i.e. when the piecewise quadratic cost JtB is used or constraints are imposed in the optimization,
this is no longer possible.

7

IV. S TABILITY ANALYSIS
The focus of this section is on the analysis of the stability properties of the state estimators obtained by solving, at
each time instant, either Problem EtA or EtB . Specifically, a complete analysis is first provided in the more involved
case of Problem EtB . This will be followed by a short discussion on the main differences in the analysis with respect
to Problem EtA . Notice that the analysis carried out in [15] for the nonlinear case cannot be directly applied in
the present context, since the binary sensors do not satisfy the observability requirement of [15]. The proofs of all
results reported in this section can be found in the Appendix.
For each sensor i and for each time instant t ≥ N, let us denote by Θti the observability matrix concerning the
set Iti of the switching instants in the observation window Wt , i.e,
Θti = col(Ci Ak−t+N )k∈Iti .

(14)

Then, the observability matrix related to the switchings in Wt of all binary sensors is
p
Θt = col(Θti )i=1
.

(15)

Please notice that the observability matrix defined in (14)-(15) is actually related to the linear subsystem (1), with
output zt , of the overall system (1)-(2) considering only those discrete-time instants at which some binary sensor
output switches.
The following uniform observability assumption is needed in order to ensure that enough information is provided
by the binary sensors in each window Wt .
A2 For any t ≥ N, rank(Θt ) = n, with n = dim(xt ).
Remark 3: The above uniform observability assumption is made in accordance with the observation that each
output switching can be associated with a linear (albeit uncertain) measurement of the form (3). Hence, each
switching instant k can be thought of as a sampling instant for the linear output zik . This means that observability of
the system depends crucially on the output switching instants in each observation window which, in turn, clearly
depend on the thresholds and of the time window length N. In practice, the threshold (or the thresholds when
multiple sensors are available) and the time window length N must be chosen taking into account the system
dynamics so as to ensure that such an irregular sampling preserves observability. For instance, when only one
binary sensor is available, clearly N should be substantially greater than 2n − 1, with n = dim(xt ), so as to ensure
that at least n output switching instants are present in each observation window. While some analytical results on
observability under irregular sampling are available [5], the simplest approach amounts to studying, for instance
by numerical simulations, how the observability measure δ varies as a function of the thresholds and of the time
window length N. See for instance Figure 2 in Section 5 concerning the considered case study. Of course, depending
on the system dynamics, time-invariant thresholds may not be sufficient to always ensure uniform observability
(think for example to the case of a constant linear output). In these cases, observability can be recovered by making
each threshold oscillate in the range of variability of the corresponding continuous output zti with a sufficiently high
frequency and by choosing N so that each observation window contains a sufficient number of threshold oscillation
periods. This latter solution is particularly convenient in case zti is a measurement collected by a remote sensor and
a time-varying threshold τti is used for transmission scheduling.
To see this, notice first that sufficient conditions relating the rank of an observability matrix under irregular sampling
to the number of samples and to the eigenvalues of the state transition matrix can be found along the lines of [5].
Specifically, for an observable sampled-data system (1)-(2), the observability matrix Θt defined in (14)-(15) has
full rank n if the number of switchings νt in the window Wt is such that
νt /N ≥ 2(n − 1)/N + ωmax /π
4

(16)

where ωmax = maxλ ∈sp(A) ∠λ , sp(A) being the spectrum (set of eigenvalues) of A and ∠λ the argument of λ ∈ C.
Notice that ωmax can be interpreted as the bandwidth of the system (in radians). Hence it turns out that, asymptotically
for large N, the condition (16) amounts to requiring that the density of switchings νt /N be greater than or equal to
ωmax /π < 1, which represents the system-bandwidth to Nyquist-bandwidth ratio. This means that, for large values
of N, uniform observability can be ensured even with a single binary sensor by making the threshold oscillations
sufficiently fast so as to ensure that the density of switching exceeds ωmax /π < 1.

8

Before stating the main stability results, some preliminary definitions are needed. Given a symmetric matrix S,
let us denote by λ (S) and λ (S) the minimum and maximum eigenvalues of S, respectively. Further, given a matrix
M, let us denote by kMk , λ (M 0 M)1/2 its norm. Given a generic subset Ψ of an Euclidean space, let us define
ρΨ , supv∈Ψ kvk. Given a generic quantity Gi related to the i−th binary sensor, let us define G , maxi kGi k and
G , mini kGi k. Finally, let us define the uniform observability measure associated to the matrices Θt as
δ = inf kΘt k = inf λ (Θt0 Θt )1/2 .
t≥N

t≥N

Notice that, under assumption A2, it can be stated that δ > 0. The following result can now be stated.
◦
Theorem 1: Let assumptions A1 and A2 hold. For each t ≥ N, let the estimate x̂t−N,t
be generated by solving
◦
B
.
Problem Et , with xt−N recursively obtained via equation (7), and consider the estimation error et−N , xt−N − x̂t−N|t
Then, the weighted norm of the estimation error can be recursively bounded as
ket−N k2P ≤ a1 ket−N−1 k2P + a2 , t = N, N + 1, . . .

(17)

b1 kAk2
,
b2
2 + c kBk2 ρ 2 + c ρ 2 + c ρ 2
c1 kA − Ik2 ρX
2
3 W
4 V
U
a2 =
,
b2
!



d1
1
δ 2R
λ (P)
4+
d2 + R , b2 =
+
b1 =
λ (P)
λ (Q)
2 4λ (P)

(18)

where
a1 =

and c1 , c2 , c3 , c4 , d1 , d2 are suitable constants (given in the proof). In addition, if the weights Q and Ri , i = 1, . . . , p,
are selected such that a1 < 1, the norm of the estimation error turns out to be asymptotically bounded in that

1/2
a2
◦
lim sup ket−N k ≤ e∞ ,
.
1 − a1
t→+∞

The reason for analyzing the estimate at the beginning of the observation window is that, due to the nature of
◦
the MHE estimation scheme, the estimate x̂t−N|t
is used to generate the prediction x̄t−N+1 used at time t + 1. This
◦
◦
makes it possible to recursively write et−N+1 = xt−N+1 − x̂t−N+1|t+1
as a function of et−N = xt−N − x̂t−N|t
. Notice
that even in the noise-free case, i.e., when the process disturbance and the measurement noise are zero and hence
ρW = ρV = 0, the asymptotic bound e◦∞ on the estimation error does not go to zero due to the presence of the term
2 + c kBk2 ρ 2 in a . Indeed, such a term accounts for the intrinsic uncertainty associated with the
c1 kA − Ik2 ρX
2
2
U
threshold-crossing instants in discrete-time as discussed at the end of Section 2 (see equation (3) and the subsequent
discussion). With this respect, it is worth recalling that, when the discrete-time system under consideration is obtained
by sampling a continuous-time system, the quantities kA − Ik and kBk vanish as the sampling interval Ts goes to
zero. This means that the smaller is the sampling interval, the smaller turns out to be the asymptotic bound on the
estimation error since the information concerning the threshold-crossing instants becomes more precise.
Another important issue concerns the solvability of the stability condition a1 < 1. In particular, the following
result can be readily proved.
Proposition 2: Let assumption A2 hold. Then, when δ > 0, it is always possible to select the weights P, Q and
Ri , i = 1, . . . , p, so that a1 < 1. In particular, for given Q and Ri , i = 1, . . . , p, the condition a1 < 1 can be satisfied
by letting P = εP, with P any positive definite matrix, provided that ε is suitably small.
Hence, if the observability measure δ is strictly positive, it is sufficient to choose P sufficiently small in order
to ensure the satisfaction of the stability condition a1 < 1. This result is in accordance with the well-known results
on stability of MHE algorithms which stipulate that stability is ensured provided that the weight on the prediction
is sufficiently small [15].
◦
Remark 4: In the statement of Theorem 1 the estimates x̂t−N,t
are generated by solving Problem EtB , in which the
◦ ∈ X for k = t − N, . . . ,t are present. From the practical point of view, including such constraints in
constraints x̂k|t

9

the optimization is useful in order to take into account bounds on the state variables in the design of the estimator. In
fact, in many contexts, for instance when the bounds correspond to some physical constraints, providing estimates
outside bounds can be meaningless. On the other hand, in other cases, the bounds on the state variables can
be unknown. The proposed approach is flexible enough to deal also with such a case since the inclusion of the
constraints is not necessary for the stability of the estimation error dynamics. In fact, in the unconstrained case,
◦ , x̂◦
the lower bound of each term ι(x̂k|t
k+1|t ) (see the Appendix) can be derived as follows


◦ 2
◦
◦
◦
ι(x̂k|t
, x̂k+1|t
) ≥ kCi x̂k|t
− τ i k2Ri − 3(Li )2 kA − Ik2 kx̂k|t
k + kBk2 ρU2 + kw◦k|t k2
◦ k2 = kl˜ k2 , where l˜ , col(x̂◦ )
◦
◦
and ∑k∈Iti kx̂k|t
t|t
t|t
k|t k∈Iti = Φt x̂t−N|t + Γt ũt + Λt w̃t . It can be readily observed that the
matrices Φt , Γt and Λt are proportional to Θti , Hti and Dti . More precisely, Θti = (Ci ⊗ In )Φt , Hti = (Ci ⊗ In )Γt and
Dti = (Ci ⊗ In )Λt .
◦
Consider now the case in which, for each t ≥ N, the estimate x̂t−N|t
is generated by solving Problem EtA , with xt−N
◦ , . . . , x̂◦
recursively obtained via equation (7). Notice that, in this case, no constraint is imposed on the estimates x̂t−N
t
which can be readily obtained as the unique global minimum of the strictly convex quadratic function JtA . A close
inspection of the proof of Theorem 1 shows that the same line of reasoning can be applied also for Problem EtA .
◦ , x̂◦
The main difference is that, when deriving the lower bound for the optimal cost, each term ι(x̂k|t
k+1|t ) can be
i
◦
i
2
A
simply replaced with the quantity kC x̂k|t − τ k in accordance with the definition of cost Jt . Then an inequality
analogous to (17) can be derived, with the important difference that, in the definition of the novel a2 , ρX can be
replaced by ρX (which is consistent with the fact that the constraint set X is not used in the solution of Problem
EtA ).
Remark 5: While the foregoing analysis does not account for the possible presence of the additional constraints
discussed in Section III-A, analogous results could be easily obtained also when the constraints (10) and/or (12)
are imposed in the determination of the state estimates. In this case, the bound on the estimation error turns out to
be smaller thanks to the additional information provided by such constraints.
Remark 6: As a final remark, it is pointed out that the extension of the stability results reported here to the case
in which the binary measurements are obtained by thresholding nonlinear output maps and/or the system dynamics
is nonlinear does not entail particular conceptual difficulties, by combining the analysis of Theorem 1 with that of
[15], [16]. On the other hand, in this case, establishing a link between the observability properties and the number
of threshold crossing instants (see Remark 2) appears more challenging. Further, for nonlinear output maps, the
resulting cost functions need not be convex.

V. N UMERICAL RESULTS
In this section, numerical examples are presented in order to show the effectiveness of the proposed MHE
algorithms for binary measurements. In particular, two different case-studies will be considered: a first simple
example concerning a 2-mass 2-spring oscillator and a single binary sensor just for the sake of testing the algorithms’
capabilities on a critically observable system, and a second example on a network of 2-mass 2-spring oscillators
with multiple binary sensors to illustrate a more realistic application of the estimators.
Example 1
Let us consider the 2-mass 2-spring mechanical system of Fig. 1. The state of the system is defined as x =
[x1 , ẋ1 , x2 , ẋ2 ]0 where x1 and x2 are the displacements of the two masses from their static equilibrium positions.
Accordingly the system is described by the continuous-time linear state equations ẋ(t) = Ac x(t) with


0
1
0
0
k2
− (k1 +k2 ) 0
0


m1
m1
Ac = 
(19)

0
0
0
1


k2
0 − mk22 0
m2
where k1 , k2 are the stiffnesses of the springs and m1 , m2 the corresponding masses. The parameters are set to
m1 = 1 = m2 = 1 [Kg], k1 = k2 = 10 [N/m], and the continuous-time model is discretized with sampling interval
Ts = 0.1 [s]. Further, it is assumed that only the displacement x2 (third state component) is measured by a single

10

Fig. 1. 2-mass 2-spring mechanical oscillator of example 1.

threshold sensor so that the output matrix turns out to be C = [0, 0, 1, 0]. In all the simulations, the initial state is
chosen so as to impose the harmonic motion condition, i.e. x0 = [0.618, 0, 1, 0]0 , making the two masses oscillate
with the same frequency but different amplitudes within the interval [−1, 1]; the initial phase of the oscillations is
a uniformly distributed random variable. The process disturbance is taken equal to zero, while the measurement
noise is a white sequence with uniform distribution in the interval [−ρV , ρV ]. In order to tune the proposed MHE
algorithms for appropriate performance, the threshold value τ of the binary sensor and the length N of the estimation
sliding window need to be properly selected. To this end, it has been analyzed by means of numerical simulations
how the observability measure δ varies as a function of N and τ, as shown in Fig. 2 with a simulation time interval
of 50 [s] and a noise level ρV = 0.05. As shown in Fig. 2, observability requires sufficiently large window size (
N ≥ 60 with τ = 0.5). Also notice that the observability measure as a function of N has a monotonically increasing
behaviour with some characteristic plateaus. Further, it is perfectly symmetric with respect to τ: if the threshold
value is outside the range [−1, 1] of the system output, then no information is provided by the binary sensor; τ = 0
also implies poor observability as sampling the sinusoid in proximity of zero provides little information about the
sinusoid amplitude. From Fig. 2, we chose N = 100 and τ = 0.5 for the forthcoming simulation results, so that
assumption A2 holds. For the weight matrices we selected Q = I4 , R = 1 and P = εI4 with ε < 10−4 in order to
satisfy the stability condition a1 < 1 according to Proposition 2. Hereafter, for the sake of brevity, the filter obtained
by solving at each time instant Problem EtA will be referred to as least-squares MHE (LSMHE) algorithm. Similarly,
piece-wise MHE (PWMHE) will indicate the filter obtained by solving Problem EtB . Figs. 3 and 4 show the time
behavior of the true state variables and of the corresponding estimates in a random simulation with a single binary
sensor by using, respectively, the LSMHE and PWMHE algorithm with measurement noise level ρV = 0.05 and the
same initialization for both algorithms. As it can be seen, although the estimators are initialized far from the true
initial state (ρX = 5) and the amount of information exploited in the cost function is limited, the estimates resulting
from both algorithms converge to the true trajectories of the systems state vector, and, as expected, the PWMHE
algorithm exhibits much better performance in the transient thanks to the additional information taken into account
in the definition of cost function JtB . In order to better appreciate the accuracy of the proposed algorithms and take
into account the timescales of the systems, Monte Carlo simulations have been performed by randomly varying the
measurement noise realization, the phase of the oscillations for the true state trajectories, and the a priori prediction
x0 , which is randomly generated with uniform distribution in [−5, 5]4 . For the sake of comparison, also a particle
filter with standard sequential importance sampling and 103 particles was tested in the same setting. The results are
not reported here because the implemented particle filter was not able to converge and to track the true state. As
performance index, in Fig. 5 we have adopted a relative error, i.e. the root mean square error (RMSE) normalized
by the Euclidean norm of the true system state, where
!1
L ke k2 2
t,l
RMSE(t) = ∑
,
(20)
L
l=1

11

Fig. 2. Example 1 - (a) Observability measure δ as a function of the length N of the estimation sliding window (with τ = 0.5). (b)
Observability measure δ as a function of the threshold value τ (with N = 100). The results in (a)-(b) have been evaluated over 100 Monte
Carlo trials.
N
1
5
20
35
50
100
150
CPU TIME ( IN [ S ])

LSMHE
PWMHE
0.50·10−3
0.25
0.56·10−3
0.42
1.79·10−3
1.11
3.23·10−3
2.07
5.30·10−3
3.19
22.83·10−3
7.53
78.90·10−3
15.70
TABLE I

PER ITERATION STEP FOR DIFFERENT VALUES OF

N.

and et,l is the state estimation error at time t in the l−th simulation run and L = 100 is the number of Monte Carlo
trials. Fig. 5 confirms the effectiveness of the MHE algorithms for state estimation with binary observations.
The computational burden of solving both Problems EtA and EtB , as a function of the length N of the estimation
sliding window, has been evaluated by means of the CPU time per iteration step (a notebook with an Intel Core
i7-2640M CPU @ 2.80 GHz has been used in simulations). The results are reported in Table I. Notice that PWMHE
is by far more computationally expensive than LSMHE (computing time three orders of magnitude larger in this
specific small-size example). As a matter of fact, the solution of Problem EtA can be found analytically by an
explicit matrix formula, while for the solution of EtB a convex mathematical programming problem has to be

12

Fig. 3. Ground truth (dashed blue line) and estimates (solid red line) of the state components versus time, by solving the estimation problem
EtA .

Fig. 4. Ground truth (dashed blue line) and estimates (solid red line) of the state components versus time, by solving the estimation problem
EtB .

solved. However, it is worth to point out that the PWMHE algorithm has been implemented by using standard
functions of the Matlab Optimization Toolbox, without resorting to ad-hoc optimization routines. Hence, we are
confident that much faster computing times can be achieved. The dependence of performance on the threshold τ
and the noise level ρV is analyzed in Fig. 6, where the ARMSE (i.e., the asymptotic RMSE defined as the average
of the RMSE after the transient computed in the time interval [25, 40] [s]) is reported for the LSMHE algorithm.
Also to compute the ARMSEs in Fig. 6, the RMSEs have been normalized by the Euclidean norm of the true
system state. As observed from Fig. 6a, when the threshold of the sensor is close to zero, performance undergoes
a substantial deterioration. Such a behavior is due to the fact that as τ goes to zero, the observability measure δ

13

Fig. 5. Example 1 - Normalized RMSEs of the LSMHE and PWMHE filters, evaluated over 100 Monte Carlo trials.

Fig. 6. (a) ARMSE as a function of the threshold τ (with ρV = 0.05 and N = 100). (b) ARMSE as a function of the measurement noise
level ρV (with τ = 0.5 and N = 100). All the results have been evaluated over 100 Monte Carlo trials.

associated with the output switching instants becomes small. Moreover, as shown in Fig. 6b, the ARMSE decreases
almost linearly as the noise level decreases, but, even when the noise is zero, the ARMSE does not go to zero due
to the intrinsic uncertainty associated with binary measurements as discussed at the beginning of Section 2. For
similar reasons, even if the process disturbance is taken equal to zero in the simulations, the second term of the
cost functions (6) and (9) does not go to zero, since the estimates need not coincide with the true state and hence
x̂k+1|t is in general different from Ax̂k|t .
Example 2
Finally, in order to numerically assess the performance of the proposed MHE algorithms when the dimensionality
of the system state and the number of binary sensors increase, the network in Fig. 7 of six coupled 2-mass 2-

14

Fig. 7. Network of six coupled 2-mass 2-spring oscillators. Each node of the network has a binary sensor, monitoring the corresponding
third state component.

Fig. 8. Normalized RMSEs of the LSMHE and PWMHE filters, evaluated over 100 Monte Carlo trials, for a network of six 2-mass 2-spring
oscillators.

spring oscillators (like the one in Fig. 1) is considered. It is assumed that each node is equipped with a binary
sensor measuring the third component of the local state vector, with threshold belonging to the range [−1, 1].
The network dynamics turns out to be described by a discrete-time linear dynamical system with matrices A =
I6 ⊗ Ad − γL ⊗ I4 and C = I6 ⊗ [0, 0, 1, 0], where: Ad = exp(Ac Ts ); L is the Laplacian matrix of the network; the
sampling interval is Ts = 0.1 [s]. For the sake of simplicity, we have chosen the same value γ = 0.02 for the
coupling constants between all the connected sites, which ensures the synchronization of the system states. Note
that synchronization is reached if γ < 0.31685. The threshold values of the six binary sensors are taken, respectively,
equal to [0.5, 0.2, −0.5, −0.8, −0.2, 0.3]0 . In all simulations, the initial state of each 2-mass 2-spring system is a
uniformly distributed random variable centred around the vector x0 = [0.618, 0, 1, 0]0 with variations of ±5 for each
component, while the measurement noise is a white sequence uniformly distributed in the interval [−0.05, 0.05].
Moreover, the validity of Proposition 2 for the network is ensured by choosing ε = 10−5 with P = I24 . The duration
of each simulation experiment is fixed to 35 [s], and the corresponding RMSE of the proposed MHE filters is
averaged over 100 Monte Carlo trials. In Fig. 8 the RMSEs, normalized by the Euclidean norm of the true system
state, of the LSMHE and PWMHE algorithms are plotted. It can be seen that, also in this case, the PWMHE filter
exhibits better performance in the transient, and that the convergence of its estimation error is slower by a factor
of approximately 4 with respect to the case of the single oscillator.

15

VI. C ONCLUSIONS
The paper considers state estimation for linear discrete-time systems, in which the available information is
provided by binary multi-sensor observations, in the presence of unknown but bounded noises affecting both the
system and the measurement. Two novel moving-horizon estimators have been introduced, resulting from the
minimization of a least-square and a piece-wise quadratic cost function, respectively, with the possible inclusion of
constraints. The stability of the estimation error dynamics for the proposed filters has been analyzed and related to
the measure of observability associated with the time instants in which the binary outputs switch. Two simulation
examples concerning respectively a single mechanical oscillator and a network of coupled oscillators, have been
worked out in order to demonstrate the effectiveness of the proposed approach.
ACKNOWLEDGEMENT
We acknowledge fruitful discussions with Daniele Mari and Nicola Forti.
A PPENDIX
EtA :

Closed-form solution of Problem
Let us consider the cost function (6). Under the assumption A1, (6) can be written as the following quadratic form:
0
0
0
0
0
Ut−N + rt−N ,
Mt−N Ŷt−N|t + 2Ŷt−N|t
Dt−N − Dt−N
Ŷt−N|t + rt−N = Ŷt−N|t
Mt−N Ŷt−N|t − Ŷt−N|t
JtA = Ŷt−N|t

where Ŷt−N|t = col(x̂t−N+i|t )Ni=0 ∈ RnN , Dt−N = −Ut−N ∈ RnN and the
defined as

P + A0 QA + ζ j,1
−A0 Q
0
0 QA + ζ
0Q

−QA
Q
+
A
−A
j,2

Mt−N = 
..
..
..

.
.
.
0

0

and


Ut−N

with
i
δ j,h

(21)

matrices Ut−N ∈ RnN , Mt−N ∈ RnN×nN are

···
0

···
0


..

···
.

0

Q + A0 QA + ζ j,N

···

A0 QBu0 − Pxt−N − π j,1
A0 QBu1 − QBu0 − π j,2
..
.









=
,
 0

A QBuN−1 − QBuN−2 − π j,N−1 
A0 QBuN − QBuN−1 − π j,N

(
1, if ∃ j ∈ Iti : j = h
=
0, else

, h = 1, . . . , N,

and
p

p

0

i
π j,h = ∑ δ j,h
Ci Ri τ i , h = 1, . . . , N,
i=1

0
rt−N = x̄t−N
Px̄t−N +

t

∑

0

i
ζ j,k = ∑ δ j,h
Ci RiCi ,
i=1

p

2

u0k B0 QBuk + ∑ hi Ri τ i ∈ R,

k=t−N

i=1

hi = dim(Iti ).
Necessary condition for the minimum of the cost function (21) is
∇Ŷt−N|t JtA (Ŷt−N|t ) = 2Mt−N Ŷt−N|t + 2Ut−N = 0,

(22)

◦
for any t = N, N + 1, . . .. Solving (22) as a function of x̂t−N|t , we obtain the optimal estimates x̂t−N|t
, t = N, N + 1, . . .
that minimize the cost function (6), namely


In 0
.{z
. . 0}
◦
−1
|
x̂t−N|t =
Mt−N
Dt−N , t = N, N + 1, . . .
(23)
∈R(N−1)n×n

16

Choosing the weighting matrices P and Q as positive semi-definite matrices and Ri > 0, the solution (23) corresponds
to a global minimum, since the Hessian matrix Mt−N of the cost function is strictly positive definite. As a final
remark, notice that there are many equivalent ways of writing the solution of Problem EtA and the particular
form presented here is a consequence of the fact that we consider as optimization variables the state estimates
x̂t−N+i|t for i = 0, . . . , N. An alternative would be to consider as optimization variables the state estimate x̂t−N|t
at the beginning of the observation interval together with the estimates of the process disturbance ŵt−N+i|t =
x̂t−N+i+1|t − Ax̂t−N+i|t − But−N+i for i = 0, . . . , N − 1. In this case, each x̂t−N+i|t would be written as a function of
x̂t−N|t and the observability matrix would explicitly appear in the solution.

Proof of Proposition 1:
For each k = t − N, . . . ,t − 1, we initially introduce the constraints for the i−th measurement equation, i = 1, . . . , p:
(
Ci x̂k|t < τ i + ρVi , if yik = −1
(24)
Ci x̂k|t > τ i − ρVi , if yik = 1
The system (24) is equivalent to the inequality
yik (Ci x̂k|t + yik ρVi ) > yik τ i .

(25)

Observing that (yik )2 = 1, ∀k = t − N, . . . ,t − 1, we obtain
yik (Ci x̂k|t − τ i ) + ρVi > 0, k = t − N, . . . ,t − 1.
p
p
If we define φk = diag(yik ) ∈ R p×p , i = 1, . . . , p, τ p = col(τ i )i=1
∈ R p and ν = col(ρVi )i=1
∈ R p , then we can write

φk (Cx̂k|t − τ p ) + ν > 0,
since φk0 φk = I p . Moreover, introducing the matrices Φt , T and V as in (13), the constraints (24) can be written in
matrix form, namely


Φt vec (CX̂t − T )0 < vec(V )
(26)

0



where X̂t = x̂t−N|t , . . . , x̂t|t . Observing that vec (CX̂t )0 ≡ (C ⊗ In )vec X̂t0 , it can be noted that (26) is equal to
(12), so that the proposition is proved.

Proof of Theorem 1:
Some preliminary definitions are needed. Notice first that, while the function ω(Ci x, y)kCi x−τ i k is not differentiable
for Ci x = τ i , for Ci x 6= τ i one has


∂
0,
if y(Ci x − τ i ) > 0 ,
ω(Ci x, y)kCi x − τ i k =
.
−yCi , if y(Ci x − τ i ) < 0 .
∂x
Hence ω(Ci x, y)kCi x − τ i k is globally Lipschitz with Lipschitz constant Li = kCi k, for i = 1, . . . , p. Further, consider
i = col(Ci x̂◦ )
for each sensor i and each sliding window Wt , the vector z̃t|t
k|t k∈Iti . Then, we can write
i
◦
z̃t|t
= Θti x̂t−N|t
+ Hti ũt + Dti w̃t◦ ,

where
ũt = col(uk )k∈[t−N,t] ,
◦
◦
w◦k|t = x̂k+1|t
− Ax̂k|t
− Buk ,

w̃t◦ = col(w◦k|t )k∈[t−N,t] ,
0

and Hti and Dti are suitable matrices. Let φ i be defined as supt≥N λ (Dti Dti )1/2 . Clearly, φ i is finite since Dti can
assume only a finite number of configurations in the estimation window.
◦ ; the aim is to find a lower and an upper bound
Let us now consider the estimation error as et−N = xt−N − x̂t−N
for the optimal cost
◦
Jt◦ = kx̂t−N|t
− xt−N k2P +

t−1

∑
k=t−N

p
◦
◦
kx̂k+1|t
− Ax̂k|t
− Buk k2Q + ∑

t

∑

i=1 k=t−N

◦
ω(zik , yik )kCi x̂k|t
− τ i k2Ri .

(27)

17

to derive a bounding sequence on the norm of the estimation error.
– Upper bound on the optimal cost Jt◦ :
For the optimality of the cost function Jt◦ , we have Jt◦ ≤ JtB
Jt◦ 6 kxt−N − xt−N k2P +

x̂k|t =xk , k∈Wt
p

t−1

kwk k2Q + ∑

∑
k=t−N

and hence

t

∑

ω(zik , yik )kzik − τ i k2Ri .

(28)

i=1 k=t−N

The discontinuous function ω(zik , yik ) is non zero if and only if zik − τ i ∈ V i , i.e. if the system output is close to the
i−th sensor threshold and the measurement noise makes the sensor detection incoherent with the system evolution.
Thus, the upper bound (28) can be rewritten as
Jt◦ ≤ kxt−N − xt−N k2P + Nλ (Q)ρW2 + p(N + 1)Rρ V2 .

(29)

– Lower bound on the optimal cost Jt◦ :
Let us consider a time instant k ∈ Iti and suppose, for the sake of notational simplicity, that yik = 1 and yik+1 = −1
(up-down threshold crossing). Note that the dual case can be analysed in a similar way. Thus, in the cost function
Jt◦ the following contribution is present:
◦
◦
◦
◦
− τ i k2Ri
− τ i k2Ri + ω(zik+1 , −1)kCi x̂k+1|t
) , ω(zik , 1)kCi x̂k|t
, x̂k+1|t
ι(x̂k|t


◦
◦
◦
− τ i k2Ri + ω(zik+1 , −1)kCi x̂k+1|t
− τ i k2Ri − ω(zik , −1)kCi x̂k|t
− τ i k2Ri ,
= ω(zik , 1) + ω(zik , −1) kCi x̂k|t
|
{z
}
=1, by definition

where
◦
◦
◦
ω(zik+1 , −1)kCi x̂k+1|t
− τ i k2Ri − ω(zik , −1)kCi x̂k|t
− τ i k2Ri ≤ (Li )2 k(A − I)x̂k|t
+ Buk + w◦k|t k2 .
◦ ∈ X for k = t − N, . . . ,t, it can be stated that each term ι(x̂◦ , x̂◦
Since x̂k|t
k|t k+1|t ) has a lower bound, such that


2
◦
◦
◦
+ kBk2 ρU2 + kw◦k|t k2 .
ι(x̂k|t
, x̂k+1|t
) ≥ kCi x̂k|t
− τ i k2Ri − 3(Li )2 kA − Ik2 ρX

Since k is switching instant, hence
yik = hi (Ci xk + vik ) = 1
and
yik+1 = hi (Ci Axk +Ci Buk +Ci wk + vik+1 ) = −1,
i.e. there exists α ∈ [0, 1] such that αzik + (1 − α)zik+1 = τ i , from which
τ i = Ci xk + ζki ,
where ζki = δki + ηki . Then,
1
◦
◦
◦
−Ci xk k2Ri − kζki k2Ri ,
kCi x̂k|t
− τ i k2Ri = kCi x̂k|t
−Ci xk − ζki k2Ri ≥ kCi x̂k|t
2
where

2
kζki k2Ri ≤ 4Ri kCi k2 kA − Ik2 ρX
+ kCi k2 kBk2 ρU2 + kCi k2 ρW2 + (ρVi )2 .
Summarizing the previous results, if we consider ∀i only the instants k ∈ Iti , we obtain
p


◦
◦
Jt◦ ≥ kx̂t−N|t
− xt−N k2P + ∑ ∑ kCi x̂k|t
−Ci xk k2Ri − βt − σt ,
i=1 k∈Iti

where
p

βt = ∑

∑i

i=1 k∈It

 i


2
2
4R kC j k2 kA − Ik2 ρX
+ kCi k2 kBk2 ρU2 + kCi k2 ρW2 + (ρVi )2 + 3(Li )2 (kA − Ik2 ρX
+ kBk2 ρU2 )

18

p

and σt = ∑

∑i 3(Li )2 kw◦k|t k2 are quantities with an upper bound. Indeed, it can be stated that:

i=1 k∈It




2
2
2
2
2
2
+C kBk2 ρU2 +C ρW2 + ρ V2 + 3p(N + 1)L kA − Ik2 ρX
+ kBk2 ρU2 = β̆t
βt ≤ 4p(N + 1)R C kA − Ik2 ρX
and


2
i
σt ≤ 3p max L
i

2

t−1

kw◦k|t k2 ≤

∑
k=t−N

i
3pL h
kx̂t−N|t − xt−N k2P + Nλ (Q)ρW2 + p(N + 1)Rρ V2 = σ̆t .
λ (Q)

To conclude the calculation of the lower bound, let us define z̃ti , col(zk )k∈Iti and R̃i , Ri I|Iti | and write
p

ψt , ∑

∑

i=1 k∈Iti

p
p


◦
i
◦
kCi x̂k|t
−Ci xk k2Ri = ∑ kz̃t|t
+ Hti ũt + Dti w̃t◦ − Θti xt−N − Hti ũt − Dti w̃t − ṽti k2R̃i ,
− z̃ti k2R̃i == ∑ kΘti x̂t−N|t
i=1

i=1

with w̃t , col(wk )k∈[t−N,t] and ṽti , col(vik )k∈Iti . Hence,

p 
1
1 i ◦
2
i ◦ 2
i
2
i 2
◦
kΘt (x̂t−N|t − xt−N )kR̃i − kDt w̃t kR̃i − kDt w̃t kR̃i − kṽt kR̃i ≥ kΘt (x̂t−N|t
− xt−N )k2R̃ − µ̆t ,
ψt ≥ ∑
4
4
i=1
where
p

µt

=

∑ Ri

 i 2


kDt k kw̃t◦ k2 + ρW2 + (ρVi )2

i=1

#
2

φ 
2
kxt−N − xt−N k2P + Nλ (Q)ρW2 + p(N + 1)Rρ V2 + φ ρW2 + ρ V2 = µ̆t ,
≤ pR
λ (Q)
"

i.e.
ψt ≥

δ 2R
4λ (P)

◦
kx̂t−N|t
− xt−N k2P − µ̆t =

In conclusion
◦
Jt◦ ≥ kx̂t−N|t
− xt−N k2P +

δ 2R
4λ (P)

δ 2R
4λ (P)

ket−N k2P − µ̆t .

ket−N k2P − β̆t − σ̆t − µ̆t .

(30)

Now we can exploit the bounds on the optimal cost Jt◦ in order to obtain a bounding sequence on the norm of the
estimation error. More specifically, combining (29) and (30), we derive the following inequality:
◦
kx̂t−N|t
− xt−N k2P +

But, noting that

δ 2R
4λ (P)

ket−N k2P ≤ β̆t + σ̆t + µ̆t + kxt−N − xt−N k2P + Nλ (Q)ρW2 + p(N + 1)Rρ V2 .

1
◦
kx̂t−N|t
− xt−N k2P ≥ ket−N k2P − kxt−N − xt−N k2P
2

and
xt−N − xt−N = Aet−N−1 + wt−N−1 ,
namely



kxt−N − xt−N k2P ≤ 2 kAk2P ket−N−1 k2P + λ (P)ρW2 ,

inequality (31) can be rewritten as
ket−N k2P ≤ a1 ket−N−1 k2P + a2 ,

(31)

19

where the coefficients a1 and a2 are defined as in (18) and
2

2

d2 = 3L φ̄ −2 ,


2
2
c1 = c2 = p(N + 1) 4R C + 3L ,
!
h
i
b1
2
2
c3 = b1 + Nλ (Q)
− 1 + pR 4(N + 1)C + φ ,
2λ (P)
!
b1
c4 = p(N + 1)R
− 1 + pR(4N + 5).
2λ (P)
d1 = 2pφ ,

Since a2 is a positive scalar, if we further impose that a1 < 1, the asymptotic upper bound e◦∞ can be easily derived,
in that
t−1

ket k2P < at1 ke0 k2P + a2 ∑ a1j ,
j=0

which tends to a2 /(1 − a1 ) as t → ∞.



Proof of Proposition 2: Notice first that the stability condition a1 < 1 can be rewritten as
!


2R

λ (P)
d1
δ
1
d2 + R kAk2 ≤
4+
+
.
λ (P)
λ (Q)
2 4λ (P)
By letting P = εP, with P any positive definite matrix, the above inequality becomes
!



λ (P)
d1
1
δ 2R
2
.
4+
d2 + R kAk ≤
+
λ (Q)
2 4 ε λ (P)
λ (P)
It can be seen that the left-hand side of such an inequality does not depend on ε, whereas the right-hand side goes
to infinity as ε goes to 0, provided that δ 2 > 0. Hence, when δ 2 > 0, it is always possible to ensure that the stability
condition a1 < 1 holds by taking any Q, Ri , i = 1, . . . , p, P, and then choosing ε suitably small.

R EFERENCES
[1] Battistelli, G. and Benavoli, A. and Chisci, L. Data-driven communication for state estimation with sensor networks. Automatica, vol.
48, pp. 926–935, 2012.
[2] Sijs, J. and Lazar, M. Event-based state estimation with time synchronous updates. IEEE Trans. on Automatic Control, vol. 57, pp.
2650–2655, 2012.
[3] Shi, D. and Chen, T. and Shi, L. Event-triggered maximum likelihood state estimation. Automatica, vol. 50, pp. 247–254, 2014.
[4] Wang, L.Y. and Xu, G. and Yin, G.G. State reconstruction for linear time-invariant systems with binary-valued output observations.
Systems and Control Letters, vol. 57, pp. 958–963, 2008.
[5] Wang, L.Y. and Li, G.G. and Guo, L. and Xu, C.-Z. State observability and observers of linear-time-invariant systems under irregular
sampling and sensor limitations. IEEE Trans. on Automatic Control, vol. 56, pp. 2639–2654, 2011.
[6] Wang, L.Y. and Zhang, J.F. and Yin, G.G. System identification using binary sensors. IEEE Trans. on Automatic Control, vol. 48, pp.
1892–1907, 2003.
[7] Wang, L.Y. and Yin, G.G. and Zhang, J.F. Joint identification of plant rational models and noise distribution functions using binaryvalued observations. Automatica, vol. 42, pp. 543–547, 2006.
[8] Koutsoukos, X.D. Estimation of hybrid systems using discrete sensors. Proceedings 42nd IEEE Conf. on Decision and Control, pp.
155–160, 2003.
[9] Aslam, J. and Butler, Z. and Constantin, F. and Crespi, V. and Cybenko, G. and Rus, D. Tracking a moving object with a binary sensor
network. Proceedings 1st ACM Conf. on Embedded Networked Sensor Systems, pp. 150–161, 2003.
[10] Djuric, P.M. and Vemula, M. and Bugallo, M.F. Target tracking by particle filtering in binary sensor networks. IEEE Trans. on Signal
Processing, vol. 56, pp. 2229–2238, 2008.
[11] Ristic, B. and Gunatilaka, A. and Gailis, R. Achievable accuracy in Gaussian plume parameter estimation using a network of binary
sensors. Information Fusion, vol. 25, pp. 42–48, 2015.
[12] Jazwinski, A.H. Limited memory optimal filtering. IEEE Trans. on Automatic Control, vol. 13, pp. 558–563, 1968.
[13] Rao, C.V. and Rawlings, J.B. and Lee, J.H. Constrained linear estimation–a moving horizon approach. Automatica, vol. 37, no. 10,
pp. 1619–1628, 2001.
[14] Rao, C.V. and Rawlings, J.B. and Mayne, D.Q. Constrained state estimation for nonlinear discrete-time systems: stability and moving
horizon approximations. IEEE Trans. on Automatic Control, vol. 48, no. 2, pp. 246–257, 2003.

20

[15] Alessandri, A. and Baglietto, M. and Battistelli, G. Moving horizon state estimation for nonlinear discrete-time systems: new stability
results and approximation schemes. Automatica, vol. 44, pp. 1753–1765, 2008.
[16] Alessandri, A. and Baglietto, M. and Battistelli, G. and Zavala, V.M. Advances in moving horizon estimation for nonlinear systems.
Proceedings 49th IEEE Conference on Decision and Control, pp. 5681–5688, 2010.
[17] Alessandri, A. and Baglietto, M. and Battistelli, G. Receding-horizon estimation for switching discrete-time linear systems. IEEE
Trans. on Automatic Control, vol. 50, no. 11, pp. 1736–1748, 2005.
[18] Guo, Y. and Huang, B. Moving horizon estimation for switching nonlinear systems. Automatica, vol. 49, no. 11, pp. 3270–3281, 2013.
[19] Farina, M. and Ferrari-Trecate, G. and Scattolini, R. Moving-horizon partition-based state estimation of large-scale systems. Automatica,
vol. 46, no. 5, pp. 910–918, 2010.
[20] Haber, A. and Verhaegen, M. Moving Horizon Estimation for Large-Scale Interconnected Systems. IEEE Trans. on Automatic Control,
vol. 58, no. 11, pp. 2834–2847, 2013.
[21] Schneider, R. and Hannemann-TamÃ¡s, R. and Marquardt, W. An iterative partition-based moving horizon estimator with coupled
inequality constraints. Automatica, vol. 61, pp. 302–307, 2015.
[22] Farina, M. and Ferrari-Trecate, G. and Scattolini, R. Distributed moving horizon estimation for linear constrained systems. IEEE Trans.
on Automatic Control, vol. 55, no. 11, pp. 2462–2475, 2010.
[23] Farina, M. and Ferrari-Trecate, G. and Scattolini, R. Distributed moving horizon estimation for nonlinear constrained systems.
International Journal of Robust and Nonlinear Control, vol. 22, no. 2, pp. 123–143, 2012.
[24] Liu, A. and Yu, L. and Zhang, W.-A. and Chen, M.Z.Q. Moving horizon estimation for networked systems with quantized measurements
and packet dropouts. IEEE Trans. on Circuits and Systems I: Regular Papers, vol. 60, pp. 1823–1834, 2013.
[25] Battistelli, G. and Chisci, L. and Gherardini, S. Moving horizon state estimation for discrete-time linear systems with binary sensors.
Proceedings 54th IEEE Conf. Decision and Control, pp. 2414–2419, 2015.
[26] Westwick, D. and Verhaegen, M. Identifying MIMO Wiener systems using subspace model identification methods. Signal Processing,
vol. 52(2), pp. 235–258, 1996.
[27] Glaria López, T.A. and Sbarbaro, D. Observer design for nonlinear processes with Wiener structure. Proceedings 50th IEEE Conf.
Decision and Control and European Control Conference, pp. 2211–2316, 2011.
[28] Blanchini, F. and Miani, S. Stabilization of LPV systems: State feedback, state estimation, and duality. SIAM Journal on Control and
Optimization, vol. 42, pp. 76–97, 2003.
[29] Louveaux, F.V. Piecewise convex programs. Mathematical Programming, vol. 15(1), pp. 53–62, 1978.
[30] Lau, K.K. and Womersley, R.S. Multistage quadratic stochastic programming. Journal of Computational and Applied Mathematics,
vol. 129(1-2), pp. 105–138, 2001.
[31] Patrinos, P., and Sarimveis, H. Convex parametric piecewise quadratic optimization: Theory and algorithms. Automatica, vol. 47(8),
pp. 1770–1777, 2011.
[32] Wu, L. Error bounds for piecewise convex quadratic programs and applications. SIAM Journal on Control and Optimization, vol.
33(5), pp. 1510–1529, 1995.

