Useful redundancy in parameter and time delay estimation for continuous-time models

TL;DR

This paper introduces an algorithm for estimating parameters and time delay in continuous-time systems using instrumental variable methods, enhanced by useful redundancy to ensure convergence to the global minimum.

Contribution

It presents a novel approach combining useful redundancy with filtering to improve the convexity of the cost function in time delay estimation.

Findings

Effective convergence to the global minimum demonstrated

Algorithm outperforms traditional methods in numerical tests

Improved robustness in parameter and delay estimation

Abstract

This paper proposes an algorithm to estimate the parameters, including time delay, of continuous time systems based on instrumental variable identification methods. To overcome the multiple local minima of the cost function associated with the estimation of a time delay system, we utilise the useful redundancy technique. Specifically, the cost function is filtered through a set of low-pass filters to improve convexity with the useful redundancy technique exploited to achieve convergence to the global minimum of the optimization problem. Numerical examples are presented to demonstrate the effectiveness of the proposed algorithm.

Tables

Table 1. Table 1: Global convergence as a percentage of 100 Monte Carlo simulations for case 1

Sampling scheme		Existing method [12]						Proposed method
	Initial delay	0s	1s	3s	5s	7s	9s	0s	1s	3s	5s	7s	9s
Regular sampled data	SNR = 5dB	65%	100%	100%	97%	100%	100%	100%	100%	100%	100%	100%	100%
	SNR = 15dB	69%	89%	92%	99%	100%	100%	100%	100%	100%	100%	100%	100%
Irregular sampled data	SNR = 5dB	1%	1%	0%	17%	100%	90%	100%	100%	100%	100%	100%	100%
	SNR = 15dB	1%	0%	7%	18%	100%	89%	100%	100%	100%	100%	100%	100%

Table 2. Table 2: Global convergence as a percentage of 100 Monte Carlo simulations for case 1 with a randomly selected initial delay

Sampling scheme	SNR	Existing method [12]	Proposed method
Regular sampled data	5dB	90%	100%
	15dB	97%	100%
Irregular sampled data	5dB	40%	100%
	15dB	41%	100%

Table 3. Table 3: Global convergence as a percentage of 100 Monte Carlo simulations for Case 2

Sampling scheme		Existing method [12]						Proposed method
	Initial delay	0s	1s	3s	5s	7s	9s	0s	1s	3s	5s	7s	9s
Regular sampled data	SNR = 5dB	1%	52%	22%	43%	76%	68%	99%	99%	99%	99%	99%	99%
	SNR = 15dB	0%	60%	6%	28%	70%	69%	100%	100%	100%	100%	100%	100%
Irregular sampled data	SNR = 5dB	8%	16%	23%	26%	66%	48%	100%	100%	100%	100%	100%	100%
	SNR = 15dB	9%	10%	19%	27%	55%	50%	100%	100%	100%	100%	100%	100%

Table 4. Table 4: Global convergence as a percentage of 100 Monte Carlo simulations for case 2 with a randomly selected initial delay

Sampling scheme	SNR	Existing method [12]	Proposed method
Regular sampled data	5dB	52%	99%
	15dB	46%	100%
Irregular sampled data	5dB	48%	100%
	15dB	21%	100%

Equations

$$x(t)=G_0(p)u(t-{\tau }_0)=\frac{B(p)}{A(p)}u(t-{\tau }_0),$$

$$B(p)$$

$$A(p)$$

$$y(t_k)=x(t_k)+e(t_k).$$

$$h_k=t_{k+1}-t_k,k=1,2,...,N-1,$$

$$\frac{A(p)}{A(p)}y(t_k)=\frac{B(p)}{A(p)}u(t_k-{\tau }_0)+e(t_k),$$

$$A(p)y_A(t_k)=B(p)u_A(t_k-{\tau }_0)+e(t_k),$$

$$y_A(t_k)=\frac{1}{A(p)}y(t_k),u_A(t_k-{\tau }_0)=\frac{1}{A(p)}u(t_k-{\tau }_0).$$

$$y_A^{(n)}(t_k)={\phi }_N(t_k,{\tau }_0)\theta +e(t_k)$$

$${\phi }_N(t_k)=[-y_A^{(n-1)}(t_k),-y_A^{(n-2)}(t_k),...,-y_A(t_k),$$

$$u_A^{(m)}(t_k-{\tau }_0),u_A^{(m-1)}(t_k-{\tau }_0)...u_A(t_k-{\tau }_0)],$$

$$\theta =[a_1,a_2,...,a_n,b_0,b_1,...,b_m{]}^T,$$

$$y_A^{(i)}(t_k)=\frac{d^i{y}_A(t_k)}{d{t}^i},u_A^{(i)}(t_k-{\tau }_0)=\frac{d^i{u}_A(t_k-{\tau }_0)}{d{t}^i}.$$

$${\hat{\phi }}_N(t_k)=[-{\hat{x}}_A^{(n-1)}(t_k),-{\hat{x}}_A^{(n-2)}(t_k),...,-{\hat{x}}_A(t_k),$$

$$u_A^{(m)}(t_k-{\tau }_0),u_A^{(m-1)}(t_k-{\tau }_0)...u_A(t_k-{\tau }_0)]$$

$${\hat{x}}_A(t_k)=\frac{1}{\hat{A}(p)}\hat{x}(t_k);\hat{x}(t_k)=\frac{\hat{B}(p)}{\hat{A}(p)}u(t_k-{\tau }_0),$$

$$u_A(t_k-\tau )=\frac{1}{A(p)}u(t_k-\tau ),$$

