Handling of time delays in system identification with regularized FIR models

3.1 Regularized FIR models

Due to its model structure, a large number of parameters is required to cover the dynamics of the given process properly. By regularization, the effective number of parameters can drastically be reduced. Therefore, a satisfying small variance error can be achieved. Thus, the major drawback of FIR models can be overcome. By introducing an additional penalty term in the estimation procedure, regularization is applied:

where R ̲ is the regularization matrix which is dependent on the hyperparameters η ̲ . The hyperparameter λ denotes the regularization strength. The second term is the regularization term which introduces a relationship between the parameters of the FIR model themselves. The linking of adjacent parameters is defined by the matrix R ̲ . Due to the fact that the parameters of an FIR model correspond to the impulse response coefficients of the modeled process, prior knowledge on the dynamic behavior can be incorporated in the matrix R ̲ . From the Bayesian perspective the model parameters θ ̲ can be interpreted as a Gaussian process. The parameter vector θ ̲ is assumed to be a random variable with a prior distribution of zero mean and the covariance matrix P ̲ [10], which is also known as the kernel matrix. It has to be noted that R ̲ corresponds to the inverse parameter covariance matrix P ̲ − 1 of the prior. The hyperparameter λ is used to adjust the strength of the regularization. With a regularization strength of zero the prior knowledge in R ̲ has no influence on the estimated parameters. In contrast, for a high regularization strength almost only the prior is used for the parameter estimation. Furthermore, this approach preserves the analytical form of the solution from LS [2, 3]:

3.2 Identification of time-delay systems with regularized FIR models

A major benefit of FIR models in general is the property that they are robust w.r.t. unknown time delays, because each parameter can be independently set to zero. By regularization neighboring parameters are coupled and this advantage is weakened or lost. This coupling leads to undesired behavior, because only after the time delay, the model parameters should deviate from zero. Thereby, regularization causes a blurring of the parameter vector, which normally changes suddenly after the time delay has passed. This behavior can be seen in Figure 2. It causes that the regularization strength can no longer be chosen high, because the regularization penalty term does not fit to the process anymore. To overcome this issue, the regularization matrix (11) has to be modified to also allow the incorporation of time delays.

Figure 2:

Model parameters θ _j from the estimation of a second-order process with time delay with an unregularized () and a regularized FIR model ().

4.1 Time-delay kernel

From the Bayesian perspective the inverse of the regularization matrix P ̲ = R ̲ − 1 can be interpreted as the covariance matrix of a Gaussian prior with zero mean. Partitioning the prior of the parameters into two parts (‘d’ and ‘dyn’) leads to

The zero off-diagonal elements account for the fact that parameters corresponding to the time delay and parameters corresponding to the dynamic behavior are not correlated with each other. Additionally, two different covariance matrices P ̲ d and P ̲ dyn represent the two parameter groups.

For the second part ( P ̲ dyn ) it can clearly be argued that all assumptions on the undelayed system, such as exponential decay and smoothness in the TC kernel case or impulse response preservation for the IRP kernel, still hold. For the first part ( P ̲ d ) , however, these assumptions do not hold. So the questions arises, how to build the covariance matrix P ̲ d for the time-delay parameters.

A natural choice can be made from the fact, that all time-delay parameters should be exactly zero. Due to the absolute certainty of this condition, all variances of the parameter prior, which populate the diagonal elements of the covariance matrix, can be set to zero. Additionally, no covariance between the parameters exist which results in zero entries at the off-diagonal elements in P ̲ d as well.

The goal of the identification from the Bayesian perspective is to find the parameter distribution conditioned on the measured outputs. The mean of this distribution corresponds to the regularized LS solution [3]

By splitting the parameter vector

as well as the regressor matrix

into a time-delay and a dynamic part and including this into the mean of the parameter distribution we get

