On multi-step prediction models for receding horizon control
Enrico Terzi Lorenzo Fagiano Marcello Farina Riccardo Scattolini

	
	
	
	
	

On multi-step prediction models for receding horizon control
    Enrico Terzi, Lorenzo Fagiano, Marcello Farina, Riccardo Scattolini The authors are with the Dipartimento di Elettronica, Informazione e Bioingegneria, Politecnico di Milano, Milano, Italy. E-mail addresses:
			{enrico.terzi| lorenzo.fagiano |marcello.farina| riccardo.scattolini}@polimi.it 
    
=======================================================================================================================================================================================================================================================================================================






§ INTRODUCTION

This manuscript contains technical details of recent results developed by the authors on learning-based model predictive control for  linear time invariant systems.




§ PROBLEM FORMULATION


Let us consider a single-input, single-output (SISO), open-loop stable, discrete-time, strictly proper linear time invariant (LTI) system with n states, input u(k)∈ℝ and output z(k)∈ℝ, where k∈ℤ is the discrete time variable. The measurement of the system output y(k) is affected by a bounded additive disturbance d(k):

    y(k)=z(k)+d(k).

 (Measurements)
	
		
  * |d(k)|≤d, ∀ k∈ℤ.
		
  * The input u(k) is measured with negligible noise.□
	 



 (Input bounds) 
The input u belongs to a compact set: u(k) ∈𝕌⊂ℝ, ∀ k∈ℤ. □



Let us denote with p∈ℕ a finite number of time steps in the future. We are interested in deriving a prediction model of the future output z(k+p), exploiting the input and output measurements collected in the time interval [k-o+1,k], where o∈ℕ is the chosen order of the model, and the future (planned) inputs in the interval [k,k+p-1].  Specifically, at any instant k let us define the regressor φ_p(k)∈ℝ^2o-1+p as: 

    φ_p(k)≐[ Y_o^T(k)    U_o^T(k)    U^T_p(k) ]^T

where ^T is the matrix transpose operation and

    Y_o(k)    ≐[
    [     y(k);   y(k-1);        ⋮; y(k-o+1);          ]]               
    U_o(k)≐[
    [   u(k-1);        ⋮; u(k-o+1);          ]] 
    
    U_p(k)    ≐[
    [     u(k);   u(k+1);        ⋮; u(k+p-1);          ]].

Then, we consider the following linear model structure:

    ẑ(k+p)=φ_p(k)^Tθ_p

where ẑ(k+p) is the predicted output and θ_p ∈ℝ^2o-1+p is the model parameter vector. We refer to models of the form (<ref>) as “multi-step”, since for each p∈ℕ the corresponding prediction model provides directly an estimate of z(k+p). This is different from the (most common) approach of deriving a one-step-ahead model and then iterating it p times to compute predictions for the subsequent future time-steps.

In addition to the prediction model, we also want to derive guaranteed bounds on its accuracy. In particular, we aim for the following global guaranteed accuracy bound:

    |z(k+p)-ẑ(k+p)|≤τ_p(θ_p).

The value τ_p(θ_p) is termed “global” since it holds for any value of the regressor φ_p within a specified set, as we further detail in the remainder of this paper.

A collection of multi-step models derived for all p∈[1,p] provides an estimated sequence of future system outputs, up to the prediction horizon p<∞, together with an associated sequence of guaranteed uncertainty intervals τ_p(θ_p). These models can be then embedded in a robust finite-horizon optimal control problem to be solved in a receding-horizon approach, thus realizing a MPC law based on multi-step predictions <cit.>. Three important reasons to consider multi-step models are:

	
  * the model identification procedure results in convex optimization problems, as opposed to the nonlinear programs arising when identifying one-step-ahead models with an output-error (i.e. simulation) criterion <cit.>;
	
  * the guaranteed bound τ_p(θ_p) pertaining to a multi-step model is less conservative than the one pertaining to a one-step-ahead model iterated p times;
	
  * by considering an independent model for each value of p, there is no need for trade-offs between model accuracy at high-frequency (i.e. short prediction horizon) vs. low-frequency (i.e. long prediction horizon) that arise when choosing the simulation horizon in the identification of one-step-ahead models with an output-error criterion <cit.>. 