$$\{\begin{array}{l}\dot{\text{z}}(t-\tau )={\text{F}}_c\text{z}(t-\tau )+{\text{G}}_{c}u(t-\tau )\\ u_A(t-\tau )={\text{H}}_c\text{z}(t-\tau ).\end{array}$$

$$\{\begin{array}{l}\text{z}(t_{k+1}-\tau )=\text{F}(h_k)\text{z}(t_k-\tau )+\text{M}(t_k)\\ u_A(t_k-\tau )={\text{H}}_c\text{z}(t_k-\tau )\end{array}$$

$$\text{F}(h_k)=e^{{\text{F}}_c{h}_k},$$

$$\text{M}(t_k)=\sum _{i=1}^{n_k-1}\prod _{j=i+1}^{n_k-1}\text{F}(h_{k,j})\text{G}(h_{k,i})u(t_{k,i}),$$

$$\text{F}(h_{k,i})=\text{F}(m_{k,i}\Delta )\text{F}({\delta }_{k,i}),$$

$$\text{G}(h_{k,i})=\text{F}(m_{k,i}\Delta )\text{G}({\delta }_{k,i})+\text{G}(m_{k,i}\Delta ).$$

$$F(p)=\frac{1}{(p+{\omega }_c^{SVF}{)}^n},$$

$${\hat{\theta }}_0=[{\Phi }_N{\Phi }_N^T{]}^{-1}{\Phi }_N{Y}_N.$$

$${\Phi }_N=[{\phi }_N^T(t_1)...{\phi }_N^T(t_N)],$$

$$Y_N=[y_A^{(n)}(t_1)y_A^{(n)}(t_2)...y_A^{(n)}(t_N){]}^T,$$

$$\hat{x}(t_k)=\frac{\hat{B}(p,{\hat{\theta }}^{j-1})}{\hat{A}(p,{\hat{\theta }}^{j-1})}u(t_k-{\tau }_0).$$

$$F_c(p)=\frac{1}{\hat{A}(p,{\hat{\theta }}^{j-1})},$$

$${\hat{\theta }}^j=[{\Psi }_N{\Phi }_N^T{]}^{-1}{\Psi }_N{Y}_N.$$

$${\Phi }_N=[{\phi }_N^T(t_1)...{\phi }_N^T(t_N)],$$

$${\Psi }_N=[{\hat{\phi }}_N^T(t_1)...{\hat{\phi }}_N^T(t_N)],$$

$${\bar{y}}_A^{(n)}(\rho )={\bar{\phi }}_N(\rho )\theta +\overline{ϵ}(\rho ),$$

$$\hat{\theta }(\tau )=[{\bar{\Psi }}_N(\rho ){\bar{\Phi }}_N^T(\rho ){]}^{-1}{\bar{\Psi }}_N(\rho ){\bar{Y}}_N,$$

$$\hat{\tau }=\underset{\tau }{\text{arg min}}\widetilde{\stackrel{ˉ}{J_N}}(\tau )=\underset{\tau }{\text{arg min}}\overline{J_N}(\rho ){|}_{\theta =\hat{\theta }(\tau )},$$

Full text

††thanks: The material in this paper was not presented at any conference.

, , ,

1 Introduction

The goal of system identification is to estimate the parameters of a model in order to analyse, simulate and/or control a system. In the time domain, there are two typical approaches to identify a system, i.e. discrete-time (DT) identification and continuous-time (CT) identification. For several decades, DT identification has been dominant due to the strong development of the digital computer. More recently, estimation using continuous-time identification methods has received much attention due to advantages such as providing insights to the physical system and being independent of the sampling time [15] [17] [25]. For example, with irregular sampling time, the DT model becomes time-varying; hence the DT system identification problem becomes more difficult while for CT identification, the system is still time-invariant. In reality, irregular sampling occurs in many cases, e.g. when the sampling is event-triggered, when the measurement is manual and also in the case of missing data [5].

In some CT identification problems, one wants to estimate the system parameters and any unknown time delay, as there exist many practical examples including, chemical processes, economic systems, biological systems, that possess time delays. It is important to estimate the delay accurately since a poor estimate can result in poor model order selection and inaccurate estimates of the system parameters.

There are many approaches to estimate a system time delay [9]. A simple approach is to consider the impulse response data, e.g. estimate the time delay by finding where the impulse response becomes nonzero [10] or by noting the delay where the correlation between input and output is maximum [10] [11]. Another approach is to model the delay by a rational polynomial transfer function using a Padé or similar approximation and then estimate the time delay as part of the system parameters [1] [2] [18]. In [27] [6] [7], the time delay and system parameters of a Multiple Input Single Output CT system are estimated in a separable way using an iterative global nonlinear least-squares or instrumental variable method.

Recently, a method of estimating the parameters and time delay of CT systems has been suggested in [12] which is based on a gradient technique. The parameters and the time delay are estimated separately, i.e. when one is estimated, the other is fixed which is then repeated in an iterative manner. In this approach, the Simplified Refined Instrumental Variable (SRIVC) method is used to estimate the parameters whilst the time delay is estimated using the Gauss-Newton method. In addition, due to the effects of multiple minima in the cost function to be minimized for the time delay [20] [23] [8], a low-pass filter is employed to increase convexity. As shown in [8] [14] [13], a suitable low-pass filtering operation on the estimation data can help to extend the global convergence region of the cost function, hence improve the accuracy of the time delay estimate.

In this paper, we adopt the idea of using a low-pass filter where instead of using only one filter, we suggest to use multiple low-pass filters and incorporate the useful redundancy technique [3] [4]. Useful redundancy is a technique that was developed to avoid local minima when solving a nonlinear inverse problem. The concept is to generate a family of cost functions that have different local minima, but share the same global minimum with the original cost function of the optimization problem. Whenever the algorithm is trapped in a minimum, the solver path is switched to another solver path in a way such that the minimum found using the new solver path corresponds to a decrease in the original cost function. This allows the algorithm to cross local minima and converge to the global minimum, hence improving the accuracy of the estimated parameters. For the algorithm described in this paper, the multiple cost functions are generated by filtering the original time delay cost function through a number of low-pass filters with different cut-off frequencies that span the system bandwidth.

The paper is organized as follows. Section 2 describes the model setting and Section 3 recalls both the SRIVC algorithm and the SRIVC-based time delay estimation using a low-pass filter. Section 4 formulates the algorithm of the new method and provides analysis on the effectiveness of the proposed method. Section 5 presents numerical experiments and results for both regular and irregular sampling schemes. Finally, the conclusion will be drawn in Section 6.

2 Model setting

Consider a continuous-time linear, time invariant, single input single output system,

[TABLE]

with

[TABLE]

where $\tau_{0}$ is the time delay, $u(t),\ x(t)$ are the input and deterministic output of the system respectively and $p$ is the differential operator, i.e. $p^{(i)}x(t)=d^{i}x(t)/dt^{i}$. In addition, the following assumptions are made:

Assumption 2.1

Polynomials $B(p)$ and $A(p)$ are coprime.

Assumption 2.2

The system is asymptotically stable.

Assumption 2.3

The high frequency gain of $G_{0}(p)$ is 0, i.e. $G_{0}(p)$ is strictly proper.

The deterministic output $x(t_{k})$ is measured as $y({t_{k}})$ in the presence of noise, i.e.

[TABLE]

Furthermore, we consider the sampling time of the input $u(t_{k})$ and output data $y(t_{k})$ as either regular or irregular. The time-varying sampling interval is denoted as,

[TABLE]

where $N$ is the length of the data.

The objective of a CT system identification problem is to estimate the time delay, $\tau_{0}$, and the parameters $a_{1},a_{2},...,a_{n},b_{0},b_{1},...,b_{m}$ of the CT model in (1), using the measured input and output data, $u(t_{k})$ and $y(t_{k})^{N}_{k=1}$ respectively.

3 Parameter and time delay identification of Continuous-time models with SRIVC and filtering

There exists a large number of continuous time identification methods see, for example, [17] [28]. In this paper, we consider the Simplified Refined Instrumental Variable Continuous Time method (SRIVC), which is developed in the literature by Young and Garnier [16][17][28]. We begin with a basic description of the SRIVC method.

3.1 Traditional SRIVC method

To simplify the description we first assume the time delay $\tau_{0}$ is known. The SRIVC method can then be summarized as follows.

From (1) and (2),

[TABLE]

which becomes,

[TABLE]

with

[TABLE]

From (5), a linear regression model can be formed as

[TABLE]

where

[TABLE]

To solve the linear regression problem in (7), the SRIVC algorithm uses an instrumental variable (IV) method. The instruments in the SRIVC method are chosen as the estimated noise free outputs, i.e. the regressor becomes,

[TABLE]

where

[TABLE]

with $\hat{B}(p)$ and $\hat{A}(p)$ estimates of $B(p)$ and $A(p)$.

3.2 Implementation of the CT filtering operation for irregular sampled data and arbitrary time delay

The SRIVC method can be used to estimate system parameters from regular sampled data as well as from irregular sampled data. However, one difficulty in implementing the SRIVC method for irregular sampled data in the presence of an arbitrary time delay is the CT filtering operation, e.g.,

[TABLE]

when $\tau$ is an arbitrary time delay.

There are two reasons for this difficulty:

$u(t_{k}-\tau)$ is not available from the measured data, 2. 2.

The digital simulation of $u(t_{k}-\tau)$ is generally performed in state-space form hence, an equivalent discrete-time state space representation of the CT state space model needs to be computed. For irregular data, as the sampling interval $h(t_{k})$ is time-varying, the computational load using standard methods, e.g. expm in Matlab, to compute the transformation matrices will be large.

As suggested in [12], the two problems mentioned above can be solved as follows,

$u(t_{k}-\tau)$ can be constructed from the neighbouring data based on the inter-sample behaviour, e.g. zero-order-hold (ZOH), first-order-hold (FOH), etc. 2. 2.

To reduce the computational load, the time-varying sampling interval is divided into two intervals: the first interval is a multiple of a constant sampling period, the second interval is the residual. The discretization matrix of the first interval is pre-computed using a standard method and stored in an array. The discretization matrix for the second interval is computed by a fast approximation, e.g. the 4th order Runge-Kutta (RK4) method. The final discretization matrix of the sampling interval $h_{k}$ is a product of the two matrices.

Specifically, convert the transfer function (TF) model in (11) to the equivalent CT state-space model,

[TABLE]

The equivalent DT state-space model will be,

[TABLE]

where $\textbf{F}(h_{k})$, $\textbf{M}(t_{k})$ are computed as follows,

1. (a)

Denote $t_{k,1},t_{k,n_{k}}$ ($n_{k}\geq 2$) as $t_{k}-\tau$, $t_{k+1}-\tau$ respectively; $t_{k,i},i=2,...,n_{k}-1$ as the transition time-instants of $u(t)$ between $t_{k}-\tau$ and $t_{k+1}-\tau$; and $h_{k,i}=t_{k,i+1}-t_{k,i}$. 2. (b)

Let $\Delta$ be a constant sampling period; $m_{k,i}$ be a positive integer; $\delta_{k,i}\geq 0$ such that:

$h_{k,i}=m_{k,i}\Delta+\delta_{k,i}<(m_{k,i}+1)\Delta$. 3. (c)

Compute $\textbf{F}(m\Delta),\textbf{G}(m\Delta)$. 4. (d)

Compute $\textbf{F}(\delta_{k,i}),\textbf{G}(\delta_{k,i})$ using, e.g. RK4. 5. (e)

Finally, compute $\textbf{F}(h_{k})$, $\textbf{M}(t_{k})$ as follows,

[TABLE]

with

[TABLE]

Now that we can generate filtered signals of the irregular sampled data, we describe the traditional SRIVC algorithm [16][17][28] in Algorithm I.

ALGORITHM I

Step 1. Initialization

Create the stable state variable filter (SVF),

$F(p)=\dfrac{1}{(p+\omega_{c}^{SVF})^{n}},$

(14)

where $\omega_{c}^{SVF}$ is chosen to be larger than or equal to the bandwidth of the system.

Filter $y(t_{k})$ and $u(t_{k}-\tau_{0})$ via the SVF to generate derivatives of the signals, i.e. $y_{A}^{(n)}(t_{k}),y_{A}^{(n-1)}(t_{k}),...,$ $y_{A}(t_{k})$ and $u_{A}^{(m)}(t_{k}-\tau_{0}),u_{A}^{(m-1)}(t_{k}-\tau_{0}),...,u_{A}(t_{k}-\tau_{0})$ by using the method described in [16].

Use the least squares method to estimate the initial parameters,

$\hat{\theta}_{0}=[\Phi_{N}\Phi_{N}^{T}]^{-1}\Phi_{N}Y_{N}.$

(15)

with

$\displaystyle\Phi_{N}=[\varphi_{N}^{T}(t_{1})\ ...\ \varphi_{N}^{T}(t_{N})],$

(16)

$\displaystyle Y_{N}=[y_{A}^{(n)}(t_{1})\ y_{A}^{(n)}(t_{2})\ ...\ y_{A}^{(n)}(t_{N})]^{T},$

where $\varphi_{N}$ is defined in (8).

Step 2. Iterative estimation

for j = 1:convergence111Convergence requires the relative error between the estimated parameters of two consecutive iterations to be less than an $\epsilon>0$ where $\epsilon$ is a small value chosen to achieve the desired accuracy.

If the estimated TF model is unstable, reflect the right half plane zeros of $\hat{A}(p,\hat{\theta}^{j-1})$ into the left half plane and construct the new estimate.

Generate the instrumental variable,

$\hat{x}(t_{k})=\dfrac{\hat{B}(p,\hat{\theta}^{j-1})}{\hat{A}(p,\hat{\theta}^{j-1})}u(t_{k}-\tau_{0}).$

(17)

Filter the input $u(t_{k}-\tau_{0})$, output $y(t_{k})$ and the instrumental variable $\hat{x}(t_{k})$ by the continuous time filter

$F_{c}(p)=\dfrac{1}{\hat{A}(p,\hat{\theta}^{j-1})},$

(18)

to generate the derivatives of these signals.

Using the prefiltered data, generate an estimate $\hat{\theta}^{j}$ using the IV method,

$\hat{\theta}^{j}=[\Psi_{N}\Phi_{N}^{T}]^{-1}\Psi_{N}Y_{N}.$

(19)

where $\Psi_{N}$ is the IV matrix generated by the instrumental variables $\hat{x}(t_{k})$ and $\Phi_{N}$ is the regression matrix,

$\displaystyle\Phi_{N}=[\varphi_{N}^{T}(t_{1})\ ...\ \varphi_{N}^{T}(t_{N})],$

(20)

$\displaystyle\Psi_{N}=[\hat{\varphi}_{N}^{T}(t_{1})\ ...\ \hat{\varphi}_{N}^{T}(t_{N})],$

where $\hat{\varphi}_{N}$ is defined in (9).

end

3.3 SRIVC-based time delay estimation with filtering

A recently developed method [12] to estimate both the time delay and parameters for the CT model (2) considers it as a separable nonlinear least squares problem. The SRIVC algorithm and the Gauss-Newton method are used in [12] to estimate the system parameters and the time delay respectively. In this problem the cost function for the time delay estimation has multiple minima [12], hence a low-pass filter is utilized to extend the global convergence region [12] [13] [14].

Let $L(p)$ be a CT low-pass filter with cut off frequency $\omega_{c}^{L}$. Applying the filter $L(p)$ to the linear regression (7), we obtain

[TABLE]

where $\rho=[\theta^{T},\tau]$ and $\bar{.}$ represents the signal $y,u$ or $e$ filtered by $L(p)$.

When the estimated parameter $\hat{\theta}$ is a function of $\tau$, i.e.

[TABLE]

then the time delay $\tau$ can be estimated as [12][21],

[TABLE]

with

[TABLE]

where $s$ is chosen to guarantee that at sample $s$, $t_{s}\geq\tau$.

By using the Gauss-Newton method, the system parameters and the time delay can be iteratively estimated,

[TABLE]

where $\mu^{j}$ is the step size and

[TABLE]

The SRIVC-based time delay estimation with filtering [12] is summarized in Algorithm II.

ALGORITHM II

Step 1. Initialization

Set the initial value $\hat{\tau}^{0}$, the boundaries222$\tau_{min}$ and $\tau_{max}$ are boundaries for the delay, which are known a priori, i.e. $\tau_{min}\leq\tau_{0}\leq\tau_{max}$ $\tau_{min},\tau_{max}$, $\Delta\tau_{min}$333$\Delta\tau_{min}$ is the limit value of the step size of the time delay estimate using Gauss-Newton method., the cut-off frequency $\omega_{c}^{L}$ of $L(p)$, the cut-off frequency of $F(p)$, $\omega_{c}^{SVF}$, and a small positive $\epsilon$ to indicate convergence.

Based on the initial value $\hat{\tau}^{0}$, use the SRIVC method to compute $\hat{\theta}^{0}$.

Compute $\nabla^{2}\tilde{\bar{J_{N}}}(\hat{\tau}^{0})$ and $\nabla\tilde{\bar{J_{N}}}(\hat{\tau}^{0})$ from (26).

Step 2. Iterative estimation

for j = 1:convergence

Compute $\Delta\hat{\tau}^{j}$ using the filtered input/output data, equation (26) and

$\Delta\hat{\tau}^{j}=-[\nabla^{2}\tilde{\bar{J_{N}}}(\hat{\tau}^{j})]^{-1}\nabla\tilde{\bar{J_{N}}}(\hat{\tau}^{j}).$

(27)

(a)

Compute $\hat{\tau}^{j+1}=\hat{\tau}^{j}+\Delta\hat{\tau}^{j}$

If $\hat{\tau}^{j+1}\not\in[\tau_{min},\tau_{max}]$, let $\Delta\hat{\tau}^{j}=\Delta\hat{\tau}^{j}/2$ and repeat this step.

If $|\Delta\hat{\tau}^{j}|\leq\Delta\tau_{min}$, break.

(b)

Estimate $\hat{\theta}^{j+1}$ using the SRIVC method and time delay $\hat{\tau}^{j+1}$ from the filtered input/output data.

(c)

Compute $\tilde{\bar{J_{N}}}(\hat{\tau}^{j+1})$. If $\tilde{\bar{J_{N}}}(\hat{\tau}^{j+1})\geq\tilde{\bar{J_{N}}}(\hat{\tau}^{j})$, let $\Delta\hat{\tau}^{j}=\Delta\hat{\tau}^{j}/2$ and return to (a).

If $\tilde{\bar{J_{N}}}(\hat{\tau}^{j})-\tilde{\bar{J_{N}}}(\hat{\tau}^{j+1})\geq\epsilon$, go to Step 1, else break.

end

Step 3. Refined parameter estimation

Repeat Step 2 with low-pass filter $L(p)=1$.

Remark 3.1

The smaller the cut-off frequency $\omega_{c}^{L}$, the larger the global convergence region. However, when $\omega_{c}^{L}$ is chosen too small, information from the data is lost, which affects the accuracy of the parameter estimation [12]. The suggestion in [12] is,

[TABLE]

where bw is the bandwidth of the system $G_{0}(p)$ (rad/sec).

4 SRIVC-based time delay estimation with multiple filtering and redundancy method

In this section, we propose a method that improves convergence to the global minimum of the time delay optimization problem using the useful redundancy technique [3][4]. The useful redundancy technique was developed to avoid local minima when solving a non-linear inverse problem. As stated in the previous section, the main problem in time delay estimation is that the cost function $J_{N}(\rho)$ possesses many local minima [12]. The filtering operation [14] described in Section 3.3 can be employed in order to increase the global convergence region of the cost function associated with the time delay estimation. However, when the data is very noisy or the initial value, $\hat{\tau}^{0}$, is located far from the global minimum, using only one filter does not guarantee that the solution of the optimization problem converges to the global minimum.

To demonstrate how useful redundancy can be utilized in our estimator, we first describe the useful redundancy technique.

4.1 The useful redundancy technique

We define the useful redundancy technique by quoting directly from [4].

Definition 1 [4]. Consider an optimization problem,

[TABLE]

Then it is called M-safely redundant if and only if the following conditions hold:

There exists a finite $M$ cost functions $J_{i}$ sharing the same global minimum $\rho^{*}\in\mathbb{R}^{n}$. 2. 2.

There exists a solver (or an iterative scheme) $\ell$ and a finite number of iterations $r^{*}\in\mathbb{N}$ such that for some $\gamma\in[0,1)$ and all $\rho\in\mathbb{P}$ the following inequality holds,

[TABLE]

where $\ell^{(r^{*})}(\rho,J_{i})$ is the candidate solution obtained after $r^{*}$ iterations of $\ell$ using the cost function $J_{i}$ and starting from the initial guess $\rho$. $\hfill\square$

The solver path $(\rho,J_{i})$ is defined as the sequence of iterates $\ell^{(j)}(\rho,J_{i})$ for the solver $\ell$ when the cost function $J_{i}$ starts from an initial guess $\rho$. Condition 2 means that for any initial $\rho$, there always exists a solver path $\ell$ that corresponds to a decrease in the original cost $J_{0}$ after at most $r^{*}$ iterations. It is proven in [4] that if an optimization problem is M-safely redundant following from Definition 1, that convergence to a global minimum is guaranteed.

An algorithm that describes the useful redundancy technique is given in Algorithm III.

ALGORITHM III

Step 1. Initialization Choose an initial value $\rho^{(0)}$.

Step 2. Iterative estimation

for $k$=0:converge

Set $i=0$ and $failure=true$

While $failure$ do

(a)

Use the solver path $(\rho,J_{i})$ starting from $\rho^{(k)}$, find the corresponding minima $\xi^{(k,i)}$.

(b)

Compute

failure=\Big{(}J_{0}(\xi^{(k,i)})>(1-\epsilon)J_{0}(\rho^{(k)})\Big{)} 444$\epsilon$ is a predefined small value that is chosen based on the desired accuracy.

(c)

If $failure$ then $i=(i+1)\ \text{mod}\ M$,

Else set $k=k+1,\rho^{(k)}=\xi^{(k,i)}$.

End while

end

4.2 Theoretical results related to the time delay cost function

To construct an M-safely redundant optimization problem for the estimation of the time delay, we need a cost function $J_{0}$ and multiple solver paths that satisfy the two conditions in Definition 1. Here, the solver paths are generated by filtering the time delay estimation error using a set of low-pass filters with different cut-off frequencies. The cost function $J_{0}$ is formulated from these filtered cost functions such that there always exists a solver path that corresponds to a decrease in $J_{0}$ after a finite number of iterations. Next we describe the set of filters and the cost function $J_{0}$ required to satisfy the conditions in Definition 1.

Consider a continuous-time, linear, time-invariant SISO system described by

[TABLE]

We make a further assumption on the noise $e(t)$.

Assumption 4.1

$e(t)$* is white random process uncorrelated with $u$ having intensity $\lambda$.*

For an estimate $G(p,\theta)$ and $\tau$, the estimation error $\epsilon(t,\theta,\tau)$ can be computed as,

[TABLE]

As mentioned in Section 3, the delay can be estimated by minimizing the cost function $J(\theta,\tau)$, where $J(\theta,\tau)=\int_{-\infty}^{\infty}\!\epsilon(t,\theta,\tau)^{2}\,\mathrm{dt}$. If the estimation error is filtered by the low-pass filter $L(p)$, then an estimate of $\theta$ and $\tau$ can be computed by,

[TABLE]

where

[TABLE]

which by Parseval’s theorem is equivalent to,

[TABLE]

with $\Psi_{u}(\omega)$ the spectral density of $u(t)$ and the spectral density $\Psi_{v}(\omega)=\lambda$ follows from Assumption 4.1. If the transfer function $G_{0}(p)$ is known; the input signal white noise, i.e. $\Psi_{u}(\omega)=1$; and $\delta\tau=\tau-\tau_{0}$, the cost function $\bar{J}(\theta,\tau)$ can be written as,

[TABLE]

Recall that we are concerned with how to choose a set of filters and the cost function $J_{0}$ such that the two conditions of an M-safely redundant problem are satisfied. To satisfy the first condition we need to ensure all the cost functions share the same global minimum. The second condition requires the filter set to be constructed such that for any initial value of $\delta\tau$, there always exists a cost function whose solver path corresponds to a decrease in the cost function $J_{0}$ after a finite number of iterations. Note that these conditions can be checked if we know the locations of the global minimum and the extrema, i.e. the minima and maxima of $\bar{J}(\delta\tau)$.

First we establish a result for the location of the global minimum of the (non)filtered delay cost function.

Theorem 4.1

Consider the system $G_{0}(p)$ as described in (30). When $\lambda=0$, $\forall\delta\tau\in\mathbb{R}$, $\forall$ low-pass filters $L(p)\neq 0$, such that $\bar{J}(\delta\tau)\geq\bar{J}(0)$, the equality $\bar{J}(\delta\tau)=\bar{J}(0)$ occurs if and only if $\delta\tau=0$.

Proof. The proof of Theorem 4.1 is provided in Appendix A.1.

Theorem 4.2 provides a result for the locations of the extrema of the filtered cost function $\bar{J}(\delta\tau)$.

Theorem 4.2

If $L(p)$ is selected such that $L(p)G_{0}(p)$ is an ideal low-pass filter with cut-off frequency $\omega_{c}$ rad/sec, then the locations of the $i^{th}$ positive extrema of the time delay cost function $\bar{J}(\delta\tau)$ can be approximated by,

[TABLE]

where $T_{c}=2\pi/\omega_{c}$ and $i\geq 1$. When $i$ is even, the corresponding extremum is a minimum and when $i$ is odd, it is a maximum.

Proof. The proof of Theorem 4.2 is provided in Appendix A.2.

Remark 4.1

From Theorem 4.2, the locations of the extrema for the filtered cost function are known, hence the global convergence region is,

[TABLE]

corresponding to the distance from the maxima closest to the global minimum.

The result of Remark 4.1 can be seen in Fig. 1, which provides the plot of $\bar{J}(\delta\tau)$ when $\omega_{c}=2\pi$ (the constant $\lambda$ of the integral in (35) is set to 0 here for simplicity).

4.3 Time delay estimation with the useful redundancy technique

Now we can specify a filter set that satisfies the two conditions in Definition 1.

Condition 1: From Theorem 4.1, the first condition is satisfied for any set of low-pass filters as the global minimum of all the filtered time delay cost functions occurs when $\delta\tau=0$.

Condition 2: The second condition depends critically on the choice of filters and cost function $J_{0}$. The filter set needs to be chosen such that for any initial value of $\delta\tau$, there always exists a cost function whose solver path corresponds to a decrease in the cost function $J_{0}$ after a finite number of iterations.

To verify that condition 2 is satisfied, consider the necessary and sufficient conditions. To satisfy the necessary condition, we need to show that for any initial value $\delta\tau_{0}$, there always exists a path that has a minima location, $\delta\tau_{1}$, closer to the global minimum w.r.t. $\delta\tau_{0}$. For the sufficient condition, we need to establish that $J_{0}(\delta\tau_{1})<J_{0}(\delta\tau_{0})$.

First, we check the necessary condition. A simple choice is a filter set $L_{k}(p),k=\overline{1,n_{f}}$ where the cut-off frequencies, $\omega_{c,k}$, of $L_{k}(p)G_{0}(p)$ span from $1/\beta$ to 1 of the system bandwidth ($10\geq\beta>1$)555From Remark 3.1, $\beta$ is chosen $\leq 10$.. The cut-off frequencies are chosen based on the corresponding periods, $T_{c,k}$, being linearly spaced, i.e.,

[TABLE]

where $T_{bw}=2\pi/bw$, $bw$ is the bandwidth of the system $G_{0}(p)$ (rad/sec). This is due to the local minima of the filtered cost function being linearly dependant on the period of the cut-off frequency (Theorem 4.2).

Finally we need to prove the existence of such a filter set that satisfies the necessary condition of Condition 2 in Definition 1.

Theorem 4.3

There exists a filter set $L_{k}(p),\ k=\overline{1,n_{f}}$ chosen such that, $\forall\ \beta\geq\dfrac{5/4-1/(5\pi^{2})}{3/4-1/(3\pi^{2})}$, where $L_{k}(p)G_{0}(p)$ are ideal low-pass filters chosen with linearly spaced periods such that the cut-off frequencies, $\omega_{c,k}$, spanning from $1/\beta$ to 1 of the system bandwidth, i.e.

[TABLE]

where $T_{bw}=2\pi/bw$ and $T_{c,k}=2\pi/\omega_{c,k}$. Then,

[TABLE]

where $\xi(L_{q}(p),\delta\tau_{0})$ is the corresponding minimum found using the filtered time delay cost function generated by $L_{q}(p)$ with the initial delay $\delta\tau_{0}$.

Proof. The proof of Theorem 4.3 is provided in Appendix A.3.

Remark 4.2

From Theorem 4.3, we see that with the choice of,

[TABLE]

there always exists an $n_{f}$, that can be computed by (56) in Appendix A.3, that ensures the filter set defined in Theorem 4.3 satisfies the necessary part of Definition 1, Condition 2.

Remark 4.3

Note that (56) provides a loose lower bound. It is possible to have a smaller value of $n_{f}$ and still obtain a filter set that satisfies the necessary condition.

Remark 4.4

If $G_{0}(p)$ does not have resonant peaks then it doesn’t matter if $G_{0}(p)$ is known or unknown, as any ideal low pass filter $L_{k}(p)$ with bandwidth smaller than the bandwidth of $G_{0}(p)$ will allow $L_{k}(p)G_{0}(p)$ to approximate the desired low pass filter behaviour. If $G_{0}(p)$ has resonant peaks and the bandwidth of $L_{k}(p)$ is chosen significantly smaller than that of $G_{0}(p)$ then $L_{k}(p)G_{0}(p)$ will approximate the desired low pass filter behaviour.

An empirical method can also be utilized to determine the number of filters. For example, to check if a set of $n_{f}$ filters with the ratio $\beta$ satisfies the necessary condition, we can plot $n_{f}$ cost functions that follow from (50) in Appendix A.2 and check the minima locations using (36) to see that for any initial delay, $\delta\tau_{0}$, there always exists a cost function where the corresponding minima $\delta\tau_{i}$ is found closer to the global minimum with respect to $\delta\tau_{0}$. Fig. 2 shows a set of six cost functions with the ratio $\beta=2$ for a system bandwidth of $2\pi$ rad/sec, i.e. the cut-off frequencies span from $1/2$ to $1$ (Hz). From the figure, we can see that for any initial value of the time delay less than 15 sec, there always exists a cost function where the minima is found closer to the global minimum with respect to the initial delay. This can be confirmed by computing the minima of all the cost functions using (36) and observing the path to the global minimum.

Note that in practice, as $G_{0}(p)$ is not known a priori, and with the difficulty of implementing an ideal low-pass filter with irregular data, the filters $L_{k}(p)$ are implemented as Butterworth filters with cut-off frequencies $\omega_{c,k}$. Due to this, the extrema locations are not exactly those specified in (36). We discuss the effect of this difference by first citing the following theorem from [14], as a lemma.

Lemma 4.1

[14]** Consider a system with transfer function $F(s)=L(s)G_{0}(s)e^{-s\tau_{0}}$, let $u$ be the input signal and $w$ the corresponding output signal of the system. Then we have the following relationship,

[TABLE]

where $\gamma_{w}(t)$ is the autocorrelation function of $w(t)$ and $\rho_{w}(\delta\tau)=\gamma_{w}(\delta\tau)/\gamma_{w}(0)$.

Proof. See proof in [14].

With respect to the model setting in (30), $\gamma_{w}(\delta\tau)$ in Lemma 4.1 is the autocorrelation of the output of the system. Hence from Lemma 4.1, for an input signal of Gaussian distributed white noise, the filtered delay cost function $\bar{J}(\delta\tau)$ will have the form,

[TABLE]

where $b_{k}$ is the absolute value of the $k^{th}$ real pole of $L(s)G_{0}(s)$; $c_{k}$ and $d_{k}$ are the absolute values of the real and imaginary parts of the $k^{th}$ complex pole of $L(s)G_{0}(s)$; and weighting coefficients $h_{k}$ are related to the magnitude of the frequency response at the corresponding frequency. Augmenting a low-pass filter $L(s)$ to $G_{0}(s)$ obviously introduces slow poles into the system and by choosing $L(s)$ to be a Butterworth filter, it inherently introduces complex poles to the system. Hence, from (43), the filtered delay cost function will exhibit a slightly underdamped response. It can then be seen that the locations of the extrema of the cost function are associated with the underdamped response of the Butterworth filter. Hence, the augmentation of $L(s)$ with $G_{0}(s)$ when $L(s)$ is implemented as a Butterworth filter provides a cost function similar in nature to Fig. 1.

Next, to establish the sufficient condition, we need to derive a cost function $J_{0}$ that shares the same global minimum with the filtered cost functions where the minimum path corresponds to a decrease in $J_{0}$.

To analyse this condition, we take into consideration the model error that is generated in the iterations when estimating the rational part of the system due to an inaccurate estimate of the delay. Consider we have an estimate for the delay, $\hat{\tau}$, then when estimating the rational component of (30), we have the following relationship in the Laplace domain,

[TABLE]

where $\hat{G}(s)$ is the estimate for the rational part of the system. It is obvious that when $\hat{\tau}\neq\tau_{0}$, $\hat{G}(s)\neq G_{0}(s)$. Hence there always exists model error in estimating the rational part of the system if $\hat{\tau}\neq\tau_{0}$.

We now quantify the model error to understand the impact of it on the cost function. Here, a model reduction technique based on an identification method [19] is used to find a rational transfer function $G_{r}(s)$. The idea is to generate a noise-free data set from the system $G_{0}(s)e^{-(\tau_{0}-\hat{\tau})s}$, then use the SRIVC method to obtain a rational system $\hat{G}_{r}(s)$ that has the same model order as the true system and satisfies (44). For example, Fig. 3 shows the graphs of ten filtered delay cost functions after including the model error for the system $G(s)=2e^{-8s}/(0.25s^{2}+0.7s+1)$.

We now define the cost function $J_{0}$ as,

[TABLE]

where $\bar{\bar{J}}_{k}(\delta\tau)=100\dfrac{\parallel y(t)-\hat{\bar{x}}_{k}(t)\parallel_{2}}{\parallel y(t)-\mathbb{E}\{y(t)\}\parallel_{2}}$, $\hat{\bar{x}}_{k}(t)$ is the estimated system output when the filter $L_{k}(s)$ is used. Note that $\bar{\bar{J}}_{k}(\delta\tau)$ is the normalized version of the filtered cost function, $\bar{J}_{k}(\delta\tau)$.

The averaging is used to obtain a smooth decreasing cost function $J_{0}$ as there are always errors in the estimation due to noise and model error. Fig. 4 shows the average normalized mean square error for the ten filtered delay cost functions shown in Fig. 3. It can be seen that $J_{0}$ has the same global minimum as the ten filtered cost functions. Also shown in Fig. 4 is the minimum path through the filtered cost function that corresponds to a decrease in $J_{0}$. Note that the ‘$\ast$’ in Fig. 4 corresponds to a switch between cost functions $J_{k}(\delta\tau)$.

Therefore, in summary, the M-safely redundant optimization problem for the time delay estimation is defined as,

[TABLE]

with $J_{0}(\tau)$ defined as in (45). The filter set consists of $n_{f}$ Butterworth filters with cut-off frequencies, chosen such that the corresponding periods are linearly spaced, that span from $1/\beta$ to 1 of the system bandwidth where $\beta\geq\dfrac{5/4-1/(5\pi^{2})}{3/4-1/(3\pi^{2})}$ and $n_{f}$ is chosen using the methodology discussed in Section 4.3.

4.4 Algorithm

The proposed algorithm utilizing useful redundancy with multiple filters is described in Algorithm IV.

ALGORITHM IV

Step 1. Initialization

Select a set of low-pass filters as suggested in Section 4.3.

Set the boundaries $\tau_{min},\tau_{max},\Delta\tau_{min},\Delta\tau_{max}$666Note that $\Delta\tau_{max}$ is used to constrain the increment $\Delta\hat{\tau}$ in the time delay estimation, i.e. when $\Delta\hat{\tau}^{j}>\Delta\tau_{max}$, set $\Delta\hat{\tau}^{j}=\Delta\tau_{max}$. This is to ensure the Wolfe conditions are satisfied at each iteration in the Gauss-Newton technique [26]., and the SVF cut-off frequency $\omega_{c}^{SVF}$.

Choose an initial value $\hat{\tau}_{0}$ of the time delay.

Step 2. Iterative estimation

for i=1:converge

for k=1:$n_{f}$

(a)

Choose the low-pass filter $L_{k}(p)$.

(b)

Set the initial delay $\hat{\tau}^{0}_{i,k}$ as $\hat{\tau}_{i-1}$

(c)

Use the SRIVC algorithm and Gauss-Newton method as described in Algorithm II (without Step 3) to estimate the time delay $\hat{\tau}_{i,k}$.

(d)

Compute $J_{0}(\hat{\tau}_{i,k})$ (follows (45)).

end

Choose $\hat{\tau}_{i}=\underset{\hat{\tau}_{i,k}}{\text{argmin}}\ \ J_{0}(\hat{\tau}_{i,k})$.

If $|(\hat{\tau}_{i}-\hat{\tau}_{i-1})/\hat{\tau}_{i}|\geq\epsilon$, go to Step 1, else break. 777$\epsilon$ is a small value selected to obtain the desired accuracy.

end

Step 3. Refine parameter estimation

Repeat Step 2 with the low-pass filter $L(p)=1$ .

5 Numerical Examples

In this section, we demonstrate the effectiveness of the algorithm through numerical examples. These examples are commonly used in the literature [12] [16] [24] and act as defacto benchmarks to compare the performance between different CT system identification algorithms.

5.1 Case 1

Case 1 considers the second order system used in [12], i.e.

[TABLE]

In this example, we consider two sampling schemes for the input and output data.

Regular (uniform) sampling: The excitation input signal is a PRBS (pseudo-random binary sequence) of maximum length. The sampling time is 50ms, the number of stages in the shift register is 10 and the number of samples, N = 4000. 2. 2.

Irregular (nonuniform) sampling: The input excitation signal is a PRBS of maximum length. The number of stages in the shift register is 10 and the clock period is 0.5s. The input and output data are sampled at an irregular time instant $t_{k}$, with a sampling interval $h_{k}$ uniformly distributed as, $h_{k}\sim U[0.01,0.09](s)$. The realization of $h_{k}$ is randomized for each run. The number of samples, $N=4000$.

The additive output disturbance is Gaussian distributed white noise with zero mean designed to give Signal to Noise Ratios (SNR) of 5dB and 15dB. The SNR is defined as,

[TABLE]

with $P_{x}$ and $P_{e}$ the average power of the noise-free output $x(t_{k})$ and the disturbance noise $e(t_{k})$ respectively.

The system order is assumed known a priori (order 2). The frequency $\omega_{c}^{SVF}$ is chosen as 1 (rad/sec) which is same as the value used in [12]. To evaluate the algorithm, the time delay is initialized from different values $\hat{\tau}^{0}=0,1,3,5,7,9$ (similar to the experiment in [12]) and also from a random value drawn from a uniform distribution $U[0,9](s)$.

For the two sampling schemes, 100 different data sets are generated for each noise level. The system parameters and time delay are estimated using both the algorithm described in [12] and the algorithm proposed in this paper:

• •

First we consider the algorithm developed in [12]: The cut off frequency of the filter is chosen as $1/10$ of the system bandwidth as suggested in [12]. For regular data, the command ‘tdsrivc’ from the CONTSID toolbox [22] is used to estimate the time delay and the system parameters.

• •

For the multiple filtering algorithm utilizing the useful redundancy technique as proposed in this paper, the number of filters is set to 10 and the cut off frequencies, chosen based on linearly spaced periods, of the filters span from $1/10$ to $1$ of the system bandwidth. The order of the Butterworth filter is set to 10.

For both algorithms, the time delay boundaries $\tau_{min},\tau_{max}$ are set to [math] and $15$ seconds respectively. The maximum number of iterations for the Gauss-Newton method (Step 2 in Algorithm II and IV) is set to 10. The threshold, $\epsilon$, to determine convergence is set to $10^{-3}$.

We consider the estimated time delay to be the global minimum when the relative error, $\epsilon_{r}$, is less than $1\%$, i.e.

[TABLE]

where $\hat{\tau}$ and $\tau_{0}$ are the estimated delay and the true system delay respectively.

Fig. 5 contains plots of the 10 cost functions corresponding to the filters used in the useful redundancy method. It can be seen that they all share the same global minimum, however they all possess different local minima.

Fig. 6 presents a plot of $J_{0}$, for a single dataset, showing the trajectory of the time delay estimate when the sampled data is regular and the SNR is 15dB. The algorithm is started with the initial value of the time delay, $\hat{\tau}^{0}=0$. The ‘*’ in Fig. 6 represents the algorithm switching to another solver path. The algorithm stops when the condition in Algorithm IV is satisfied for the threshold $\epsilon=10^{-3}$, i.e. when the relative error in the time delay between consecutive iterations is $10^{-3}$.

Next we demonstrate the effectiveness of the proposed algorithm by showing the global convergence of 100 datasets with different initial values for the delay. The results are provided in Tables 1 and 2 where they are also compared to results from the existing method [12]. Note that the results presented here for the existing method are different to those achieved in [12] as we chose a larger bound on the delay. In [12], $\tau_{max}$ was set to 10 seconds while here, it is set to 15 seconds.

The results presented in Tables 1 and 2 show that the proposed method, which utilizes the useful redundancy technique, performs much better than the existing method [12] irrespective of the initial delay. When the initial delays are close to the true system delay, the global convergence percentage is high, i.e. approximately 100% for both algorithms. However, when the initial delays are poorly selected, the global convergence percentage is low for the existing method, e.g. with the initial delay $=0$, the percentage convergence for the existing method [12] is less than 1% while the proposed method achieves 100% global percentage convergence for both SNRs and the two sampling schemes.

5.2 Case 2: Rao-Garnier test system

In this section, we consider a system based on the Rao-Garnier continuous time benchmark [16],

[TABLE]

This system is linear, non-minimum phase with complex poles and a time delay. As with case 1, the experiment is conducted using both regular and irregular sampling schemes. For each sampling scheme, the input excitation signal is a PRBS of maximum-length. In the regular sampling case, the sample time is 10ms, the number of stages in the shift register is 10, the number of samples, N = 8000. For the irregular sampling, the number of stages in the shift register is 10 and the clock period is 0.5s. The input and output data are sampled at an irregular time instant $t_{k}$, where the sampling interval $h_{k}$ is uniformly distributed as, $h_{k}\sim U[0.01,0.05](s)$. The number of samples, $N=4500$.

We consider two SNR levels: 5dB and 15dB. The frequency $\omega_{c}^{SVF}$ is chosen as 25 rad/sec, which is approximately the bandwidth of the system. Again, the initial values for the time delay, $\hat{\tau}^{0}$, are selected as $0,1,3,5,7,9$ as well as a random value from the uniform distribution $U[0,9](s)$. The delay boundaries $\tau_{min},\tau_{max}$ are set to [math] and $15$ sec respectively.

For the existing method, the cut off frequency of the filter is chosen as $1/10$ of the system bandwidth. For the proposed method, the cut off frequencies of the filters, chosen based on linearly spaced periods, span from $1/10$ to $1$ of the system bandwidth and the number of filters used is 15. The order of the Butterworth filters is set to 10. The maximum number of iterations for the Gauss-Newton method (Step 2 in Algorithm II and IV) is 10. The threshold, $\epsilon$, for convergence is $10^{-3}$.

For comparison, 100 different data sets are generated for each noise level of each sampling scheme.

A graph showing the trajectory of a time delay estimate is plotted in Fig. 7 for a single data set. Switching between cost functions is indicated by ’$*$’. Note that, if only one filter were to be used, there is a high chance of convergence to a local minima. However, by using the useful redundancy technique, the proposed algorithm can avoid the local minima and converge to the global minimum.

The numerical results of the experiments are provided in Table 3 and 4. We can see that for both sampling schemes, the proposed method based on the useful redundancy technique performs much better as compared to the existing method. For all the initial values of delay used in this experiment, the existing method [12] never convergences 100% to the global minimum. However the proposed method utilizing useful redundancy still achieves a very high global convergence percentage, i.e. mostly 100% with any initial delay for both SNRs and sampling schemes.

6 Conclusion

The paper presents a new algorithm to estimate the parameters and time delay of a continuous-time system from regularly and irregularly sampled data. The idea is based on Instrumental Variable methods and employing the useful redundancy technique to enhance the global convergence by generating multiple cost functions by filtering the data several times. The paper also develops theoretical results related to the minima locations of the filtered delay cost function and the choice of filters to ensure the algorithm can converge to a global minimum. Numerical results show a significant improvement in the global convergence rate of the time delay estimation as compared to existing methods irrespective of the SNR.

Appendix A Appendix

A.1 Proof of Theorem 1

Proof A.1

From (35), when $\lambda=0$, $\bar{J}(\delta\tau)$ becomes,

[TABLE]

$\forall\delta\tau\in\mathbb{R},\ \omega\geq 0$, we have

[TABLE]

Therefore, $\bar{J}(\delta\tau)\geq 0,\ \forall\delta\tau\in\mathbb{R},\ \omega\geq 0$. Note that $\bar{J}(0)=0$, hence we have $\bar{J}(\delta\tau)\geq\bar{J}(0),\ \forall\delta\tau\in\mathbb{R},\ \omega\geq 0$. The equality occurs when $L(j\omega)=0,\ \forall\omega$ or $G_{0}(j\omega)=0,\ \forall\omega$ or $\cos(\omega\delta\tau))=1,\ \forall\omega$. As the system $G_{0}(p)$ and the filter $L(p)$ are not 0 (condition of Theorem 1), the equality only occurs when $\cos(\omega\delta\tau))=1,\ \forall\omega$, hence $\delta\tau=0$. $\square$

