On-line parameter and state estimation of an air handling unit model: experimental results using the modulating function method

This paper considers the on-line implementation of the modulating function method, for parameter and state estimation, for the model of an air-handling unit, the central element of HVAC systems. After recalling the few elements of the method, more attention is paid on issues related to its on-line implementation, issues for which we use two different techniques. Experimental results are obtained after implementation of the algorithms on a heat flow experiment, and they are compared with conventional techniques (conventional tools from Matlab for parameter estimation, and a simple Luenberger observer for state estimation) for their validation.


Introduction
The challenges of the clear evidence of climate change [13] bring in the front-line new regulations which aim to reduce the green-house gas emissions in all sectors. A notable amount of emissions as well as energy consumption is accounted by buildings, both in the industrial and residential sector. According to [34] the need for heating is going to decrease while the cooling demand in buildings is going to increase in the coming decades due to the climate change. In this respect, heating ventilation and air conditioning (HVAC) systems with higher energy efficiency and better building designs are continuously researched and developed.
To facilitate the transition towards reducing the energy consumption and the green-house gas emissions, mathematical models of buildings have been intensively used. The simulation tools available today [10] and the ongoing research offer a high variety of possibilities when it comes to design optimization [21], renewable integration [41], retrofitting [14], on-line fault detection diagnosis [7] and even optimized real-time control.
The HVAC system is one very important part of the building, which in Europe is estimated to share 76% of the total energy use [17]. The duty of this system is to ensure a comfortable indoor environment by regulating the temperature and the air quality. The well-functioning of the HVAC system will contribute also to an optimal energy consumption. Several subsystems are considered when modeling an HVAC system, as explained in [36]. Among them, the air handling unit (AHU) is the subsystem whose role is to condition the air circulated through the building. Different models have been proposed in the literature for modeling an AHU: physical models described by partial differential equations used by most of the simulation tools available [10], transfer functions [4], data-driven models [2], and hybrid models or reduced-order models (ROM) [3].
The latter lead to many possibilities as they easily link to estimation and control scenarios where real-time capabilities are important. In this paper, we investigate the use of the so-called modulating function method to perform both This paper addresses the question of real-time deterministic parameter and state estimation of an air handling unit using the modulating function method. The corresponding algorithms are implemented in Matlab/Simulink on an actual heat flow model and experimental results based on noisy measurements are presented. While the model itself is easy to develop and implement, performing estimation of its parameters and states is still a challenge given the actual distributed nature of these parameters, which are also subject to possible changes over time. In doing so, we also look at practical aspects related to implementation issues for using the modulating function method in an on-line situation, and which we believe has been less considered in the literature. Preliminary results were reported in [18].
After this introduction, we start this paper by briefly explaining the principle of an HVAC system and that of an AHU, and give a simple reduced-order model (ROM) of the latter (section 2). Then, the basics of the modulating function method are recalled for both parameter and the state estimation in a state-space context (section 3). A following section is more specifically addressed to issues that are more prone to arise in an on-line scenario, i.e. possible singularities coming from the choice of the modulating function, and ease of implementation of the moving-horizon version in a block diagram environment such as Simulink. For this, we use two different techniques, each one dedicated to these two specific issues (section 4). We validate the proposed techniques on a heat flow model and discuss our experimental results (section 5). Brief concluding remarks end this paper (section 6). The HVAC system within a building can be configured in many different ways with regards to the buildings' layout and purpose, application or users' needs. It usually contains a variety of components such as chillers, cooling towers, condensing boilers, AHUs, heat pumps, heat recovery units, fans, control valves, pipes, etc. as described, for example, in [36]. The air flow in a typical compact HVAC system can be represented by the schematic diagram shown in Fig.  1. The heating and cooling coils are part of separate networks with complex units that deliver cold or heat based on requirements. Similar representations can be seen in [12] or [8].
An important component of the HVAC system is the AHU (marked with green dashed box in Fig. 1). The role of this component is to condition the air that is used for the ventilation of the rooms in order to maintain a comfortable indoor environment in terms of temperature and air quality. In this unit, the fresh air is mixed with air coming from the rooms in a mixing chamber, then is circulated by a constant air volume flow fan over a cooling/heating coil. Based on the needs, the air is heated or cooled accordingly, so that the temperature at the exit of the unit achieves prescribed values.