With the previous choice of P ̲ d = 0 ̲ , from (18) it can be reasoned that θ ̲ d = 0 ̲ since P ̲ d X ̲ d T = 0 ̲ . Additionally, due to X ̲ d P ̲ d X ̲ d T = 0 ̲ the regressor part X ̲ d is not influencing the estimation of the dynamic parameters θ ̲ dyn . Which is intended by the choice of P ̲ d .

From the implementation perspective as well as the filter perspective however, a choice of P ̲ d = 0 ̲ is questionable. For the filter interpretation the inverse of the covariance matrix

is needed. This inversion is not possible with a zeros matrix for P ̲ d since it leads to a rank-deficient matrix P ̲ . Also for small values of σ ² the inversion in (18) is poorly conditioned. This problem is also known from the IRP kernel and can be overcome executing a SVD, for details see [24].

A second possibility is exploiting an approximation, that preserves the desired property of forcing all parameters corresponding to the time delay to zero. It is given by

With β the variance of the prior of all time-delay parameters can be adjusted. Compared to the previously explained special case β = 0, increasing β results in higher variance of the parameter prior. Still, with small values for β similar results can be achieved as with β = 0. This choice coincidences with the well-known ridge regression but only for the parameters describing the time delay.

For the filter perspective of the kernel the inverse of P ̲ d is needed which can be calculated by

This results in an additional penalty term in the LS problem for θ ̲ ̂ = arg min θ ̲ J with

As explained in Section 3.1, we intend using the IRP kernel for the dynamic parameters. With the time-delay extension we propose the novel kernel named impulse response and time-delay preserving (IRDP) kernel which leads to the following regularization matrix

with R ̲ d ∈ R n d × n d and R ̲ IRDP ∈ R n × n .

While deriving the time-delay kernel, the assumption of exactly known time delay is made. This, however, in most practical situations does not hold. Therefore, a mechanism has to be developed identifying the time delay along with the FIR parameters, to properly build the time-delay kernel. This problem is addressed in the following sections and solved in two different ways.

4.2 Single kernel method

The first approach for identifying RFIR models with unknown time delays is based on a single kernel method. As explained in Section 4.1 the incorporation of the time delay into the regularization matrix is done by strict distinction between the time-delay and the dynamics part. Thus, choosing the point of separation n _d in the regularization matrix is an integer optimization problem. Solving it along with the other real hyperparameters ( a ̲ , λ) results in a mixed-integer optimization problem which usually is not easy to solve.

Therefore, we modify the kernel by fuzzifying the separation. First, we build two separate regularization matrices of size n × n. The first purely consists of the time-delay kernel regularization matrix. The second is the regular dynamic kernel regularization matrix, which can be decomposed into the filter matrix as shown in (9). Then both kernels are superimposed in a weighted manner. A sigmoid function and its remainder to one are used for this purpose. The combined regularization matrix is given by

with R ̲ d ∈ R n × n . It has to be noted that the rows of F ̲ dyn are weighted with diag ( 1 ̲ − μ ̲ ) and therefore, the model parameters in the corresponding rows with zero weighting are not penalized. The values for μ ̲ are calculated by evaluating the sigmoid function

at different query points. Here y defines the inflection point of the sigmoid and z the slope. The query points are chosen corresponding to the number of FIR parameters:

Figure 3 visualizes this approach. With it, the mixed-integer optimization problem of the time-delay selection is changed into a continuous optimization problem w.r.t. the hyperparameter d ̂ . This may result in a slightly different solution, but it is much easier to solve. Additionally, the fuzzification strength can be adjusted by the second hyperparameter σ _d . By increasing σ _d the overlapping region of both regularization matrices is reduced. It is worth noticing that the case of strict distinction between both regularization matrices explained in Section 4.1 is included in this optimization problem for σ _d → ∞.

Figure 3:

Sigmoid functions for fuzzification of kernel superim-position.

Now, the hyperparameter optimization from Section 3.1 can be utilized to optimize d ̂ and σ _d along with the other kernel specific hyperparameters. An intermediate version of this is approach is achieved by presetting σ _d and only optimizing d ̂ in the hyperparameter optimization. To distinguish between these two versions, in the following we call the full optimization SK₂ and the intermediate optimization SK₁.