The use of a model structure that is linear in the parameters is justified both by point 1) above and by the observation that, in the case of zero measurement noise, the true output z(k+p) is indeed a linear function of the regressor, provided that the following assumption holds.


 (Observability, reachability and model order) 
	The system at hand is completely observable and reachable, and o ≥ n. □



Under Assumption <ref>, with straightforward manipulations one can show that the state at time k is in general a linear function of the past n input-output values, and that the output at time k+p is linear in the state and in the planned input values, hence resulting in the model structure (<ref>). From the practical standpoint, Assumption <ref> can be relaxed to account only for the observable and controllable sub-space of the system state. Moreover, this assumption can be satisfied by estimating the system order n (e.g. based on physical considerations) and/or by deriving models with growing order o and by monitoring the magnitude of the corresponding accuracy bounds, as we describe more in detail in section <ref>.

Assumption <ref>, combined with the asymptotic stability of the system under study, results in bounded sets where the regressor φ_p evolves in time. Specifically, for a given horizon p we consider a compact set Φ_p containing the regressor values of interest:

    φ_p(k)∈Φ_p⊂R^2o-1+p, Φ_p compact, ∀ p∈ℕ, ∀ k∈ℤ.

Φ_p needs not to be known explicitly and is in general a complicated set that depends on the system input-output trajectories. Rather, we assume that for a given value of p a finite batch of experimental data is available:

    (φ̃_p(i), ỹ_p(i)), i=1,…,N_p.

In (<ref>), the notation ·̃ indicates specific measured values of a quantity. φ̃_p(i) are the measured instances of the regressor φ_p∈Φ_p, and  ỹ_p(i) the corresponding measured values of noise-corrupted outputs p steps in the future, i.e. ỹ_p(i)=z(i+p)+d(i+p).  N_p<∞ is the total number of available pairs (φ̃_p(i), ỹ_p(i)) composing the data-set. These pairs can be easily built from a given data-set of measured input-output values, collected e.g. in a preliminary experiment on the system. We further define the vectors of sampled variables ṽ_p(i)≐[φ̃_p(i)^T ỹ_p(i)]^T, i=1,…,N_p and the corresponding countable set containing them:

    𝒱̃_p^N_p={ṽ_p(i) = [[ φ̃_p(i);  ỹ_p(i) ]],  i=1,… N_p}⊂ℝ^2o+p.

For any given value of φ_p∈Φ_p, the corresponding measured system output y_p is not uniquely determined a priori, rather it belongs to a set Y_p(φ_p)⊂R, due to the measurement noise d (both in the regressor and in the corresponding output). In view of Assumptions <ref> and <ref> and of the compactness of Φ_p, the set Y_p(φ_p) is also compact for any φ_p∈Φ_p. We can then define the following continuous counterpart of the set 𝒱̃_p^N_p:

    𝒱_p={v_p ≐[[ φ_p; y_p ]]: y_p∈Y_p(φ_p), ∀φ_p∈Φ_p}⊂ℝ^2o+p.

Namely, the set 𝒱_p contains all possible regressors φ_p in the compact Φ_p and, for each value of  φ_p , all possible output values in the corresponding compact set  Y_p(φ_p). We consider the following assumption linking the sets 𝒱̃_p^N_p and 𝒱_p.
 (Data-set) 
	For any β>0, there exists a value of N_p<∞ such that:
	
    d_2(𝒱_p,𝒱̃_p^N_p)≤β

	where d_2(𝒱_p,𝒱̃_p^N_p)≐max_v∈𝒱_pmin_u∈𝒱̃_p^N_pu-v_2 is the distance between sets 𝒱_p and 𝒱̃_p^N_p. □

Assumption <ref> is equivalent to assuming that lim_N_p→∞d_2(𝒱_p,𝒱̃_p^N_p)=0, i.e. that by adding more points to the measured data set, the underlying set of all trajectories of interest is densely covered. This is essentially an assumption on the persistence of excitation of the inputs used for the preliminary experiments, together with an assumption of bound-exploring property of the additive disturbance d, such that the bound d in Assumption <ref> is actually tight. 

We are now in position to formulate the problem that we will address in the remainder of this paper:


Problem 1 Under Assumptions <ref>-<ref>, for a given p<∞, p∈ℕ and any p=1,…,p, use the available data (<ref>) to identify a multi-step model of the form (<ref>) and estimate the associated guaranteed bounds (<ref>).□



