title: Data-driven filtering for linear systems using Set Membership multistep predictors
authors: Marco Lauricella Lorenzo Fagiano
text: 
  
  
  
  
  
  
  Data-driven filtering for linear systems using
  Set Membership multistep predictors
      Marco Lauricella  and 
  	Lorenzo Fagiano
  	The authors are with the Dipartimento di Elettronica, Informazione e Bioingegneria, Politecnico di Milano, Piazza Leonardo da Vinci
  32, 20133 Milano, Italy (e-mail: marco.lauricella@polimi.it; lorenzo.fagiano@polimi.it). Corresponding author: Marco Lauricella.This is the extended version of an accepted paper submitted to the 2020 IEEE Conference on Decision and Control.
      
  =================================================================================================================================================================================================================================================================================================================================================================================================================================
  
  
  
  
  This paper presents a novel data-driven, direct filtering approach for unknown linear time-invariant systems affected by unknown-but-bounded measurement noise. The proposed technique combines independent multistep prediction models, identified resorting to the Set Membership framework, to refine a set that is guaranteed to contain the true system output. The filtered output is then computed as the central value in such a set. By doing so, the method achieves an accurate output filtering and provides tight and minimal error bounds with respect to the true system output. To attain these results, the online solution of linear programs is required. A modified filtering approach with lower online computational cost is also presented, obtained  by moving the solution of the optimization problems to an offline preliminary phase, at the cost of larger accuracy bounds. The performance of the proposed approaches are evaluated and compared with those of standard model-based filtering techniques in a numerical example. 
  
  
  
  
  § INTRODUCTION
  
  In this paper, we address the problem of output filtering for the case of linear time-invariant systems subject to unknown-but-bounded measurement disturbances. Our goal is to obtain an accurate filtering of the output from the noise, and, at the same time, to provide tight bounds on the error between the filtered output and the true one. 
  The problem of process variables filtering has been widely studied over the years, and the related literature includes a large number of contributions, see, e.g., <cit.> and the references therein. Among them, the most famous filtering approaches stem from the seminal works of Kolmogorov-Wiener <cit.>, and of Kalman-Bucy <cit.>. These works, as a vast majority of the filtering methods proposed during the years, mainly focus on the case of perfectly known system model, and stochastic disturbances with known probabilistic properties, for which they achieve optimal filtering properties. Other approaches assume that the accuracy is measured by suitable noise to error norm gains, where optimal performance are achieved by ℋ_2, ℓ_1, ℋ_∞ filters, see e.g. <cit.>. A different filtering methodology is also given by Moving Horizon Estimators, see e.g. <cit.>. Most of these filtering approaches are mainly based on the assumption of exact model knowledge as well. Since this is rarely the case in real world applications, where the system model is often unknown, a two step procedure is commonly adopted, where a model of the system at hand is identified from a noise-affected set of data collected from an experiment performed on the system, which is then used to define the filter. In this way, no optimality properties can be guaranteed anymore.
  
  Here, we resort to Set Membership (SM) identification methods to design a data-based, direct (i.e., without estimating a model of the plant) output filtering approach, able to obtain a `small' estimation error, measured by suitable guaranteed accuracy bounds. 
  SM identification methods usually operate under the unknown-but-bounded disturbance and uncertainty framework, see e.g., <cit.>, and can provide strong results in terms of optimality of the accuracy bounds in a worst-case sense, while requiring rather mild assumptions on the system and the disturbances <cit.>.
  
  In the SM literature, there are several contributions addressing the problem of data-driven filter design, e.g., <cit.>, where a filtered version of the variables of interest is provided, along with guaranteed bounds on their accuracy. Inspired by the results of a previous work (<cit.>), where the theoretical properties of the error bounds of multistep and onestep iterated models are investigated, here, we propose a novel filtering algorithm, based on the combined use of independent multistep prediction models, that is able to achieve tight accuracy bounds. The main idea is to derive a set containing the true system output at a given time, as the intersection of the accuracy regions provided by several SM multistep predictors, each one with a different prediction horizon. Then, the wanted filtered output is computed as the center of such a set. This quantity is locally optimal, in a worst-case error sense, i.e. it attains the minimal achievable error bounds under the considered assumptions. This is possible thanks to the use of independent multistep predictors, allowing one to achieve good prediction accuracy for different steps ahead, see e.g., <cit.>, and to derive less conservative bounds on the prediction accuracy, when compared to a single prediction model iterated in simulation <cit.>. Moreover, independent p-steps ahead models also allow one to achieve the theoretical minimum p-steps ahead prediction error bound, as defined under the SM framework, see <cit.>. By combining the output of different multistep predictors, all attaining their theoretical minimum error, we can further refine the set containing the system output, thus obtaining even tighter accuracy bounds.
  
  The approach requires the online solution of linear programs (LPs), which can be afforded in many real world applications of interest, thanks to the computational power available nowadays. For  cases where the available computational power is limited, and/or the sampling time is too short for the online solution of optimization problems, we also present a second filtering approach based on the same principle, which uses global accuracy bounds computed in a preliminary offline identification phase, at the cost of higher conservativeness.
  
  The filtering approaches proposed in this paper are finally showcased on a numerical example, where their performance are compared to those of a standard robust Kalman filter based on the Riccati equation approach.
  
  
  
  § PROBLEM STATEMENT AND PRELIMINARY RESULTS
  
  
  
  
   §.§ Problem statement
  
  Consider a discrete time, linear time invariant system, with input u(k)∈ℝ^m, and output z(k)∈ℝ, where k∈ℤ denotes the discrete time variable, modeled using an input-output auto-regressive with exogenous input (ARX) representation of order o (see <cit.> for details regarding the order choice)
  
      z(k+1)=ψ(k)^Tθ^0,
  
  where θ^0 is the vector of the real system parameters, and the regressor ψ(k) is defined as
  
      ψ(k)=   [ Z_o^T(k)   U_o^T(k) ]^T∈ℝ^o+mo,
      
      Z_o(k)=   [ z(k)   z(k-1)    z(k-o+1) ]^T∈ℝ^o,
      
      U_o(k)=   [u(k)^T   u(k-1)^T   u(k-o+1)^T ]^T∈ℝ^mo.
  
  The output measurement y(k)∈ℝ is affected by an additive noise d(k)∈ℝ:
  
      y(k)=z(k)+d(k).
  
  For the sake of notational simplicity, and without loss of generality, here we consider a multiple input, single output system. The proposed approach can be easily generalized to the multiple output case, by applying it to one output at the time. 
  
  Let us make the following assumption on the system at hand:
  
  The measurement noise and the system input are bounded. In particular:
  
  
    * |d(k)| ≤d̅_0, ∀ k∈ℤ,  d̅_0 ∈ℝ,  d̅_0>0.
  
    * u(k)∈𝕌⊂ℝ^m,  ∀ k ∈ℤ,  𝕌 compact.
  
  
  Assumption <ref> is common for system identification problems based on the framework of unknown-but-bounded disturbance, and it is valid for many practical applications as well.
  Here, we consider independent p-steps ahead prediction models in the ARX form
  
      ẑ(k|k-p)=φ_p(k-p)^Tθ_p,
  
  where θ_p is the vector of the predictor parameters to be identified, and the noise-affected regressor φ_p(k-p) is defined as
  
      φ_p(k-p)=   [ Y_o^T(k-p)   U_p,o^T(k-p) ]^T∈ℝ^o+m(o+p-1),
      
      Y_o(k-p)=   [ y(k-p)   y(k-p-1)    y(k-p-o+1) ]^T∈ℝ^o,
      
      U_p,o(k-p)=   [u(k-1)^T    u(k-p)^T    u(k-p-o+1)^T ]^T∈ℝ^m(o+p-1).
  
  Notice that a p-steps ahead model of the system can be readily obtained by recursion of (<ref>), leading to
  
      z(k|k-p)=ψ_p(k-p)^Tθ_p^0,
  
  where ψ_p(k-p) is the noise-free counterpart of φ_p(k-p), obtained substituting Y_o(k) with Z_o(k) in (<ref>), and the parameter values in θ_p^0 are polynomial functions of the entries of θ^0.
  
  Let us assume that an experiment is performed on the system at hand to collect a finite number of sampled data (ỹ(k),ũ(k)) to be used for the identification task. Here, ·̃ is used to denote a sampled and stored data point of a variable. For each prediction step p, these data form the following set, composed of N sampled regressors and N corresponding p-steps ahead output measurements:
  
      𝒱̃_p^N ≐{ṽ_p(k)=[     ; ỹ(k) ],   k=1,,N },
  
  where 𝒱̃_p^N ⊂ℝ^1+o+m(o+p-1). For simplicity and without loss of generality, we consider that the number N of data points is the same for any considered value of p. Moreover, as customary for identification problems based on sampled data, we consider the following assumption on the available dataset:
  
  As N→∞, the sampled dataset 𝒱̃_p^N is such that the set of all the system trajectories of interest is densely covered, i.e., the input u is persistently exciting and the disturbance d is bound-exploring.
  
  
  Prediction models (<ref>) are also known as multistep predictors, as they are designed to directly provide the output prediction p-steps ahead, and do not require the integration of an underlying simulation model. Note that, here, we consider prediction models having the same autoregressive order of (<ref>). In case the order o is not known, it can be estimated resorting to the procedure proposed in <cit.>. Moreover, in the remainder of this paper, we assume that either the disturbance bound d̅ is known, or that its estimate, obtained as described in <cit.>, is available.
  
  The problem addressed in this paper is the following: given the data set (<ref>), derive a filtering algorithm that returns an estimate ẑ(k)≈ z(k) of the system output, together with guaranteed bounds on the error |ẑ(k)- z(k)|.
  
  
  
   §.§ Preliminary results
  
  The error between the true system output z(k) of (<ref>) and the p-steps ahead prediction ẑ(k|k-p), given by (<ref>), originates from two sources:  the difference between the real parameter values θ^0 and the predictor ones, and the noise affecting the measurements used in the regressor . Under Assumption <ref>, this p-steps ahead error can be upper bounded by a worst-case error bound , which can be estimated resorting to the SM framework, see <cit.>, as
  
      [                           = α·min_θ_p,λ∈ℝ^+λ;                                   subject to; |ỹ-^T θ_p |≤λ +d̅,   ∀(,ỹ): [  ; ỹ ]∈𝒱̃_p^N, ]
  
  with α>1. For a fixed value of p,  represents the global error bound related to all possible p-steps ahead predictors of the form of (<ref>), and it is used to define the Feasible Parameter Set (FPS):
  
      Θ_p ≐{θ_p : |ỹ - ^T θ_p |≤+d̅,  ∀(,ỹ) : [  ; ỹ ]∈𝒱̃_p^N }.
  
  The FPS Θ_p is the set of all possible p-steps ahead predictor parameters that are consistent with the sampled data and the disturbance bound  d̅.
  Under Assumption <ref>, Θ_p is a polytope having at most 2N facets, represented in (<ref>) with an inequality description. If Θ_p happens to be unbounded, it is and indication that the data collected from the system are not informative enough and/or that N is not big enough, invalidating Assumption <ref>, thus new data should be acquired. Moreover, the bound  includes a factor α>1 (see (<ref>)) to account for the uncertainty due to the usage of a finite dataset. More details regarding the error bounds, the FPS, their theoretical properties, and the usage of scaling parameters can be found in <cit.>.
  
  
  
  § MULTISTEP FILTERING WITH LOCAL BOUNDS
  
  Under the SM framework, it is possible to associate guaranteed accuracy bounds to multistep prediction models in the form of (<ref>), resorting to the bounds , and the FPSs Θ_p defined in Section <ref>, see e.g., <cit.>. 
  
  The guaranteed accuracy bound  is defined, for a p-steps ahead predictor with parameters θ_p and for a given regressor φ̃_p(k-p), as
  
      =γ(max_θ∈Θ_p|φ̃_p(k-p)^T(θ-θ_p) |) + ,
  
  with γ>1. Here, γ has the same role of the scaling parameter α in (<ref>), thus accounting for the uncertainty due to the usage of a finite dataset. The bound (<ref>) is named “local”, since it pertains to a specific regressor value φ̃_p(k-p). Under Assumptions <ref> and <ref>, it holds that, by construction
  
      |z(k)-φ̃_p(k-p)^Tθ_p | ≤,
  
  i.e., the system output is guaranteed to lie inside a set defined by an interval centered at the p-steps ahead prediction ẑ(k|k-p)=φ̃_p(k-p)^Tθ_p, with amplitude equal to the corresponding accuracy bound.
  
  Resorting to this property, the novel data-driven filtering approach we propose is able to refine the true system output uncertainty range. This result is obtained by intersecting the uncertainty intervals of the outputs given by a group of p̅ multistep models  (<ref>), all providing a prediction of the output z(k), ranging from ẑ(k|k-1) to ẑ(k|k-p̅). Thus, the set Z_p̅(k) containing the true output at time k is derived as:
  
      z(k)∈ Z_p̅(k)={z̅   :  φ̃_p(k-p)^Tθ_p-τ̂_p(φ̃_p(k-p),θ_p) ≤z̅≤φ̃_p(k-p)^Tθ_p+τ̂_p(φ̃_p(k-p),θ_p), 
         ∀θ_p∈Θ_p, ∀ p=1,,p̅}.
  
  The boundaries of Z_p̅(k) are:
  
      z^max(k)   =max_z∈ Z_p̅(k) z = min_p=1,,p̅ζ_p^max, 
      
      z^min(k)   =min_z∈ Z_p̅(k) z = max_p=1,,p̅ζ_p^min,
  
  where
  
  
      ζ_p^max=min_θ_p∈Θ_p(φ̃_p(k-p)^Tθ_p+),
  
  
      ζ_p^min=max_θ_p∈Θ_p(φ̃_p(k-p)^Tθ_p-).
  
  
  Note that problems (<ref>) and (<ref>) can be recast as 2 Linear Programs (LP), which requires the previous solution of 2 LPs. Take for example the original optimization problem:
  
      ζ_p^max=min_θ_p∈Θ_p( φ̃_p(k-p)^Tθ_p + γmax_θ∈Θ_p|φ̃_p(k-p)^T (θ-θ_p) |+).
  
  This can be formulated as a minimization problem, where the absolute value is removed:
  
      ζ_p^max=   min_τ,θ_p∈Θ_pφ̃_p(k-p)^Tθ_p+ γτ+
         subject to
         max_θ∈Θ_p( φ̃_p(k-p)^T(θ-θ_p) )≤τ
         max_θ∈Θ_p( -φ̃_p(k-p)^T(θ-θ_p) )≤τ
  
  Then, the maximization problems in (<ref>) can be solved as independent LPs, since
  
      max_θ∈Θ_pφ̃_p(k-p)^T(θ-θ_p)=max_θ∈Θ_p(φ̃_p(k-p)^Tθ) - φ̃_p(k-p)^Tθ_p,
   leading to:
  
      ζ_p^max=   min_τ,θ_p∈Θ_p( φ̃_p(k-p)^Tθ_p+ γτ)+
         subject to
         c_1_p_k-φ̃_p(k-p)^Tθ_p≤τ
         c_2_p_k+φ̃_p(k-p)^Tθ_p≤τ
  
  where c_1_p_k=max_θ∈Θ_pφ̃_p(k-p)^Tθ, and c_2_p_k=max_θ∈Θ_p-φ̃_p(k-p)^Tθ. Similarly, problem (<ref>) becomes
  
      ζ_p^min=-   min_τ,θ_p∈Θ_p(-φ̃_p(k-p)^Tθ_p+ γτ)-
         subject to
         c_1_p_k-φ̃_p(k-p)^Tθ_p≤τ
         c_2_p_k+φ̃_p(k-p)^Tθ_p≤τ
  
  The filtered value of z(k) is then obtained using a central algorithm on the set Z_p̅(k), since the center of the local uncertainty interval attains the minimal achievable local accuracy bound, in a worst-case error sense, see <cit.>:
  
      ẑ_f_p̅(k)≐1/2(z^max(k)+z^min(k) ).
  
  In summary, the computation of (<ref>) implies the solution of 4p̅ LPs having 2N constraints at each sampling time, and requires the offline computation of the error bounds , for p∈[1,,p̅], obtained solving p̅ LPs. Since the FPSs are defined using only data pertaining to the identification dataset, and they are not updated online, it is possible to reduce the computational effort needed to solve (<ref>) and (<ref>) by performing an offline redundant constraints removal procedure on the FPSs, as proposed also in <cit.>. This produces a significant decrease in the number of constraints, reducing the computational effort needed for the online solution of the LPs. Examples of the complexity of Θ_p in terms of number of inequalities, and of computational times for the local approach, are provided in Section <ref>.
  
  Procedure <ref> summarizes the described local filtering algorithm.
  
  	Multistep filtering with local bounds
  	
  	
  		
    * Carry out the offline estimation of , and use it to define the FPSs for the considered steps p∈[1, p̅].
  		
    * For every time sample k, compute the online solution of the LPs (<ref>) and (<ref>), for p∈[1, p̅].
  		
    * Use the central algorithm (<ref>) to obtain the filtered system output, and compute its guaranteed accuracy bound for the given time sample.
  	
  
  
  The following result pertaining to the guaranteed accuracy of the local approach holds.
  
  	Considering the output filtering algorithm (<ref>), the following properties hold:
  	
  		
    * The guaranteed local accuracy bound τ_f_p̅(k) is
  		
      |ẑ_f_p̅(k)- z(k)|≤τ_f_p̅(k)=1/2| z^max(k)-z^min(k) |.
  
  		
    * The accuracy bound (<ref>) is the smallest worst-case error bound that can be achieved by any predicted output value.
  		
    * The accuracy bound (<ref>) is smaller than any single bound , ∀ p≤p̅.
  	
  
  
  	Straightforward consequence of (<ref>)-(<ref>), (<ref>). 
  
  
  
  The variable p̅, which represents the maximum prediction horizon length used in the filtering algorithm, is a tunable variable, whose choice corresponds to a trade-off between the tightness of the accuracy bound, and the resulting computational effort (which increases with p̅).
  
  
  
  
  § MULTISTEP FILTERING WITH GLOBAL BOUNDS
  
  The multistep filtering approach proposed in Section <ref> is able to achieve tight accuracy bounds, since it is based on the guaranteed local error bound (<ref>). To do so, the filter requires the online solution of 4p̅ LPs for each time sample. Even though this can be easily done in numerous real world applications, thanks to the computational power available nowadays, there are still a number of applications where the sampling time is too short, or the available computational power is limited, and the filtering approach based on local bounds would be inapplicable online. In all these cases, it is possible to define an output filtering algorithm similar to (<ref>), but based on a global version of the guaranteed accuracy bound. The main advantage is that fixed prediction models and their bounds are computed offline during the identification phase, together with  and the FPSs, and no optimization problem needs to be solved online. On the other side, the global version of bound τ_p(θ_p) is designed to hold for any possible regressor value, corresponding to any possible system trajectory, thus it is more conservative than its local counterpart, which is designed specifically for a given regressor instance. Moreover, the multistep predictors are identified only once, so their performance could clearly not be better than that of the locally optimal prediction models described in Section <ref>. To select such fixed predictors, a suitable optimality criterion needs to be defined. 
  
  Multistep prediction models can be identified resorting to several criteria, among which least-squares estimation is probably the most famous <cit.>. Since our goal here is to obtain tight uncertainty intervals for the unknown system output, even in the global case, it seems reasonable to seek for the multistep predictors that minimize the corresponding p-steps ahead global error bound, thus obtaining the theoretical minimum global uncertainty interval for the corresponding horizon length p. For given predictor parameters θ_p, the corresponding guaranteed global accuracy bound is estimated as:
  
      τ̅̂̅_p(θ_p)=γ̅( max_∈𝒱̃_p^Nmax_θ∈Θ_p|^T (θ - θ_p ) |) + ,  γ̅>1.
  
  Under assumption <ref>, it can be shown that the global bound  is such that <cit.>:
  
      |z(k)-φ̃_p(k-p)^Tθ_p | ≤τ̅̂̅_p(θ_p),
  
  i.e. it holds for any regressor belonging to the set of system trajectories covered by data (compare with (<ref>)). Then, we search inside the FPS Θ_p for the parameter vector θ̂_p that minimizes the global bound τ̅̂̅_p(θ̂_p), i.e. 
      θ̂_p=argmin_θ_p∈Θ_pτ̅̂̅_p(θ_p).
  
  This corresponds to
  
      θ̂_p = argmin_θ_p∈Θ_p max_k=p+1,,N max_θ∈Θ_p|φ̃_p(k-p)^T (θ-θ_p) |
  
  Problem (<ref>) can be reformulated as an LP, which requires the previous solution of 2(N-p) LPs, see <cit.>, as
  
      [                         θ̂_p = argmin_ζ,θ_p∈Θ_pζ;                                       subject to; c_k_p-φ̌_p(k)^Tθ_p≤ζ,   k∈ [p+1, N]∪[N+p+1, 2N], ]
  
  where c_k_p≐max_θ∈Θ_pφ̌_p(k)^T θ, and
  
      φ̌_p(k)=[    φ̃_p(k-p)           if         k≤ N; -φ̃_p(k-p-N)           if        k> N. ]
  
  The multistep model identified using (<ref>) is tailored for a specific horizon length p, for which it achieves optimal prediction properties, in terms of guaranteed global error bound amplitude. Then, we can intersect again the uncertainty intervals defined by the global accuracy bounds, obtaining the set Z^g_p̅(k) such that:
  
      z(k)∈ Z^g_p̅(k)   ={z̅: φ̃_p(k-p)^T-τ̅̂̅_p(θ̂_p)≤z̅≤φ̃_p(k-p)^T+τ̅̂̅_p(θ̂_p),  ∀ p=1,,p̅}.
  
  The filtered output value at time k is then given by the center of the set Z^g_p̅:
  
      ẑ^g_f_p̅(k)≐1/2(z^max_g(k)+z^min_g(k) ),
  
  where
  
      z^max_g(k)   =max_z∈ Z^g_p̅(k) z = min_p=1,,p̅ζ_g_p^max, 
      
      z^min_g(k)   =min_z∈ Z^g_p̅(k) z = max_p=1,,p̅ζ_g_p^min,
  
  and
  
      ζ_g_p^max   =φ̃_p(k-p)^T+, 
      ζ_g_p^min   =φ̃_p(k-p)^T-.
  
  The guaranteed accuracy bound of ẑ^g_f_p̅(k) is
  
      τ^g_f_p̅(k)=1/2| z_g^max(k)-z_g^min(k) |.
  
  The global filtering algorithm is sketched in Procedure <ref>.
  
  
  Multistep filtering with global bounds
  
  
  
    * Perform the offline estimation of , and define the FPSs, for p∈[1, p̅]. Then, perform the offline identification of p̅ independent p-steps ahead prediction models of the type of (<ref>) by solving (<ref>) for each p=1,…,p̅.
  
    * At each time k, compute online the boundaries of the uncertainty interval Z_p̅^g(k) using (<ref>).
  
    * Use the central algorithm (<ref>) to obtain the filtered system output, and compute its guaranteed accuracy bound for the given time sample with (<ref>).
  
  
  
  
  
  § SIMULATION RESULTS
  
  The performance of the proposed filtering approaches have been tested on a numerical example, and compared to those of the well known Kalman filter, see e.g., <cit.> for details. We consider a sampled, single input single output linear time-invariant system in continuous time, whose transfer function is
  
      G(s)=160/(s+10)(s^2+0.8s+16),
  
  where s is the Laplace variable. The input provided in the experiments is a three levels signal taking values in the set {-1, 0, 1} randomly every 4 time instants. We collected 12000 input-output samples with a sampling time T_s=0.1.
  
  
  
  
  
  
  We performed the numerical simulations considering three different disturbance scenarios: a) measurement disturbance signal given by a uniformly distributed random noise belonging to the interval [-0.2,   0.2]; b) measurement disturbance signal given by a Gaussian noise d(t)∼𝒩(0,σ_d^2), with σ_d^2=0.01; c) process and measurement disturbances given by the Gaussian noises w(t)∼𝒩(0,σ_w^2), and d(t)∼𝒩(0,σ_d^2) respectively, with σ_w^2=0.001, and σ_d^2=0.01. The measured output corresponds to ỹ(k)=z(k)+d(k), while, for scenario c), the system input is u(k)=ũ(k)+w(k), where ũ(k) is the measured input. The Kalman filter performance, used as a benchmark for comparison, are assessed under two different conditions: the first one is given by a Kalman filter based on a perfect knowledge of the system model, tuned to give the best filtering performance starting from unknown initial conditions, and using the noise affected measurements collected under scenarios a)-c); the second setup is given by a Kalman filter based on a model of the system identified from the noise affected measurements resorting to the simulation error method, see e.g. <cit.>, tuned to give the best filtering performance starting from unknown initial conditions. In both cases, we adopt the robust Kalman filter based on the Riccati equation approach.
  
  
  The first half of the available dataset is then used to perform the identification of the bound , of the FPSs, of the multistep predictors, and of their guaranteed accuracy bounds , while the second half is used in validation to compute the filtered output using the local and the global filtering approaches, given by (<ref>) and (<ref>) respectively, and to obtain the corresponding uncertainty bounds (<ref>) and (<ref>). By applying the redundant constraints removal procedure, the number of FPS constraints is reduced from 12000 to an average of 331, for each Θ_p, with p∈[1, 60], with a minimum of 48, and a maximum of 598 constraints. Fig. <ref> depicts the estimated global accuracy bounds  for the independent multistep prediction models, obtained with α=1.2 and γ=1.1, for the three disturbance scenarios considered here.
  
  Table <ref> reports the guaranteed accuracy bounds, the Root Mean Squared Error (RMSE), defined as
  
      RMSE=√(∑_k=1^N ( z(k)-ẑ_f(k) )^2/N),
  
  and the maximum filtering error, given by max_k e(k), with e(k)=z(k)-ẑ_f(k), for the Kalman filter, and for the local and global output filtering approaches, considering different values of the prediction horizon p̅, for the case of uniformly distributed measurement noise. The accuracy bounds of the filtering approaches are then compared to that of the p-steps ahead model achieving the lowest  for each considered prediction horizon length. Tables <ref> and <ref> show the same performance comparison for the case of Gaussian measurement noise, and for the case of Gaussian process and measurement noises, respectively.
  
  
  
  These results show that the local filtering approach achieves the lowest filtering error, both in terms of RMSE and maximum error, and the lowest accuracy bounds, both on average and in worst-case, with respect to the global filtering approach, and to the usage of a single p-steps ahead prediction model. At the same time, the local filtering approach achieves good performances in terms of RMSE also when compared to the Kalman filter. It is worth to point out that, for p̅≥ 15, the local filtering algorithm outperforms the Kalman filter based on an estimated system model, both in terms of RMSE and of maximum filtering error, for all the disturbance scenarios considered here. Moreover, for the case of process and measurement Gaussian disturbances, the local filtering approach is able to achieve a RMSE comparable to that of the Kalman filter based on the exact system model, and much smaller maximum error values. Clearly, in all the cases where a process disturbance signal is not present, the Kalman filter based on the exact system model achieves almost zero filtering error, since the input and the system model are perfectly known, and the only source of uncertainty is given by the unknown initial conditions.
  
  The global filtering approach achieves intermediate results, closer to that obtained using a single multistep prediction model, with the advantage of a considerable reduction of the necessary online computational effort over the local approach. Table <ref> shows a brief comparison of the computational time required by the local and global approaches to provide the filtered version of one output sample. The simulations were performed on a Laptop equipped with Intel i7 dual-core processor with 2.4 GHz clock speed and 8 GB of RAM, using MatLab R2018b  solver, based on the Dual Simplex algorithm. Fig. <ref> depicts the filtered output obtained using the local and global multistep filtering approaches, together with their accuracy bounds, for different prediction horizon lengths, and compare them to the real system output z, and to its measurement ỹ. Fig. <ref> depicts an example of the true output uncertainty intervals, obtained as intersection of the multistep predictors accuracy regions, for the case of uniform measurement noise. Finally, Fig. <ref> reports the comparison between the filtered system output obtained using the local filtering approach, and that obtained from the Kalman filter based on the exact model of the system, and from the Kalman filter based on an identified system model, for the case of Gaussian process and measurement noise. 
  
  
  
  
  
  
  
  
  
  
  
  § CONCLUSIONS
  
  The data-driven direct filtering approach presented in this paper allows one to address the problem of filtering the output of linear time-invariant systems subject to unknown-but-bounded uncertainties without deriving a model of the system. The proposed filtering algorithm is able to achieve good filtering accuracy, quantified in terms of average filtering error, and of guaranteed accuracy bounds, and it allows one to compute tight guaranteed uncertainty intervals for the  true system output. The proposed method relies on the intersection of the uncertainty regions of different SM multistep predictors, which are used together to provide a filtered version of the system output, while refining the output uncertainty intervals thanks to the combined use of different prediction horizon lengths. The local filtering algorithm, and its accuracy bounds, are then obtained by means of linear programming, which makes them suitable for online implementation in real world applications. A more conservative global filtering approach was presented to lower the online computational effort, moving the solution of the optimization problems to an offline phase. Numerical simulations illustrate the performance of the proposed local filtering approach, which is able to outperform a Kalman filter based on an estimated system model, and to score similar performance to that of a Kalman filter based on the exact system model, while retaining the ability to provide optimal guaranteed accuracy bounds for the filtered output.
  
  
  IEEEtran
  10
  
  gelb1974applied
  A. Gelb, Applied optimal estimation.1em plus 0.5em minus
    0.4emMIT press, 1974.
  
  anderson2012optimal
  B. D. Anderson and J. B. Moore, Optimal filtering.1em plus 0.5em
    minus 0.4emCourier Corporation, 2012.
  
  saberi2007filtering
  A. Saberi, A. A. Stoorvogel, and P. Sannuti, Filtering theory: with
    applications to fault detection, isolation, and estimation.1em plus
    0.5em minus 0.4emSpringer Science & Business Media, 2007.
  
  kolmogorov1941stationary
  A. N. Kolmogorov, “Stationary sequences in Hilbert space,” Bull.
    Math. Univ. Moscow, vol. 2, no. 6, pp. 1–40, 1941.
  
  wiener1949smoothing
  N. Wiener, “Extrapolation, interpolation, and smoothing of stationary time
    series with engineering applications,” Wiley, New York, 1949.
  
  kalman1961new
  R. E. Kalman and R. S. Bucy, “New results in linear filtering and prediction
    theory,” Trans. ASME, D, vol. 83, pp. 95–108, 1961.
  
  Geromel2001
  J. Geromel and M. D. Oliveira, “ℋ_2 and ℋ_∞
    robust filtering for convex bounded uncertain systems,” IEEE
    Transactions on Automatic Control, vol. 46, no. 1, pp. 100–107, 2001.
  
  Sun2005
  K. Sun and A. Packard, “Robust ℋ_2 and ℋ_∞
    filters for uncertain LFT systems,” IEEE Transactions on Automatic
    Control, vol. 50, no. 5, pp. 715–720, may 2005.
  
  Nagpal1991
  K. Nagpal and P. Khargonekar, “Filtering and smoothing in an
    ℋ_∞ setting,” IEEE Transactions on Automatic
    Control, vol. 36, no. 2, pp. 152–166, 1991.
  
  Abedor1996
  J. Abedor, K. Nagpal, and K. Poolla, “A linear matrix inequality approach to
    peak-to-peak gain minimization,” International Journal of Robust and
    Nonlinear Control, vol. 6, no. 9-10, pp. 899–927, nov 1996.
  
  Voulgaris1995
  P. Voulgaris, “On optimal ℓ_∞ to ℓ_∞ filtering,”
    Automatica, vol. 31, no. 3, pp. 489–495, mar 1995.
  
  Colaneri2002
  P. Colaneri and A. Ferrante, “A J-spectral factorization approach for
    ℋ_∞ estimation problems in discrete time,” IEEE
    Transactions on Automatic Control, vol. 47, no. 12, pp. 2108–2113, dec
    2002.
  
  Shaked1992Hinf
  U. Shaked and Y. Theodor, “ℋ_∞-optimal estimation: a
    tutorial,” in [1992] Proceedings of the 31st IEEE Conference on
    Decision and Control.1em plus 0.5em minus 0.4emIEEE, 1992.
  
  alessandri2008moving
  A. Alessandri, M. Baglietto, and G. Battistelli, “Moving-horizon state
    estimation for nonlinear discrete-time systems: New stability results and
    approximation schemes,” Automatica, vol. 44, no. 7, pp. 1753–1765,
    2008.
  
  rao2003constrained
  C. V. Rao, J. B. Rawlings, and D. Q. Mayne, “Constrained state estimation for
    nonlinear discrete-time systems: Stability and moving horizon
    approximations,” IEEE transactions on automatic control, vol. 48,
    no. 2, pp. 246–258, 2003.
  
  milanese2013bounding
  M. Milanese, J. Norton, H. Piet-Lahanier, and E. Walter, Bounding
    approaches to system identification.1em plus 0.5em minus 0.4emSpringer Science & Business Media, 2013.
  
  walter1990estimation
  E. Walter and H. Piet-Lahanier, “Estimation of parameter bounds from
    bounded-error data: a survey,” Mathematics and Computers in
    simulation, vol. 32, no. 5-6, pp. 449–468, 1990.
  
  garulli1999robustness
  A. Garulli, A. Tesi, and A. Vicino, Robustness in identification and
    control.1em plus 0.5em minus 0.4emSpringer, 1999, vol. 245.
  
  milanese2011unified
  M. Milanese and C. Novara, “Unified set membership theory for identification,
    prediction and filtering of nonlinear systems,” Automatica, vol. 47,
    no. 10, pp. 2141–2151, 2011.
  
  milanese2010direct
  M. Milanese, F. Ruiz, and M. Taragna, “Direct data-driven filter design for
    uncertain LTI systems with bounded noise,” Automatica, vol. 46,
    no. 11, pp. 1773–1784, 2010.
  
  yang2009set
  F. Yang and Y. Li, “Set-membership filtering for systems with sensor
    saturation,” Automatica, vol. 45, no. 8, pp. 1896–1902, 2009.
  
  garulli1997conditional
  A. Garulli, A. Vicino, and G. Zappa, “Conditional central algorithms for
    worst-case estimation and filtering,” in Proceedings of the 36th IEEE
    Conference on Decision and Control, vol. 3.1em plus 0.5em minus
    0.4emIEEE, 1997, pp. 2453–2458.
  
  lauricella2020set
  M. Lauricella and L. Fagiano, “Set membership identification of linear systems
    with guaranteed simulation accuracy,” IEEE Transactions on Automatic
    Control, 2020.
  
  shook1991identification
  D. S. Shook, C. Mohtadi, and S. L. Shah, “Identification for long-range
    predictive control,” in IEE Proceedings D-Control Theory and
    Applications, vol. 138.1em plus 0.5em minus 0.4emIET, 1991,
    pp. 75–84.
  
  haber2003long
  R. Haber, U. Schmitz, and R. Bars, “Long-range optimal model and
    multi-step-ahead prediction identification for predictive control,”
    IFAC Proceedings Volumes, vol. 36, no. 16, pp. 477–482, 2003.
  
  farina2011simulation
  M. Farina and L. Piroddi, “Simulation error minimization identification based
    on multi-stage prediction,” International Journal of Adaptive Control
    and Signal Processing, vol. 25, no. 5, pp. 389–406, 2011.
  
  lauri2010pls
  D. Laurí, M. Martínez, J. V. Salcedo, and J. Sanchis, “PLS-based
    model predictive control relevant identification: PLS-PH algorithm,”
    Chemometrics and Intelligent Laboratory Systems, vol. 100, no. 2, pp.
    118–126, 2010.
  
  shook1992control
  D. S. Shook, C. Mohtadi, and S. L. Shah, “A control-relevant identification
    strategy for GPC,” IEEE Transactions on Automatic Control, vol. 37,
    no. 7, pp. 975–980, 1992.
  
  sjoberg1995nonlinear
  J. Sjöberg, Q. Zhang, L. Ljung, A. Benveniste, B. Delyon, P. Glorennec,
    H. Hjalmarsson, and A. Juditsky, “Nonlinear black-box modeling in system
    identification: a unified overview,” Automatica, vol. 31, no. 12, pp.
    1691–1724, 1995.
  
  terzi2019learning
  E. Terzi, L. Fagiano, M. Farina, and R. Scattolini, “Learning-based predictive
    control for linear systems: a unitary approach,” Automatica, vol.
    108, p. 108473, 2019.
  
  traub1988information
  J. Traub, G. Wasilkowski, and H. Wozniakowski, Information-Based
    Complexity.1em plus 0.5em minus 0.4emAcademic Press, New York,
    1988.
  
  soderstrom1989system
  T. Söderström and P. Stoica, System identification.1em
    plus 0.5em minus 0.4emPrentice-Hall, 1989.
  
  ljung1999system
  L. Ljung, System identification: theory for the user (2nd ed.).1em plus 0.5em minus 0.4emPTR Prentice Hall, 1999.
  
  tomita1992equation
  Y. Tomita, A. Damen, and P. V. D. Hof, “Equation error versus output error
    methods,” Ergonomics, vol. 35, no. 5-6, pp. 551–564, 1992.
  
  
  