4.3 Multiple kernel method

The second method to incorporate time delay into the procedure of identifying RFIR models is based on a multiple kernel method. Here, multiple time-delay kernels are used for the identification of time-delay processes. Thereby, the mixed-integer optimization problem discussed in Section 4.2 can be overcome. The time-delay kernel R ̲ IRDP in Section 4.1 with different time delays δ ̲ = d min , d min + 1 , … , d max is used. The number of n _MK = d _max − d _min + 1 regularization matrices R ̲ IRDP are each multiplied with a separate λ _i (all gathered in λ ̲ ) and summed up

The considered time delays and therefore also the number of multiple kernels n _MK as well as the parameter β have to be set a priori. The model parameters can be calculated with

In comparison to (8), a vector of λ ̲ is used instead of a scalar λ. This leads to more hyperparameters. To determine them, a hyperparameter optimization as described in Section 3.1.2 is performed. For the calculation of the scalar estimated time delay d ̂ the time delays of the different kernels are weighted with the corresponding regularization strengths λ _i

4.4 Comparison of the methods

In this section a kernel, which incorporates time delay in the structural assumptions of the regularization matrix, was proposed. To include an unknown time delay in the RFIR estimation and avoiding a mixed-integer optimization problem, two methods were developed. The complexity of the utilized hyperparameter optimization increases by the proposed methods. Compared to the estimation with the regular IRP kernel, method SK₁ requires one additional hyperparameter (inflection point) and SK₂ requires two additional hyperparameters (inflection point and slope). The number of hyperparameters required for the MK method depends on the number of kernels which has to be chosen a priori. This requires the selection of a time interval the time delay is expected in. The higher the uncertainty of the unknown time delay, the more hyperparameters have to be optimized. Therefore, the MK method requires n _MK − 1 more hyperparameters than the estimation with the regular IRP kernel. The kernels as well as the hyperparameters used in the different methods are summarized in Table 1.

5.1 Simulation setup

To investigate the different methods regarding robustness of the time-delay estimation and the model quality, different Monte Carlo studies with 100 runs are performed. In Table 2 three different process types are defined. The pole(s) of the different processes are chosen uniformly distributed in the given range for each run. The range is determined to cover typically present poles with a well chosen sample time.

White Gaussian noise with N = 1000 data samples is used as excitation signal. In each run other noise realizations are used for the excitation. White Gaussian noise is added to the output data such that a signal-to-noise-ratio (SNR) of 10 dB results. The FIR model order n is set to 100 in order to capture all significant impulse response coefficients. For the MK method n _MK = 21 kernels with δ ̲ = 5,6 , … , 25 are used. To evaluate the model quality, the fit between the true impulse response coefficients θ ̲ * , which we assume to be the first n coefficients of g ̃ ̲ , and the estimated model parameters θ ̂ ̲ is calculated by using

5.2 Kernel specific hyperparameter choice

The time-delay kernel explained in Section 4.1 contains the hyperparameter β which accounts for numerical stability in the calculation of the parameters θ ̂ ̲ . To make a reasonable choice on this parameter, a Monte Carlo study as described in Section 5.1 for each process order from Table 2 is performed. The parameter β is varied in the heuristically chosen range of

The order of the IRDP kernel, which is part of the time-delay kernel, is varied from first to third order.

For the SK₁ method the results show, that the mean model fit is not sensitive w.r.t. the choice of β. The same holds for the estimation accuracy of the time-delay. However, it could be observed that the smaller β is chosen the more outlier are present.

Table 3 summarizes the results of the investigation on the SK₂ and MK method w.r.t. the choice of β. It can be seen that the optimal choice of β for the SK₂ method is dependent on the process order. For first-order processes G ₁(z ⁻¹) higher values for β yield the best results both, for fit and accuracy of the time-delay estimation. For higher-order processes such as G ₂(z ⁻¹) and G ₃(z ⁻¹) lower values for β show better results. Against that, the order of the IRDP kernel does not influence the best choice of β which can be seen from identical best values for β in each column.