We propose in the next section an approach to solve Problem 1, based on a Set Membership (SM) identification methodology <cit.>, <cit.>, which guarantees convergence of the derived error bounds to suitably defined optimal values.



§ LEARNING MULTI-STEP PREDICTION MODELS:
A SET MEMBERSHIP APPROACH
 

For the sake of simplicity, in the following we consider a single value of p∈[1,p], the extension to any other value is straightforward. The proposed approach consists of the following steps.

	
  * Define an optimality criterion to evaluate the model estimates, corresponding to an optimal (i.e. minimal) error bound.
	
  * Derive a procedure to estimate the optimal error bound.
	
  * Based on the available data and the error bound estimate, build the set of all parameters that are consistent with this information (Feasible Parameter Set, FPS).
	
  * Using the information summarized in the FPS, for any given model of the form (<ref>) compute the related guaranteed error bound τ_p.
	
  * Select a nominal model with minimal guaranteed error bound.

We describe next each step in detail, followed by a discussion on tuning and extensions of the approach.



 §.§ Optimal parameters and optimal error bound


Assume that a value of θ_p has been fixed. Then, under the considered working assumptions, ∀ k∈ℤ:

    y(k+p)   =z(k+p)+d(k+p)
       =ẑ(k+p)+ϵ_p(θ_p,φ_p(k))+d(k+p)
       =φ_p(k)^Tθ_p+ϵ_p(θ_p,φ_p(k))+d(k+p)

where ϵ_p(θ_p,φ_p(k)) is the error between the true system output and the estimated one:

    ϵ_p(θ_p,φ_p(k))=z(k+p)-φ_p(k)^Tθ_p

The quantity ϵ_p(θ_p,φ_p(k)) accounts for the quality of the chosen parameter values, for the model order mismatch (o>n, compare Assumption <ref>) and for the noise in the regressor measurements. Since the underlying system dynamics are time-invariant,  ϵ_p depends inherently on the model parameter vector θ_p and on the specific  regressor φ_p(k). ϵ_p(θ_p,φ_p(k)) is bounded because both  y(k+p) and θ_p^Tφ_p(k) are, due to the stability of the system and compactness of the set containing the input values. From (<ref>) and Assumption <ref> we have:

    [  |y(k+p)-φ_p(k)^Tθ_p| ≤ |ϵ_p(θ_p,φ_p(k))|+d;                                ≤ϵ̅_p(θ_p)+d ]

where ϵ̅_p(θ_p) is the global error bound with respect to all possible regressors of interest in the compact Φ_p:

    [                                         ϵ̅_p(θ_p)=                                          min_ϵ∈ℝ ϵ;                                                                                            subject to;                                                    |y_p-φ_p^Tθ_p|≤ϵ+d, ∀ (φ_p,y_p):[[ φ_p; y_p ]]∈𝒱_p ]

The quantity ϵ̅_p(θ_p) is the tightest bound on the global (i.e. worst-case) estimation error that a given parameter vector  θ_p can produce. We can now define the  optimal parameter values (i.e. optimal models) as those that minimize such a bound. As a technical assumption, we consider all the parameters within a compact set Ω⊂ℝ^2o-1+p. Ω can take into account application-specific prior information on the model parameters or, if no such information is available, it can be chosen as a large-enough set (e.g. by considering box constraints of ±10^15 on each element of the parameter vector). This assumption allows us to use maximum and minimum operators instead of supremum and infimum. The set Θ_p^0 of optimal parameter values is the following:

    Θ_p^0={θ_p^0: θ_p^0=min_θ_p ∈Ω ϵ̅_p(θ_p)},
 
and we denote with ϵ̅_p^0 the corresponding optimal error bound:

    ϵ̅_p^0=min_θ_p ∈Ω ϵ̅_p(θ_p).
 
Considering (<ref>)-(<ref>) we can alternatively write:

    Θ_p^0={θ_p : |y_p-φ_p^Tθ_p|≤ϵ̅_p^0+d, ∀ (φ_p,y_p):[[ φ_p; y_p ]]∈𝒱_p}.
 



 §.§ Estimating the optimal error bound