A.2 Proof of Theorem 2

Proof A.2

If $L(s)$ is selected such that $L(s)G_{0}(s)$ is an ideal low-pass filter with cut-off frequency $\omega_{c}$, then $\bar{J}(\delta\tau)$ becomes,

[TABLE]

where C=\dfrac{1}{2\pi}\int_{-\infty}^{\infty}\!\lambda\big{|}L(j\omega)\big{|}^{2}\mathrm{d}\omega.

The minima and maxima of $\bar{J}(\delta\tau)$ occur at the roots of,

[TABLE]

From (51), we can see that the extrema of $\bar{J}(\delta\tau)$ are also the extrema of the function $sinc(\omega_{c}\delta\tau)$. However, the minima of $\bar{J}(\delta\tau)$ will be the maxima of $sinc(\omega_{c}\delta\tau)$ and the maxima of $\bar{J}(\delta\tau)$ will be the minima of $sinc(\omega_{c}\delta\tau)$. Note that for the function $sinc(\omega_{c}\delta\tau)$, the locations of the $i^{th}$ ($i\geq 1$) positive extremum $\tilde{\delta}\tau$ can be approximated by,

[TABLE]

where $T_{c}$ is the corresponding period of the cut-off frequency $\omega_{c}$. The minima occurs when $i$ is odd and the maxima occurs when $i$ is even.