Exemplary, in the upper part of Figure 4 the parameter fit calculated with (30) is shown for 100 input and noise realizations for second-order processes, for the different choices on β and the SK₂ IRDP₂ kernel. The distribution of the fit values shows that the estimation method for the given process is robust for β ≥ 10⁻⁴. For lower values the estimation becomes less robust which can be substantiated by the increasing number of outlier models. Similar results are also observed for the other process orders.

Figure 4:

Boxplot of parameter fit for different choices of β for second-order processes and the methods SK₂ and MK with the IRDP₂ kernel for 100 realizations.

The MK method performs almost independently from the process order and the order of the IRDP kernel and yields the best results with high values for β such as 10⁻¹ or 10°, respectively. From the lower part of Figure 4 it can be reasoned that the robustness and performance is significantly reduced for β ≤ 10⁻². This behavior can also be observed for the other process and IRDP kernel order combinations.

To summarize the results, we advise to choose β = 10° for the SK₁ method. A choice of β dependent on the process is advised for the SK₂ method. If this knowledge is not available, β = 10⁻¹ achieves a good trade-off between robustness and performance of the estimation. For the MK method choosing β = 10⁻¹ generally yields good estimation performance and robustness simultaneously and is recommended.

5.3 Criteria for hyperparameter optimization

The choice of the criterion for the hyperparameter optimization is crucial, since the model quality strongly depends on the hyperparameter. Therefore, the marginal likelihood (ML) function and generalized cross-validation (GCV) error explained in Section 3.1.2 are investigated with a Monte Carlo study described in Section 5.1 and the resultant parameter fit as well as the time-delay estimation are compared. For this purpose only the results of the second-order process G ₂(z ⁻¹) from Table 2 are shown, since the results of the other process orders are similar.

For the investigation of the different hyperparameter optimization criteria an interior-point optimization algorithm is utilized. It is chosen based on a comparison of different optimization algorithms, in which it outperformed the competitive algorithms.

In Figure 5 a boxplot of the parameter fit and the time-delay estimation for the three different methods (SK₁, SK₂, and MK) with the IRDP₂ kernel is shown. The parameter β is chosen according Table 3. Here, both the mean and the variance of the parameter fit are significantly better using the GCV error instead of ML for optimization. Similar to the fit, the GCV error outperforms the ML and consequently is the preferred choice for the hyperparameter optimization objective. For the multiple kernel (MK) approach with ML, it can be seen that there is a mean error of −0.5 for the time-delay estimation. This means that by using ML in the optimization not only one kernel is weighted highly, but two consecutive kernels with a similar regularization strength are chosen. This undesirable choice also leads to a worse parameter fit.

Figure 5:

Boxplot of parameter fit and time-delay estimation for the three different methods SK₁, SK₂, and MK with the IRDP₂ kernel for 100 realizations of G ₂(z ⁻¹). Here, the two criteria – generalized cross-validation (GCV) error and the marginal likelihood (ML) function are compared.

In summary, for all methods the models with GCV error optimized hyperparameters show better and more robust performance. Therefore, in the following only the GCV error is used for hyperparameter optimization. This result deviates from previous investigations on the TC or IRP kernel for linear dynamic system without time-delay. In that case the ML objective usually outperforms the GCV objective [25, 26].

5.4 Comparison of all proposed methods

In this section, all proposed methods for identifying time-delay processes (SK₁, SK₂, MK) are compared. To show the improvement of the novel methods developed for time-delay processes, also unregularized FIR models and RFIR models with the regular IRP kernel are compared. The order of the IRP and IRDP kernel, respectively are chosen identically to the order of the process. The hyperparameter optimization is performed with the GCV error and β is set according to Table 3. All three process orders of Table 2 are investigated. A Monte Carlo simulation with the setup of Section 5.1 is performed and the parameter fit as well as the deviation of the estimated time delay d ̂ to the true time delay d is compared. It has to be noted, that for unregularized FIR models as well as for RFIR models with the regular IRP kernel no specific time-delay estimation procedure is carried out and therefore no explicit information can be gained.