The optimal models and optimal error bound cannot be computed in practice, since the solution to (<ref>) would imply the availability of an infinite number of data and the solution to an infinite-dimensional optimization program. However, we can compute an estimate λ_p≈ϵ̅_p^0 from the available experimental data, by solving the following linear program (LP):

    [                                          λ_p=min_θ_p,λ λ;                                               subject to; |ỹ_p-φ̃_p^Tθ_p|≤λ+d, ∀ (φ_p,y):[[ φ̃_p;  ỹ_p ]]∈𝒱̃_p^N_p;                                                     λ≥ 0 ]

The last inequality constraint in (<ref>) is required to enforce a positive estimate of the error bound: without this constraint, the obtained estimate could result to be negative, especially in presence of small amount of data and output disturbance realizations with much smaller magnitude than the considered bound d.

The following result shows that, under the considered working assumptions, the estimate (<ref>) converges to the optimal one, ϵ̅_p^0.



	Let Assumptions <ref>-<ref> hold. Then:
	
		
  * λ_p≤ϵ̅_p^0;
		
  * ∀ρ∈(0,ϵ̅_p^0] ∃  N_p<∞  : λ_p≥ϵ̅_p^0-ρ□
	
	




	See the appendix.

Theorem <ref> implies that lim_N_p→∞(ϵ̅_p^0-λ_p)=0^+, i.e. that the estimate (<ref>) converges to the optimum from below. This is unavoidable with the considered problem settings, where limited information is available. In practice, one can increase the number N_p of experimental data and observe the behavior of λ_p, which typically (if the preliminary experiments are informative enough) quickly converges to a limit. Then, a practical approach to compensate for the uncertainty caused by the use of a finite number of measurements  is to inflate the value λ_p:

    ϵ̂̅̂_p=αλ_p, α>1.

With sufficiently large N_p, the coefficient α can be chosen very close to 1. We show an example of such a procedure in section <ref>, and provide more comments on the tuning of α in section <ref>. We consider the following assumption in the remainder of this paper:


 (Estimate of the optimal error bound)

	ϵ̂̅̂_p≥ϵ̅_p^0.□






 §.§ Feasible Parameter Set


We can now exploit the estimated optimal error bound to construct the tightest set of parameter values that are consistent with all the prior information, i.e. the FPS Θ_p:

    Θ_p={θ_p : |ỹ_p-φ̃_p^Tθ_p|≤ϵ̂̅̂_p+d, ∀ (φ̃_p,ỹ_p):[[ φ̃_p;  ỹ_p ]]∈𝒱̃_p^N_p}
 

The set Θ_p is non-empty by construction, since under Assumption <ref> we have (compare (<ref>) and (<ref>)):

    Θ_p^0⊆Θ_p.
 
If the FPS is bounded, it results in a polytope with at most N_p faces. If it is unbounded, then this is a sign that the employed measured data are not informative enough to derive a bound on the worst-case model error (as we show next) and that the number of available data N_p shall be increased until a bounded FPS is obtained. This situation usually occurs when very few data points are used (e.g. N_p<2o-1+p) or the preliminary experiments are not informative enough (e.g. when only steady-state data are used).




 §.§ Error bound computation for a generic model

Having defined the FPS, we can proceed to derive a bound on the prediction error achieved by any  model of the form (<ref>). Let us consider a generic value of θ_p to derive the prediction model (<ref>). Then, considering any θ_p^0∈Θ_p^0 and any φ_p(k)∈Φ_p, using (<ref>) and (<ref>) we have:

    [                                                                               | z(k+p) - ẑ(k+p) |;                                                = |φ_p(k)^Tθ_p^0 - φ_p(k)^Tθ_p + ϵ(θ_p^0,φ_p(k)) |;                                                ≤                     |φ_p(k)^T(θ^0_p-θ_p)|+ϵ̅_p^0 ]

The tightest bound we can derive on (<ref>) is based on the knowledge that θ^0_p∈Θ_p (see (<ref>)) and that ϵ̅_p^0≤ϵ̂̅̂_p (Assumption <ref>):

    |φ_p(k)^T(θ^0_p-θ_p)|+ϵ̅_p^0≤max_θ∈Θ_p|φ_p(k)^T(θ-θ_p) | + ϵ̂̅̂_p.