Therefore, for $\bar{J}(\delta\tau)$, the positive extrema can also be approximated using (52), the only difference is that the minima of $\bar{J}(\delta\tau)$ occurs when $i$ is even and the maxima occurs when $i$ is odd. $\square$

A.3 Proof of Theorem 3

Proof A.3

Here we use the positive extrema locations of the time delay cost function in (36). Since the filtered time delay cost function $\bar{J}(\delta\tau)$ is an even function, we only need to prove the theorem for $\delta\tau_{0}>0$.

Denote $\tilde{\delta}\tau_{i,k}^{min}$ as the $i^{th}$ positive minimum of the filtered time delay cost function $\bar{J}(\delta\tau)$ generated by the filter $L_{k}(p)$; $\tilde{\delta}\tau_{i,k}^{max}$ as the $i^{th}$ positive maximum of the filtered time delay cost function $\bar{J}(\delta\tau)$ generated by the filter $L_{k}(p)$. From Theorem 4.2, they can be computed by,

[TABLE]

From (53), for any filter $L_{k}(p)$, we have the following property,

[TABLE]

*Consider the two following cases:

Case 1: $\delta\tau_{0}\leq\tilde{\delta}\tau_{1,1}^{max}$, then from (54), we have, $\xi(L_{1}(p),\delta\tau_{0})=0$, which is smaller than $\delta\tau_{0}$ as $\delta\tau_{0}\neq 0$. Therefore, $\xi(L_{1}(p),\delta\tau_{0})<\delta\tau_{0}$. $\square$