Mathematical modeling
Many models have been developed to represent the dynamic heat flow behavior of an AHU. In the existing literature, the models that have been developed range from the use basic physics [37] to advanced data mining algorithms such as artificial neural networks and other machine learning techniques [1]. An overview of the variety of these models can be found in [4]. The heat flow which is considered here is schematically represented in Fig. 2. It simply consist of a fan, a heating coil, and a duct where temperature sensors can be placed. As expected from basic considerations, the temperature profile of the flowing medium, as well as the heat transfer coefficients typically varies with the length L of the duct (see for example [6]). However, in the present application, only the temperature at the exit of the unit is of interest as it is typically the one which one wants to regulate. Hence we will, in the following, consider the measured temperature T m,L , hereafter simply shortened as T m . In order to model this heat flow, we resort to second-order linear ROM. A similar albeit first-order model is described in [3]. The proposed 2-nd order linear model uses the resistance capacitor (RC) network analogy depicted in Fig. 3. In order to develop the model, we first write a flow balance equation for each node. Hence, we get for node 1, representing the air inside the unit, where T m is the measured temperature in the AHU and T e is the envelope/surface temperature of the AHU, while Q h is the heating/cooling load added in the AHU. The constant parameters of this flow balance equation are C m , the summed heat capacity of the sensor and of the air inside the unit and R ms , the thermal resistance between the temperature sensor and the envelope of the unit. For node 2, the surface/envelope of the unit, we have the differential relation where T r is the temperature of the room surrounding the unit. In equation (2), constant parameters are C s , the heat capacity of the envelope of the AHU, R ms , the thermal resistance between the sensor and the surface of the unit, and R sr , the total thermal resistance between the surface of the unit and the room.
From (1)-(2), it is simple to obtain the state-space representatioṅ where, defining the state vector and the input vector as and while assuming we only measure the temperature at the end of the unit, i.e. y = T m , we have the matrices 3 The Modulating Function Method: parameter and state estimation In this section, we briefly recall the basics of the modulating function method. For the sake of clarity, we do so on single-input single-output systems, with the multiple-input case as a direct extension.
The modulating function method being primarily based on an input-output description of a dynamical system, we start, as a preliminary by transforming our state-space presentation (3)-(4). Hence, considering the n-dimensional case where x ∈ R n , u ∈ R and y ∈ R, we differentiate the output equation (4) n − 1 times and get the well-known expression whereȳ = [y,ẏ, ..., y (n−1) ] T andū = [u,u, ..., u (n−1) ] T . Matrix O is the well-known observability matrix of the Kalman criterion and matrix T is the system specific Toeplitz matrix given by Then, differentiating equation (4) one more time and isolating x in expression (10), we get where the reversed controllability matrix C R is given by Note that an obvious condition for expression (12) to be defined is the observability of state-space representation (3)-(4). We thus have the input-output form where a T = CA n O −1 = [a 0 , a 1 , ..., a n−1 ] Defining then the vector Y as Y = [ȳ T ,ū T ] T , the following parametric form is obtained where the unknown parameter vector θ ∈ R 2n is given by θ = [−a T , b T ] T . In case a constant and unknown disturbance d is impacting the system's behavior at the same level as the input, expression (15) can be simply modified Before proceeding, let us first recall the definition of a modulating function (see [23], [19]). Definition 1. The sufficiently smooth function ϕ : [0, T ] → R is called a modulating function (of order k) if at least one of its boundaries and its derivatives up to order k equals zero, i.e. if A modulating function where ϕ (i) (0) = 0 and ϕ (i) (T ) = 0 (for i = 0, k − 1) is called a left modulating function, while a modulating function with ϕ (i) (0) = 0 and ϕ (i) (T ) = 0 is called a right modulating function. If a modulating function is such that we have ϕ (i) (0) = ϕ (i) (T ) = 0, then it is called a total modulating function.
Note a more general definition of a modulating function (see [19]) allows to include expanding horizons, allowing thus to include alternative integral approaches [28] [35].
In this paper, we are interested in on-line estimation. Regarding parameter estimation, we hereby use the following receding-horizon integral operator on the signal y(t), using a total modulating function given by [30] L where i = 0, k − 1 and T > 0. Due to the fact that ϕ(t) is a total modulating function and by simple integration by parts, operator (17) has the important property that Then, applying operator L 0 [·] on each time-varying signal of (15), and using property (18), we are able to avoid both explicit time-derivatives of measured signals and unknown initial conditions to get where Finally, we can obtain an estimate of parameter vector θ by either proceeding to a conventional receding-horizon Gramian-based estimator [32] [33] over a horizon T > 0, or aggregate a sufficient number of expressions (19), each one obtained with a different total modulating function, so that with m t ≥ 2n total modulating functions, we get the set of linear equations In this case, we can in principle directly get the parameter vector estimate by computingθ = WW T −1 Wz (23) or, alternatively, proceed by adding another receding-horizon stage similar to the Gramian-based estimator alluded to above in order to remove noise further, and define the estimate aŝ Turning now to state estimation, we resort to an integral operator based, this time, on a left modulating function, and given by where ϕ l (t) ∈ [0, T ] is a left modulating function. Because of the fact that ϕ l (T ) = 0, we have Proceeding then by applying operator L 0 l [.] on each signal of (14), we get the expression where and (and similarly for L l [u]).
Noticing then that expression (27) can be put into a form similar to (19), i.e. we have where and Combining then m l ≥ n left modulating functions, we obtain, similarly to (22), the expression which, defining the state corresponding to an observability canonical form as leads to its estimatex with an expression similar to (23). Alternatively, we can also estimate the original state of the system (3)-(4) by either using (10) or the simple transformation 4 Implementing the modulating function method for on-line applications 4