The latter bound is the tightest local error bound (i.e. valid for a given value of φ_p) for the model given by the considered parameter value θ_p. As stated in Problem 1, for the sake of using the model within a robust MPC framework, we want to derive a global bound τ_p(θ_p) holding ∀φ_p∈Φ_p. Considering (<ref>)-(<ref>) leads to:

    [                                                            max_φ_p∈Φ_p| z(k+p) - ẑ(k+p) |;                                            ≤ max_φ_p∈Φ_pmax_θ∈Θ_p|φ_p^T(θ-θ_p) | + ϵ̂̅̂_p;                                            =                                     τ_p(θ_p) ]

Such a bound cannot be derived exactly with finite data under the considered assumptions, and its computation would be intractable also if the set Φ_p were known precisely (unless some additional assumption is made, e.g. polytopic set Φ_p). However, we can approximate it by computing the maximum of (<ref>) over the data-set 𝒱̃_p^N_p:

    τ_p(θ_p)=max_i=1,…,N_pmax_θ∈Θ_p|φ̃_p(i)^T(θ-θ_p) | + ϵ̂̅̂_p

The following result shows convergence of τ_p(θ_p) to τ_p(θ_p).

	Let Assumptions <ref>-<ref> hold. Then, for any θ_p∈ℝ^2o-1+p:

	
  * τ_p(θ_p)≤τ_p(θ_p);
	
  * ∀ρ∈(0,τ_p(θ_p)] ∃  N_p<∞  : τ_p(θ_p)≥τ_p(θ_p)-ρ□






Straightforward extension of the proof of Theorem <ref>.

Considerations similar to those of Theorem <ref> hold also for the bound τ_p(θ_p): it is possible to monitor its behavior with increasing N_p value in order to evaluate its convergence. As done in (<ref>), we inflate this bound to account for the uncertainty deriving from our finite data-set:

    τ̂_p(θ_p)=γτ_p(θ_p), γ>1,

and we assume that the resulting estimate is larger than the true bound:
 (Error bound estimate for a given θ_p)

	τ̂_p(θ_p)≥τ_p(θ_p), ∀θ_p∈ℝ^2o-1+p.□




 §.§ Nominal model selection

The last step in the proposed approach is to select a nominal model. The most common approach is probably least-squares estimation, in which case the results of section (<ref>) can be anyway applied to obtain an estimate of the resulting global error bound. Since the final goal is to employ the model in a MPC algorithm, we rather seek the model that minimizes the uncertainty bound. Specifically, considering that the tightest set that contains the optimal parameter values (i.e. with minimum error, see section (<ref>)) is the FPS Θ_p, we search within this set for the parameter value that minimizes the resulting bound τ̂_p(θ_p):

    θ_p^*=min_θ_p∈Θ_pτ̂_p(θ_p).

The resulting nominal model is ẑ(k+p)=φ_p(k)^Tθ_p^*, and the associated error bound is:

    τ̂_p(θ_p^*)=γ(min_θ_p∈Θ_pmax_i=1,…,N_pmax_θ∈Θ_p|φ̃_p(i)^T(θ-θ_p) |) + ϵ̂̅̂_p.

Note that term ϵ̂̅̂_p in (<ref>) does not depend on θ^*_p and it converges to the optimal error bound ϵ̅_p^0 as N_p increases (Theorem <ref>). Moreover, it can be shown that, under the considered working assumptions, if θ_p^0 is unique then the difference |θ_p^0-θ_p^*| tends to zero, and the associated error bound (<ref>) tends to the minimum value ϵ̅_p^0 as well. Furthermore, for any value of N_p it can be shown that τ̂_p(θ_p^*) corresponds to the radius of information, i.e. the minimum guaranteed error value that can be attained with the given prior information and data <cit.>. Finally, it also holds that, for each p>1 the bounds τ̂_p(θ_p^*) computed in our approach are less conservative than the worst-case bounds obtained by iterating any 1-step-ahead prediction model of the form (<ref>). All these derivations are omitted here for the sake of brevity. 
  The optimization problem to be solved to compute θ_p^* (<ref>) takes the form 

    θ_p^*=min_θ_p∈Θ_pmax_i=1,…,N_pmax_θ∈Θ_p|φ̃_p(i)^T(θ-θ_p) |.

This problem can be solved by reformulating it as 2 N_p+1 LPs <cit.>, as reported in the appendix for completeness. 




 §.§ Learning algorithm, tuning aspects and extensions