Case 2: $\delta\tau_{0}>\tilde{\delta}\tau_{1,1}^{max}$. Denote $i_{0}$ such that,*

[TABLE]

Now consider a filter set $L_{k}(p),k=\overline{1,n_{f}}$ defined as in Theorem 4.3, with,

[TABLE]

where

[TABLE]

We need to prove the filter set satisfies the condition in (40). First, we prove that the filter set $L_{k}(p),k=\overline{1,n_{f}}$ has the following property,

[TABLE]

*where $i=\overline{1,i_{0}-1},k=\overline{1,n_{f}-1}$.

From (56),*

[TABLE]

Note that, as $\beta>1$,

[TABLE]

From (58) and (59),

[TABLE]

where $k=\overline{1,n_{f}-1}$, hence, combining (57) and (60),

[TABLE]

*where $k=\overline{1,n_{f}-1}.$

Next, consider the function*

[TABLE]

Taking the derivative w.r.t. $x$ of (62),

[TABLE]

Now $\forall x\geq 1$,

[TABLE]

Hence, $f(x)$ is an increasing function for $x\geq 1$, so we have, $f(i_{0}-1)\geq f(i),i=\overline{1,i_{0}-1}$, or,

[TABLE]

From (61) and (63),

[TABLE]

where $i=\overline{1,i_{0}-1},k=\overline{1,n_{f}-1}$, which can be rewritten as,