.1 Using time-varying modulating functions
Early works on the use of the modulating function method for offline identification of unknown parameters (see for example [31]) start with the definition of a set of specific and predefined modulating functions, which can take different forms (ie Hartly, splines, trigonometric functions, etc). Once a sufficient number of modulating functions are used, one can get a parameter estimate directly by simple inversion, as in (22). However, when considering an on-line and receding-horizon situation, this can generally lead, especially with the presence of noise on the measurements, to singularity issues. For example, considering the trivial system applying the modulating function method would simply lead tô where the denominator of (39) could create difficulties. As an illustration, we have simulated the sinusoidal signal y(t) = 15 sin 2t, corrupted with a uniform noise. As can be seen in Fig. 4, the estimate of a 0 , obtained with the fixed total modulating function ϕ(t) = t 2 (t − T ) 2 , shows singularities around the time when the signal L 0 [y](t) has its zero-crossings (see also similar spiking behaviors in Figure 3 of [9]). Seeing L 0 [y](t) =< ϕ, y > as a dot product between two functions, one has therefore an orthogonality issue between these functions.
While one way to go around this issue typically consists of using a gradient or recursive least-squares algorithm (or even using more modulating functions), another possibility is to continuously redefine the modulating functions based on the current batch of the measured signals y(t) and u(t). Indeed, setting the number of total modulating functions to m t = 2n, we impose a normalizing constraint represented as W = I mt .
(40) In this way, we are simply bypassing the matrix inversion of (22) or (24). Realizing this constraint hence consists of finding the set of m t modulating functions that will fulfil (40). Since W is composed of m t vectors (21), where L i [y] is given by (17), we have a set of m t integro-differential equations, where the unknowns are the m t total modulating functions ϕ k 's, which are also constrained at both their boundaries due to their total modulating function nature.
Let us then transform this set of integro-differential equations into integral equations by considering the n-th order derivative of ϕ(t) arising in expression (20), and define the new function α : [0, T ] → R as Then, using the anti-derivative notation the right boundaries of the derivatives of a total modulating function ϕ(t) can simply be re-written as while the zero left boundary conditions are simply fulfilled thanks to the successive integration of α. Thus, integral operator can now be re-written as so that w in (21) only consists of integrals of α. Hence, expression (40) is now a set of integral equations. Grouping the right boundary conditions (43) similarly, we also get where Γ = [α 1 , α 2 , ..., α mt ], with each vector α j composed by all boundary conditions (43) for this specific modulating function, i.e. α j = [α Combining finally (40) with (45), we have a system of linear integral equations, which is simple to solve numerically (details of this numerical method are given in Appendix A). Once a function α(·) is obtained for each instant t, constraint (40) is fulfilled, and Gramian-like expression (24) is simply replaced witĥ As for state estimation, it is also possible to normalize W l similarly so that we have the state estimatê The most notable difference with parameter estimation, is that, since W l in (35) does not depend directly on the measured signals, the m l = n left modulating functions ϕ l,j (·) can be chosen once and for all t. Note that a similar kind of normalization was also mentioned in the context of least-squares observers (see [29]).

Direct continuous-time implementation
When needing to choose a particular estimation method, one can obviously focus on performance in terms of the best possible match between the output of the considered plant and its corresponding predicted output using the parameter estimates. However, other factors can also be under consideration, such as ease of implementation, portability of the code, etc. Regarding the on-line use of modulating functions, the latter are usually first discretized (as they would be in the previous subsection or as in [9] and [30]). However, modern tools for simulation and control such as Matlab/Simulink allows one to program systems and algorithms directly in a continuous-time setting, thus allowing for the additional advantage of being able to consider non-regular samplings. In this section, we present one way to do that, with in mind direct continuous-time implementation, as opposed to looking primarly at performance. To do so, we first begin by introducing a function ψ : [0, T ] → R, which we will refer to as reversed modulating function, and where ψ(·) is such that ψ(t) := ϕ(T − t) (48) where ϕ(·) is our usual modulating function of the very basic following definition 1. In this case, note that, because of reversal (48), a left reversed modulating function corresponds to a right modulating function, and vice-versa. Then, replace operator (17) with The advantage of (49), is that, besides its slightly simpler expression than (17), it is directly put under a usual convolution form. Note, also similarly to (18), we have the property that M 0 [y (i) ] = M i [y].
Next, we notice that many modulating functions defined in the literature can be expressed by the solution of a differential equation. For example, the total modulating function ϕ(t) = t 2 (t − T ) 2 used for example (38) is the solution of differential equation ϕ (4) (t) = 0. Hence, we introduce the following state-space representatioṅ where the state vector χ ∈ R n ψ .Thus, integral operator (49) can easily be rewritten using some advantages of state-space representations. For example, M 0 [y] can be rewritten as and similarly for the other M i [y]'s. Defining now the new vector ξ y ∈ R n ψ as then it can be shown that where vector ξ y (t) is obtained as the solution of differential equatioṅ ξ y = Λξ y + ly.
Interestingly, assuming now that system (50) is stable means that filter (55)-(54) can be used to implement M 0 [y] in software tools such as Matlab/Simulink without having to proceed to a preliminary discretization (the same is of course valid for M 0 [u], with a state vector ξ u ). Note that computing M i [y] is not more difficult, and we would simply have Gathering now the M i [y] terms similarly to (30), we have (61) In case one wants to favor speed over precision, it is possible to use a single total (reversed) modulating function and use a Gramian-based expression in order to get the parameter estimate. Another advantage of using expressions such as (52) or (54), where the modulating function is generated by a stable system, is that similar expressions can also be used to obtain the parameter estimates themselves. Indeed, as in [32], we can use a generalized Gramian expression where the kernel is generated thanks to a reversed modulating function. Indeed, multiply each term of (19) by ψ(t − τ )w(τ ) and integrate to get which can be rewritten as where ξ h ∈ R 2nn ψ , and where Λ = Λ ⊗ I 2n , l = l ⊗ I 2n and Σ = Σ ⊗ I 2n (with ⊗ for the Kronecker symbol). For M 0 [G], we have same kind of filter, this time with a matrix differential equatioṅ where Ξ G ∈ R 2nn ψ ×2n . Finally, an estimate of parameter vector θ is obtained by simple inversion aŝ Interestingly, and moving now to state estimation, the method simply consists in repeating the steps (32) to (37) (note, therefore, that steps (60) through (68) are obviously not necessary).

Case study
As a case study, a heat flow experimental chamber produced by Quanser company is considered. This plant reproduces the thermodynamics of an Air Handling Unit of an HVAC system. A few technical details related to the experimental set-up will be first introduced, while results of on-line parameter and state estimation and their validation will be presented afterwards.

Quanser heat flow chamber
The heat flow experiment (Fig. 5) consists of a fiber-glass chamber equipped with a fan blowing over an electric heating coil. Both the fan and the heater can be controlled externally with a 0 − 5V input signal. The air temperature inside the box is measured by three temperature sensors positioned equidistantly. Additionally, we make use of the Quanser Q8-USB Data Acquisition Board in order to enable the communication between the computer and the experimental chamber.

Initial set-up and simulation environment
In the literature, the considered system model for control applications is usually of first order model (with or without time delay) [27] [40]. However, in practice, a second-order model can have its importance, as the surface does not only represent the envelope/box storing heat but also other components (fans, dampers, filters). In a faulty situation where the box /components inside get overheated and can be damaged, estimating the box temperature can generate important information, and help in fault detection. The set-up will be used both for parameter estimation and state estimation scenarios where the latter uses the results of the former. Because the chamber is located in a closed indoor environment and the experiment is conducted on short time interval, we will assume that the room temperature (T r ) is unknown but constant during the experiment. In order to stay close to an actual operation of a constant flow AHU, the fan is operating with a constant speed during the whole experiment.

On-line parameter estimation
Using steps (10) to (15) on state-space representation (3)-(9), we obtain the ordinary differential equation where y(t) = T m (t) is the measured output, u(t) is the input Q h , and a 1 , a 0 , b 1 , b 0 and d are the unknown coefficients of the equation, the last one, d representing the impact of the unknown input T r mentioned above, and considered here as a disturbance. These coefficients are related to the parameters of system (1)-(2) as follows: To perform the parameter estimation we have used the integration period T = 2000 sec with a sampling time T s = 2 sec, thus giving a total of 1000 samples. The additional receding horizon interval for obtaining the estimates is T = 2000 sec.

Input profile and persistence of excitation
A key point in parameter estimation is the persistence of excitation of the input signal. Reliable estimates will be obtained if the input signal is sufficiently rich so that the observed response contains the required information to perform the estimation process. As specified in [24], if a specific matrix characteristic of the input signal is non-singular, the input is considered to be persistent. It is well-known that, when estimating the parameters of a system with an input signal having enough persistence of excitation, the estimated parameters approach their true values [20]. However, it is worth mentioning that, unfortunately, not enough attention has been given to the consideration of persistent inputs in building models and thus not many studies can be found on this topic [22].  The considered input signal for parameter estimation is a pulse function shown in Fig. 6. This signal is informative enough to capture the dynamics of the heating chamber and obtain the estimates. The amplitude of signal is selected to 1.5V to change the temperature of chamber up to around 30 •C . The period of the signal is 1200 sec which is long enough for the plant to reach at least 80% of its final value.

Estimates obtained using the modulating function
The results of parameter estimation using the modulating function method described in section 3 are presented in Fig.  7. The results are given for both methods: the time-vayring modulating function technique (continuous lines) and the direct continuous-time technique (dashed lines). As expected, the direct continuous-time technique, using only one modulating function, does not perform as well as the time-varying technique (although it is computationally more efficient). The latter shows quite good results and converges to an almost constant value after the additional receding horizon interval T .
Relating to the time-vayring modulating function method presented in section 4.1, Fig. 8 shows how one of the total modulating functions (corresponding in this case to a 1 ) evolves over time. The modulating function is a signal of length T = 2000sec, and the figure shows its evolution from t = 2000 to t = 3700 (the curve highlighted in red is MF at time obtained t = 2000sec). As it can be seen, the amplitude of the MF is constantly adapting to new incoming data/measurements that are fed on-line to the estimator.

Validation
Since there is no upfront knowledge about the real values of the parameters, we will compare the results of the modulating function method with the results of two well-known estimation methods available in the System Identification Toolbox for Matlab [25]. These two methods, applicable to continuous-time systems, are the Transfer Function method and the Grey-box estimation method. Moreover, for further validation of the results, the system output is reconstructed using the estimated parameters compared with the real measurements.

Continuous Transfer Function Estimation
The first method used for comparison is the "tfest" function of Matlab which estimates a continuous transfer function for the given data. In this method the input and output of the model are filtered using a pre-filter L(s), and the differential equation of the system is written in a regression form with respect to the filtered signals. The parameters are then estimated to minimize the prediction error.
Matlab uses different numerical search methods to estimate the parameters iteratively. In case of selecting the default setting, a combination of four line search algorithms 'Subspace Gauss-Newton least squares search', 'Levenberg-Marquardt least squares search', 'Adaptive subspace Gauss-Newton search', and 'Steepest descent least squares search' methods will be applied and the first descent direction leading to a reduction in estimation cost is used.
In order to initialize the parameters, different algorithms are available in Matlab. The algorithm used in this study is the instrument variable approach in which the parameters are estimated using the least-squares method [16] and [26].

Grey-Box Model Estimation
The second method used for comparison is the grey-box model estimation using the "idgrey" function of Matlab. This method directly uses the parametric state-space representation of the system given in (3)-(4) and minimizes the prediction error. Matlab uses the same numerical search methods as explained in section 5.4.1. However, contrary to the transfer function estimation technique, the initial value for the parameters should be given by the user.  Table 1 shows the estimated parameters using the MF method together with those obtained by the transfer function and grey-box approaches. As it can be seen, the results of the MF method are very close to the results of the transfer function method, while there is a discrepancy in some parameters with when compared to the grey-box method. The comparison is more clear in Fig. 9 where the output of the estimated models are plotted together and compared with the measured data. The goodness of the fit is 96.95%, 96.51% and 97.42%, for the modulating function, transfer function and gray-box methods, respectively, which confirms the validity of the proposed MF method.
The main reason for achieving a slightly different (better) result by grey-box estimation method is that the reported goodness of fit was based on non-filtered measured data which matches the structure of grey-box method. The grey-box method works on the original input-output signals whereas the other two methods, i.e. modulating function and transfer function methods, use a pre-filter for the input-output signals (L(s) in transfer function method and the modulating integral operators in the MF method). Therefore the optimization algorithm in the MF and transfer function methods minimizes the filtered error while the reported goodness of fit was based on the original error. In case we would define the goodness of fit based on filtered signals, a higher fit would be achieved for MF and transfer function methods.
Another important point that should be mentioned is that our implementation of the MF method is on-line while the other two methods are iterative and off-line. In particular, the convergence of the grey-box method highly depends on the initial guess for the parameters (in this study, the estimated parameters from MF method are considered as initial guess for the grey-box method). Therefore, achieving a result in the order of off-line methods in finite time further attests the acceptable performance of our results.

On-line state estimation
Given the estimated parameters from previous section, it is straightforward to apply the method described in section 3 to obtain the state estimates. Here, both T m and T s are estimated on-line using an integration interval T = 50sec with sampling period T s = 1sec.
In a normal operation set-up, the requirements for the temperature at the exit of the AHU might not vary too much in a short time interval. However, to test the proposed state estimation algorithm in different work conditions, a pseudo random input profile shown in Fig. 10 is applied to the system. This is not a common operating profile, but it helps to visualize the state estimator and evaluate its behavior.   In Fig. 11, it can be observed that the estimated temperature inside the chamber is very close to the noisy measured temperature. Also, as expected, the value for the state representing the envelope temperature of the chamber is lower than the temperature inside, and follows the same heating profile which is in perfect concordance with the heating input profile. The estimated states are comparably good with the ones obtained from a standard observer. In this specific case, the convergence of the states is faster due to the fixed integration interval, which ensures convergence after the first receding horizon window.

Conclusion
A 2nd order model is presented in this paper to represent the heat flow in the air handling unit present in compact HVAC systems. To facilitate the real-time applications of the model we have proposed an unconventional method, the modulating function, with two different implementation approaches.  Figure 11: Measured vs. estimated temperature inside the heating chamber using the modulating function and simple observer.
As a case study, a heat flow experiment by Quanser, whose thermodynamics are very close to those of an AHU, is considered. The on-line parameter and state estimation algorithms were implemented in Matlab/Simulink. The results underline the potential of this approach, showing a good match between the measured and estimated parameters and states. A good convergence of the parameters was observed having as a baseline the parameters estimated with standard estimation methods: the transfer function method and the grey-box method. Moreover, using several (as opposed to only one) modulating functions gives a better estimation of the parameters compared with an approach that uses only one, at the cost of speed although the latter one mentioned requires less resources to implement and run. The method can be successfully used for on-line state estimation as well, the results being comparable good with a standard observer.
and for i = 3, n: The operators L i [y] of w in (44) can be similarly approximated, so that  (40) can thus be approximated by where α = [α 1 , α 2 , ..., α mt ] and the matrix K is given by Proceeding similarly with the discrete approximation of Γ in (45), it can be expressed where B =      T n s 1 T Q n−1 . . .
The discrete approximation of z is given by z ≈ z = T s y T α.
The parameter estimate vector is finally obtained after applying the simple discretized version of receding-horizon expression (46).