In the closing part of this section we report the overall identification algorithm, to be repeated for each prediction step p=1…p:

	
  * Collect N_p measured regressor values and the corresponding measured output instances (see (<ref>))

	
  * Solve the optimization problem (<ref>) and compute ϵ̂̅̂_p=αλ_p, that is needed to define the FPS (<ref>)

	
  * Solve (<ref>) to derive the nominal model θ^*_p and select a value of γ>1 to compute the worst case prediction error estimate τ̂_p(θ_p^*) (<ref>).

The main tuning parameters in the approach are the model order o, the disturbance bound d, and the scalars α and γ. In practical applications, a good estimate or even exact knowledge of o and d is often available from considerations on the plant and the available sensors. If not, an order selection procedure consists, as anticipated in section <ref>, in monitoring the estimate ϵ̂̅̂_p with increasing values of o: typically there is a clear convergence to a constant or slowly decreasing bound when o≥ n. We provide an example of this procedure in the next section. Regarding the choice of d, over- or under-estimating this value leads,  respectively, to either a too optimistic estimate of the prediction errors (since part of the model uncertainty is then hidden in the disturbance bound) or to a higher conservativeness (since the prediction error bound is then accounting also in part for the measurement noise). In a similar way, too large values of α and γ can lead to conservative error bounds, while values too close to 1 might give error bounds that are too tight and could be violated by new, previously unseen data. As a matter of fact, both conditions (over- and under-estimation) can be easily monitored on-line, and the derived tuning parameters d, α and γ suitably adjusted. This feature can be exploited in view of employing the proposed approach in an adaptive framework, which is a planned extension of this work. Other planned extensions that can be developed relatively easily are to consider nonlinear systems, since we can still employ the same model parametrization and the described approach by embedding the effect of nonlinearities in the error term ϵ_p(θ_p,φ_p) (see (<ref>)), and the use of other choices of linearly-parametrized models, e.g. truncated series of functional bases, since the same approach and main results would still hold, referred to the specific chosen class of models.




§ SIMULATION RESULTS


The proposed algorithm has been tested on the benchmark example already considered in <cit.>. In particular, the data are generated according to the continuous-time transfer function:

    G(s)=2229/(s+1)(s^2+30s+229)


The input-output samples are collected with a sampling time T_s = 0,2 s according to the output equation y(k) = z(kT_s)+d(kT_s), where d(t) is a colored noise obtained by low-pass filtering a randomly generated number, with a first-order filter with time constant of 0.2s. The disturbance d(t) is bounded in the interval [-0.2 0.2], ∀ t, and the prediction horizon is p̅=10.

The collected dataset is made of 500 input-output data samples. The input is a three level signal taking values randomly each 4s in the set {-1,0,1}.

In order to select the order of the model, we solve problem (<ref>) for different orders o, as reported in Fig. <ref>. The trend of the solution is monotonically decreasing, and this can be explained thanks to the noisy regressors that increase in number in vector φ_p as o grows (see (<ref>)). The effect of noise contained in them vanishes over time due to compensation of terms, eventually reducing the bound λ_p ∀ p=…p̅ for growing order o of the model. The final choice has been o=3, consistently with Assumption <ref>, that matches the order of the system, even though its true dynamics are dominated by the pole in s̅=-1. 



Once λ_p is computed ∀ p=1 …p̅, it is inflated according to (<ref>) with the fixed coefficient α=1.2, i.e. with a 20% additive margin. 
The conservativeness of such a choice can be evaluated from the trend that λ_p exhibits with respect to the number of data. In view of Theorem <ref>, in fact, lim_N_p→∞(ϵ̅_p^0-λ_p)=0^+ and the convergence rate gives a qualitative idea of how reliable ϵ̂̅̂_p^0 is with respect to the real unknown value ϵ_p^0. Fig. <ref> presents the values of λ_p for the fixed steps p=3,6,9 with an increasing percentage of the dataset used for the computation. Using the whole dataset, the curves approach an asymptotic value which presumably correspond to the theoretical bound ϵ̅_p^0.