[TABLE]

*for $i=\overline{1,i_{0}-1},k=\overline{1,n_{f}-1}$.

From (64) and (53), we conclude,*

[TABLE]

Now, we will prove that $\tilde{\delta}\tau_{i,1}^{max}\geq\tilde{\delta}\tau_{i,n_{f}}^{min}$, $i=\overline{1,i_{0}}$. Similar to the previous analysis, it’s trivial to show that,

[TABLE]

As $\beta\geq\dfrac{5-4/(5\pi^{2})}{3-4/(3\pi^{2})}$, therefore, for $\forall i\in\mathbb{Z}^{+}$,

[TABLE]

Then by multiplying with $T_{bw}$, we have,

[TABLE]

Lastly, we show that for any filter set $L_{k}(p),k=\overline{1,n_{f}}$ defined as in Theorem 4.3 with $n_{f}\geq(1/M+\beta-2)/(1/M-1),n_{f}\in\mathbb{N}$,

[TABLE]

Case 2-A:

[TABLE]

From (54) and (67), it is obvious that $\xi(L_{1}(p),\delta\tau_{0})=\tilde{\delta}\tau_{i_{0}-1,1}^{min}$, and $\tilde{\delta}\tau_{i_{0}-1,1}^{min}<\delta\tau_{0}$, hence $\xi(L_{1}(p),\delta\tau_{0})<\delta\tau_{0}$.

$\square$**

Case 2-B:

[TABLE]

From (53) we see that, $\forall i\in\mathbb{Z}^{+}$,

[TABLE]

hence,

[TABLE]

Now combining (66) and (68),

[TABLE]

we have,

[TABLE]

Considering (70) and (71), there always exists a value $q\geq 2$ such that $\tilde{\delta}\tau_{i_{0}-1,q}^{min}<\delta\tau_{0}<\tilde{\delta}\tau_{i_{0}-1,q-1}^{min}$. Following from (65), we have $\tilde{\delta}\tau_{i_{0}-1,q-1}^{min}\leq\tilde{\delta}\tau_{i_{0},q}^{max}$. Therefore,

[TABLE]

Combining (54) and (72),

[TABLE]

therefore, $\xi(L_{q}(p),\delta\tau_{0})<\delta\tau_{0}$. $\square$