In Figure 6 the investigation of first-order processes is shown. It can be seen that unregularized FIR models and RFIR models with IRP kernel show a similar performance. Due to the inappropriate regularization matrix, the regularization strength has been chosen low. The specifically developed methods for time-delay processes show a significant improvement of the estimation. Especially, the multiple kernel (MK) method outperforms the other methods in the parameter fit as well as in the accuracy of the time-delay estimation. For the time-delay estimation SK₂ shows a smaller 75 % quantile, but yields more outlier models in comparison to SK₁.

Figure 6:

Boxplot of parameter fit and time-delay estimation for the five different methods FIR, IRP, SK₁, SK₂, and MK with the IRP₁ and IRDP₁ kernel for 100 realizations of G ₁(z ⁻¹). Note: There is no time-delay estimation for the methods FIR and IRP.

The results for a second-order process are shown in Figure 7. It again shows the superiority of the MK method. For second-order processes the performance of the IRP method is superior in comparison to the unregularized FIR models. This can be explained by the smoothness assumption of the IRP kernel since the impulse response of a first-order process is not continuous. However, the impulse response of a second-order process is continuous and therefore the smoothness assumption is not violated. Consequently, this allows to choose the regularization strength λ higher. For this process order, the method with variable slope SK₂ obtains better results than with a fixed slope. Here, all deviations of the time-delay estimation of the MK method lie in a range of 0– and −1.

Figure 7:

Boxplot of parameter fit and time-delay estimation for the five different methods FIR, IRP, SK₁, SK₂, and MK with the IRP₂ and IRDP₂ kernel for 100 realizations of G ₂(z ⁻¹).

For third-order processes the results of the parameter fit as well as the time-delay estimation are similar to the second-order process. Only the time-delay estimation of method SK₁ is slightly inferior and is not as robust as for the second-order process.

In Figure 8 estimated impulse responses of a second-order process G ₂(z ⁻¹) with the fixed poles a = 0.9, b = 0.7, and the time delay d = 10 are shown. The process is excited and disturbed as explained in Section 5.1. The parameter variance within the different sets of FIR models corresponding to the chosen estimation methods show the superiority of the novel methods. Especially in the time-delay part of the impulse response almost every deviation from zero can be suppressed. Additionally, improvements in the smoothness are also visible in the dynamics part. Since the results of the methods SK₁ and SK₂ are similar, only SK₂ is shown for improved visualization. Moreover, the MK method outperforms the SK₂.

Figure 8:

Impulse responses of a second-order process G ₂(z ⁻¹) with poles a = 0.9, b = 0.7, and time delay d = 10 for 50 different noise realizations estimated with the methods FIR, IRP, SK₂, and MK with the IRP₂ and IRDP₂ kernel.

In summary, the multiple kernel (MK) method outperforms all other methods in both, robustness and performance of the estimation. However, the number of hyperparameters increases with the number of time delays assumed a priori. Nevertheless, the hyperparameter optimization achieves a sparse configuration of λ ̲ even with n _MK = 21, i.e., most λ _i = 0, typically only one λ _i ≠ 0, for i = 1, 2, …, n _MK. Moreover, the incorporation of the time delay in the structure of the kernel significantly improves the model estimation of time-delay processes in comparison to the standard techniques. The methods SK₁ and SK₂ show similar model qualities, although for second- or third-order processes optimizing the slope of the sigmoids improves the estimation performance. In contrast to regular regularization methods like the TC or IRP kernel, for the proposed methods the GCV error should be used as criteria for the hyperparameter optimization instead of ML to get the best model performance.