The multi-step predictors are eventually built by solving (<ref>) to compute the nominal model for each step and the related guaranteed worst-case error bound τ̂_p(θ^*_p). The resulting bounds are dramatically less conservative then those coming from the iteration of the 1-step-ahead predictor, as shown in Fig. <ref>. 
We also carried out a comparison with the worst-case bounds associated to Least Square (LS) models, one for each p, reported in Fig. <ref> too. The result  confirms the optimality of our proposed approach, that minimizes the uncertainty bounds. 
Note also that it may happen that a model computed with LS does not belong to the FPS. 
Finally, the descending trend of τ̂_p(θ^*_p) is due to the increasing number of noise-free input values u(k) in the regressor as p increases, combined with the fading effect of past input and output data.



The final outcome of the algorithm consists in p̅ nominal models and as many guaranteed worst-case errors bounds, that are then tested with validation data. A detail of the simulation is reported in Fig. <ref>.






§ CONCLUSIONS

We presented an approach to derive multi-step prediction models, and the related uncertainty bounds, for an unknown linear system affected by additive measurement disturbance. Under suitable assumptions, we demonstrated convergence of the bounds to their theoretical minimum based on the available information. The approach requires the solution to linear programs only. The derived models are particularly suited to robust model predictive control design, since they can predict the future trajectory of the system on a finite horizon and the related uncertainty intervals.





§ APPENDIX





 §.§ Proof of Theorem <ref>


Problem (<ref>) can be rewritten as:

    [                                                                                 min_θ_p ∈Ω, ϵ_p ∈ℝϵ_p;                                                                                            subject to;                                                    |y_p-φ_p^Tθ_p|≤ϵ+d, ∀ (φ_p,y_p):[[ φ_p; y_p ]]∈𝒱_p ]
 

At first, note that the solution to (<ref>) must imply a regressor, denoted with φ_p,0, and a corresponding output value, y_p,0, satisfying the constraint in (<ref>) with equality:

    ϵ̅^o_p=| y_p,0-φ_p,0^Tθ_p^0 | - d̅ =min_θ_p∈Ω| y_p,0-φ_p,0^Tθ_p  | - d̅


Proof of claim 1). From the definition of 𝒱_p and 𝒱̃^N_p_p it holds that 𝒱̃^N_p_p ⊂𝒱_p; thus, from (<ref>) neglecting the trivial case λ_p=0, we have:

    λ_p(N_p) =    min_θ_p ∈Ωmax_[[ φ̃_p;  ỹ_p ]]∈𝒱̃_p^N_p|ỹ_p - φ̃_p^Tθ_p | - d̅≤
       min_θ_p ∈Ωmax_[[ φ_p; y_p ]]∈𝒱_p|ỹ_p - φ̃_p^Tθ_p | - d̅


In view of (<ref>) we have then

    λ_p(N_p)  ≤min_θ_p ∈Ω| y_p,0- φ_p,0^Tθ_p  | - d̅

which implies λ_p(N_p)  ≤ϵ̅_p^0.



Proof of claim 2). Starting from (<ref>), with standard properties of absolute values, we compute

    λ_p(N_p)=   max{0,min_θ_p ∈Ωmax_[[ φ̃_p;  ỹ_p ]]∈𝒱̃_p^N_p|ỹ_p-φ̃_p^Tθ_p| - d̅}≥
       min_θ_p ∈Ω{|y̅_p (N_p)-φ̅_p(N_p)^T θ_p| - d̅}


where 

    [ φ̅_p(N_p); y̅_p(N_p) ]= min_[[ φ̃_p;  ỹ_p ]]∈𝒱̃_p^N_p[[ φ_p,0; y_p,0 ]]-[[ φ̃_p;  ỹ_p ]] _2



Now adding and subtracting y_p,0 and φ_p,0^Tθ_p from (<ref>) and neglecting the trivial case λ_p=0:

    λ_p(N_p) ≥   min_θ_p ∈Ω{| (y̅_p(N_p)-y_p,0) + 
       (φ_p,0-φ̅_p(N_p))^Tθ_p + y_p,0 - φ_p,0^Tθ_p | - d̅}
    ≥   min_θ_p ∈Ω{|  y_p,0 - φ_p,0^Tθ_p - (-y̅_p(N_p)+y_p,0) + 
       (φ_p,0-φ̅_p(N_p))^Tθ_p | - d̅}
    ≥   min_θ_p ∈Ω{|  y_p,0 - φ_p,0^Tθ_p |}-
       min_θ_p ∈Ω{| (-y̅_p(N_p) + y_p,0) +
        (-φ_p,0+φ̅_p(N_p))^Tθ_p | - d̅}
    
     =   ϵ̅_p^0 + d̅-
    min_θ_p ∈Ω{| (-y̅_p(N_p) + y_p,0) +
        (-φ_p,0+φ̅_p(N_p))^Tθ_p | - d̅}


Then, simplifying d̅, we compute

    λ_p(N_p) ≥ϵ̅_p^0 - min_θ_p ∈Ω   {| (-y̅_p(N_p) + y_p,0) 
       + (-φ_p,0+φ̅_p(N_p))^Tθ_p |}



Considering Assumption <ref>, we know that for N_p →∞

    ∀β>0, ∃N̅_p(β): φ̅_p(N̅_p)-φ_p,0≤β, |y̅_p(N̅_p)-y_p,0|≤β


and thus, using the inequality | a^Tb|≤a_2b_2

    [                                      λ_p(N̅_p(β)) ≥ϵ̅_p^0 - min_θ_p ∈Ω{| (-y̅_p(N̅_p(β)) + y_p,0) |;                                                                   +θ_p _2  φ_p,0-φ̅_p(N̅_p(β)) _2 };                                                                ≥ϵ̅_p^0 - β(1+ sup_θ_p ∈Θ_p^0θ_p_2 ) ]


Finally, choosing β≤ρ/(1+ sup_θ_p ∈Θ_p^0θ_p_2) concludes the proof.





 §.§ Reformulation of problem (<ref>)

The problem (<ref>) can be stated as  
    min_θ_p∈Θ_pmax_i=1,…,N_pmax_θ∈Θ_p|φ̃_p(i)^T(θ-θ_p) |


The absolute value in (<ref>) can be separated into two terms:

    min_θ_p∈Θ_pmax_θ∈Θ_pmax_i=1… N_pmax[  φ̃_p(i)^T(θ-θ_p); -φ̃_p(i)^T(θ-θ_p) ]


Let us now introduce the quantity φ̌_p(j)=φ_p(i)     if     i≤ N_p

-φ_p(i)    if    i > N_p and recast (<ref>) accordingly:

    min_θ_p∈Θ_pmax_θ∈Θ_pmax_j=1… 2N_pφ̌_p(j)^T(θ-θ_p) = 
       min_θ_p∈Θ_pmax_j=1… 2N_pmax_θ∈Θ_pφ̌_p(j)^T(θ-θ_p) = 
       min_θ_p∈Θ_pmax_j=1… 2N_pmax_θ∈Θ_p(φ̌_p(j)^Tθ -φ̌_p(j)^Tθ_p)

Let us define:

    c_j=max_θ∈Θ_pφ̌_p(j)^Tθ, j=1,…,2N_p

Then, (<ref>) is equivalent to:

    min_θ_p∈Θ_pmax_j=1… 2N_p (c_j -φ̌_p(j)^Tθ_p)


Problem (<ref>) corresponds to a single LP that reads:

    min_θ_p∈Θ_pζ
    
    s.t.      c_j-φ̌_p(j)θ_p ≤ζ,    j=1… 2N_p
 


plain  

1
	
	boyd2004convex
	S. Boyd and L. Vandenberghe.
	Convex optimization.
	Cambridge university press, 2004.
	
	Ljung99
	L. Ljung.
	System Identification - Theory For the User (2nd ed.).
	PTR Prentice Hall, Upper Saddle River, N.J., 1999.
	
	MNPW96
	M. Milanese, J. Norton, H. Piet-Lahanier, and E. Walter (Eds.).
	Bounding Approaches to System Identification.
	Springer US, Boston, MA, 1996.
	
	shook1991identification
	D. Shook, M. Coorous, and S.L. Shah.
	Identification for long-range predictive control.
	In IEE Proceedings D (Control Theory and Applications), volume
	138, pages 75–84. IET, 1991.
	
	TFFS18b
	E. Terzi, M. Farina, L. Fagiano, and R. Scattolini.
	Robust predictive control with data-based multi-step prediction
	models.
	In European Control Conference 2018.
	
	Traub80
	J.F. Traub and H. Woźniakowski.
	A General Theory of Optimal Algorithms.
	Academic Press, New York, 1980.
	



