Consistent identification of continuous-time systems under multisine input signal excitation

For many years, the Simplified Refined Instrumental Variable method for Continuous-time systems (SRIVC) has been widely used for identification. The intersample behaviour of the input plays an important role in this method, and it has been shown recently that the SRIVC estimator is not consistent if an incorrect assumption on the intersample behaviour is considered. In this paper, we present an extension of the SRIVC algorithm that is able to deal with input signals that cannot be interpolated exactly through hold reconstructions. The proposed estimator is generically consistent for any input reconstructed through zero or first-order-hold devices, and we show that it is generically consistent for continuous-time multisine inputs as well. The statistical performance of the proposed estimator is compared to the standard SRIVC estimator through extensive simulations.


Introduction
System identification involves using measured input and output data for building mathematical models that characterise a system's behaviour. Different approaches to system identification have been developed depending on whether a discrete-time (DT) or continuous-time (CT) model is needed. Continuous-time system identification has applications in many areas of science and engineering such as economics, biology, physics and control, with comprehensive literature written on the subject [21,9,29]. Although the system identification community has focused mainly in DT setups, as it has been investigated during a predominantly digital era, there are many reasons why CT system identification has had a renewed interest during the last decades [6]. For example, model coefficients are directly linked to physical parameters, and more parsimonious models can be obtained as knowledge of the relative degree of the CT system can be accommodated. Also, contrary to the DT system identification using the forward shift This paper was not presented at any IFAC meeting. Corresponding author: R. A. González. Email addresses: grodrigo@kth.se (Rodrigo A. González), crro@kth.se (Cristian R. Rojas), siqi.pan@uon.edu.au (Siqi Pan), james.welsh@newcastle.edu.au (James S. Welsh).
operator, irregular and fast sampling can be easily handled, since the associated parameters remain invariant with respect to the varying sampling period and the model poles do not become statistically ill-defined as the sampling period decays to zero.
One of the main difficulties in CT system identification is the treatment of time derivatives. Since the goal is to obtain an estimate of a CT system, knowledge of the derivatives of the input and output are, either explicitly or implicitly, required. However, these derivatives are not exactly computable when only sampled input-output data is obtained. To overcome this problem, many algorithms have been suggested (see, e.g., [24,21] and the references therein). One of the most popular algorithms is the Simplified Refined Instrumental Variable method for Continuous-time systems (SRIVC), which was first presented in [31]. This method has been suggested for general use due to its robustness and accuracy in practical applications [11]. Many further extensions of this method also exist in the literature, for example, to handle non-uniformly sampled data [14] or multi-input systems [8]. Extensions to output error (OE) and Box-Jenkins (BJ) models [3], unification of DT and CT transfer function estimation [30], and a comprehensive consistency analysis [20] have also been presented.
The SRIVC algorithm uses interpolation of the input and output data in order to compute filtered regressor and instrument vectors in an iterative estimation procedure. This reconstruction of the CT input and output signals is usually implemented through simple interpolation schemes like zero-order hold (ZOH) or first-order hold (FOH) devices, independently of the nature of the true signals [10]. For inputs that can be described exactly with these reconstruction schemes, the SRIVC estimator has recently been shown to be generically consistent [20]. However, when the intersample behaviour assumption on the model input does not match that of the system input, continuous-time estimation methods can deliver large estimation errors if the sampling period is large [23], and in particular, the SRIVC estimator is known to be generically inconsistent in this case. Important input signals for identification that cannot be described by holds are band limited signals such as multisines. These signals are advantageous due to their flexibility regarding power spectrum design, time domain averaging possibilities, simplification of the model validation step and finite sample estimation performance [22]. For these input signals, the complete CT input signal is known to the practitioner, but the SRIVC procedure only performs simple interpolations of the input, which impact its consistency regardless of the sampling period.
In summary, in this paper, • we present a refinement of the SRIVC method that is shown to yield generic consistency of the estimated model parameters for CT multisine input signal excitations; • we prove that, given knowledge of the CT multisine input signal and measured output samples, the exact computation of the input regressors is necessary and sufficient for a generically consistent estimate of the CT system; • we propose a computationally efficient algorithm for computing the regressors under the multisine case, and introduce an extension of the SRIVC algorithm for arbitrary CT inputs, which are not necessarily constructed through hold devices; and • we exemplify the consistency of the proposed estimator through extensive Monte Carlo simulations.
The remainder of this paper is organised as follows. The identification problem is formulated in Section 2. Section 3 provides a description of the SRIVC estimator and its consistency properties. The proposed SRIVC-type method is presented and analysed in Section 4, and Section 5 illustrates this method with extensive numerical examples. Finally, conclusions are drawn in Section 6.

Problem formulation
Consider a linear and time-invariant (LTI), causal, stable, proper, single-input single-output, CT system where p is the Heaviside operator, i.e., pg(t) := dg(t)/dt, and the numerator and denominator polynomials are coprime and given by Suppose that the CT input u(t) is known from t = t 1 to t = t N , and that N noisy measurements of the output x(t) are obtained at the instants {t k } N k=1 . In other words, the output observations are given by where it is assumed that the sampled noise sequence {v(t k )} can be described as a zero-mean and finite variance random process. Due to the nature of the sampled signals and the difficulty of computing the timederivative of CT white noise, which does not have finite variance [1], we only consider DT noise in this paper.
To identify the system, we propose the model structure where the parameter vector θ := a 1 , a 2 , . . . , a n , b 0 , b 1 , . . . , b m needs to be estimated. The goal is to obtain an accurate model of the CT system G * (p) := B * (p)/A * (p) given the knowledge of N samples of the output measurements and the CT input signal. Note that in this framework the input signal is not limited to hold reconstructions. Hence, the description includes the standard framework where u(t) is assumed to be obtained through a ZOH or FOH and extends to more general inputs, such as multisines and band limited signals [22].
The identification of the system G * (p) can be done by obtaining the data points {u(t k ), y(t k )} and applying a method for CT system identification, such as in [3], or as in [31,17,12] for regular sampling schemes. In most of these algorithms, however, the hold reconstructions of the input and output are assumed, and they are independent of the exact nature of the signals. In this work, we show that the knowledge of the exact intersample behaviour of the input can provide further insights for a better design of the identification procedure.
3 The Simplified Refined Instrumental Variable method for Continuous-time systems (SRIVC) The SRIVC estimator is an adaptive instrumental variable algorithm where parameter-dependent CT filters are updated iteratively. In each step, the instruments are computed using the parameter estimate obtained in the previous iteration until the model parameters have converged. The iterative procedure of the SRIVC algorithm is designed so that the sum of squares of the residuals (also called the generalised equation errors or GEEs) ε(t k ), is minimised. The residuals are written as where Note that in (2) and (3) we have adopted a mixed notation of CT operators and DT data. Since this dichotomy is repeatedly encountered in this paper, we formalise it in the following remark.
Remark 1. In this paper, G(p)x(t k ) means that the DT signal x(t k ) is interpolated in some manner, e.g., using a ZOH or FOH, and the resultant output through the CT filter G(p) is sampled at t = t k . On the other hand, and later sampled at t = t k .
The SRIVC method is described in Algorithm 1, where we denote ϕ f (t k ) as the filtered regressor vector,φ f (t k ) as the filtered instrument vector, and y f (t k ) as the filtered output. Note that line 8 of Algorithm 1 requires the DT signals to be prefiltered by CT transfer functions. This is usually done by assuming a ZOH or FOH reconstruction for the input and output signals and then simulating the response by using, for example, the lsim command in MATLAB. Although this approach has provided a quick procedure to compute the filtered regressor and instrument vectors, it is prone to approximation errors that can jeopardise the statistical properties of the method.

Remark 2.
In the SRIVC method, the user has several choices regarding the intersample behaviour assumptions. In particular, the intersample behaviour of the input in both (4) and (5) can be chosen, as well as the reconstruction of the output signal for the filtering steps in (4) and (6). Usually the output is selected to have a FOH behaviour, since it is argued that it typically gives rise to a satisfactory approximation if the sampling period is small [4].

Consistency Analysis of the SRIVC estimator
Previous works [28,30] have suggested that the SRIVC estimator uses the optimal instrumental variable terms, , y(t k ))} N k=1 , model order (n, m), initial vector estimate θ 1 , tolerance and maximum number of iterations M 2: Using θ 1 , form the estimated system polynomials A 1 (p) and B 1 (p) 3: j ← 1, flag ← 1 4: while flag = 1 and j ≤ M do 5: if Reflect the unstable poles of 1/A j (s) into the stable region of the complex s-plane 7: end if 8: 9: Compute the parameter estimate 10: if θ j+1 − θ j θ j < then j ← j + 1 14: end while 15: Output: θ j and its associated model B j (p)/A j (p). and that it minimises the prediction error and maximises the likelihood function, but they lack rigorous theoretical analysis regarding the influence of the interpolation of the input and output for the prefiltering step. Only recently [20] has the intersample behaviour of the signals been taken into account for the consistency analysis. In [20,Theorem 1], the generic consistency of the SRIVC estimator was proven for inputs that can be exactly interpolated by FOH or ZOH devices. More precisely, under mild assumptions regarding the sampling period and persistence of excitation of the input, the following statements are true for an input that is exactly reconstructible with FOH or ZOH interpolation: singular. 1 2. The true parameter θ * is the unique converging point. 3. As the sample size N approaches infinity, θ j+1 in (7) converges to θ * for j ≥ 1.
Also, the effect of choosing a different intersample behaviour than that of the system input was also analysed in [20,Corollary 3]. In the following, we say that a correct intersample behaviour in the input signal is assumed whenever the intersample behaviour of such signal in the SRIVC algorithm matches that of the system input.
In [20] it was shown that the SRIVC estimator 1. remains generically consistent if an incorrect assumption on the intersample behaviour is used for generating the filtered signals in the instrument vectorφ f (t k ), and 2. is generically inconsistent if an incorrect assumption on the intersample behaviour is used for filtering the input signal in the regressor vector ϕ f (t k ).
This result indicates that the intersample behaviour of the input signal needs to be correctly taken into account for the consistency of the SRIVC estimator. In particular, it implies that if the system input is a signal that is not produced by a hold mechanism, the estimator will be generically inconsistent. This argument holds regardless of whether the additive noise v(t k ) is white or coloured.

Consistent SRIVC-type method
As mentioned in the previous section, a correct assumption of the intersample behaviour of the input (ZOH of FOH) in the regressor vector ϕ(t k ) guarantees generic consistency under mild conditions. The extension of this principle constitutes our main contribution. In this work, we propose an extension of the SRIVC method that computes the filtered regressors exactly for a wide class of input signals whose intersample behaviour is known and prove its consistency for multisine input excitations.
The generalised equation error for the proposed approach is where u f (t) = 1 A(p) u(t). In (8), the predicted output measurement is explicitly calculated by first computing the underlying CT signal, and later evaluating it at t = t k . The proposed estimator follows the procedure described in Algorithm 1, but the filtered regressor and instrument vectors in Equations (4) and (5) now become and Note that the t k in (10) follows the notation in Remark 1.

Remark 3.
The proposed estimator is an extension of the standard SRIVC estimator that uses the complete CT input signal for identification. For input signals that are reconstructed exactly through a ZOH or FOH, this estimator is equivalent to the SRIVC estimator. Thus, the SRIVC-type estimator with prefiltering stage given by (9) and (10) is generically consistent under the same assumptions as in [20].
In order to further analyse the asymptotic properties of the proposed estimator, we first study its consistency for multisine inputs. Later, an extension for arbitrary input excitations is presented.

Multisine input signal
We consider multisine input signals of the form where m u , {α l } mu l=0 , {ω l } mu l=0 and {ψ l } mu l=0 are input parameters. The frequencies ω l are assumed to be positive and distinct, and without loss of generality we assume that the weights α l are positive as well. It is well known that the output in stationary regime of an asymptotically stable LTI filter H(s) when u(t) is applied is also a multisine, given by y(t) = H(0)α 0 + mu l=1 α l |H(iω l )| cos(ω l t + ψ l + ∠H(iω l )). (12) This property of LTI systems provides a natural way to obtain exact values for the signal evaluations in (9) and (10), and it is of low computational cost, since the prefiltering is directly obtained by evaluating (12) with the corresponding filter. Another advantage of this approach is that it extends naturally to non-uniformly sampled data. For such type of sampling, the proposed method is not as computationally intensive as the standard SRIVC method, since the algorithm only obtains approximations of the filtered output p i A −1 j (p)y(t k ), i = 0, . . . , n, instead of computing approximations of the filtered values of both u(t k ) and y(t k ). The filtered output computations can be carried out by, e.g., an adaptive Runge-Kutta method (as in [3]), or by any oversampling technique with intersample behaviour assumptions.
We now prove the consistency of the proposed estimator for the multisine input. The assumptions we use during the analysis are the following: and asymptotically stable with A * (p) and B * (p) being coprime. Since deterministic inputs will be considered in conjunction with stochastic noise processes, our analysis uses the standard definition of expectation for quasi-stationary signals [16, pp. 34], which is Theorem 4. Consider the SRIVC-type estimator with a fixed sampling period h and filtered regressor and instrument vectors given by (9) and (10) respectively, and suppose that Assumptions (A1) to (A5) hold. Then, the following statements are true: 1. There exists a maximum sampling period h * > 0 such that, if h ≤ h * , the matrix E{φ f (t k )ϕ f (t k )} is generically non-singular. 2. If h ≤ h * , the true parameter θ * is the unique converging point. 3. As the sample size N approaches infinity, θ j+1 converges to θ * for j ≥ 1.

Proof.
Proof of Statement 1. By substituting and On the other hand, we also havê and S(−B j , A j ) is the Sylvester matrix associated with the polynomials −B j (p) and A j (p), whose nonsingularity follows from the same analysis done in [20], where Assumption (A4) is used. With this, we compute Thus, for showing that E{φ f (t k )ϕ f (t k )} is generically non-singular for a small enough sampling period h, it is sufficient to show that Ψ = 0 and Φ is generically nonsingular for a small enough h. The difference between the analysis in [20, Theorem 1] and the proof in the current paper is that the signals of interest are hybrid in nature: some are evaluations of CT signals, whereas others are DT signals interpolated with a reconstruction device, such as a FOH.
The proof of Ψ = 0 can be found in Lemma 8 in the Appendix. Regarding the invertibility of Φ, we will conve- and with S(−B * , A * ) being the Sylvester matrix associated with the polynomials −B * (p) and A * (p), which is nonsingular since A * (p) and B * (p) are coprime. Hence, we can write the expected value of interest as where It is shown in Lemma 9 that Φ 1 is generically nonsingular, which means that the first summand of the right hand side of (18) is generically non-singular.
Finally, as h tends to zero, the infinity norm of the difference between the direct evaluation of a CT signal and its interpolated counterpart also tends to zero. Thus, ∆(t k ) → 0 as h → 0. This, together with the fact that (generic) non-singularity of a matrix is preserved under small-enough matrix perturbations [13,Chap. 6], leads to the first statement of the theorem.

Statement 2.
Suppose thatθ is a limiting point of the iteration in (7), where ϕ f (t k ) andφ f (t k ) are defined as in (9) and (10) respectively, and the corresponding polynomials are denoted byĀ(p) andB(p). These polynomials are coprime by Assumption (A4  (7), at the converging point and as N tends to infinity, as (19) where ε(t k ,θ) is the GEE (8) evaluated at the converging point. Since the matrix inverse in (19) is assumed to be non-singular, the second expectation in (19) must be zero, i.e., where r = max(n + m * , n * + m) = n + m. Then, the GEE in (8) can be rearranged as where H = h 0 , h 1 , . . . , h n+m . Now, note that the instrument vectorφ f (t k ) can be written aŝ where S(−B,Ā) is a Sylvester matrix associated with the polynomialsB(p) andĀ(p), which again is nonsingular. So, we can express (20) as Following a similar approach as in Lemma 8, we conclude thatΨ = 0, and by Lemma 9,Φ is generically nonsingular. Thus, for (21) to hold we need H = 0, which implies thatĀ i.e., θ * is the unique limiting point.

Statement 3. The proof follows from the analysis made for proving Statement 3 of Theorem 1 in [20]. 2
Note that if the commonly used FOH (or ZOH) were chosen as the intersample behaviour of the signals when discretising the prefilters, the reconstruction of u(t) would suffer from high frequency distortion, which usually leads to inaccuracies in the computation of ϕ f (t k ) andφ f (t k ). As stated next, only an inaccurate computation of the regressor vector ϕ f (t k ) causes generic inconsistency of the proposed method under CT multisine input excitation.
Corollary 5. Assume that an incorrect intersample behaviour in the input signal is assumed, but nevertheless satisfies G(p)u(t k ) = {G(p)u(t)} t k as h → 0. The SRIVC-type estimator with filtered regressor and instrument vectors given by (9) and (10) respectively 1. remains generically consistent if an incorrect assumption on the intersample behaviour is used for generating the filtered signals in the instrument vectorφ f (t k ), and 2. is generically inconsistent if an incorrect assumption on the intersample behaviour is used for filtering the input signal in the regressor vector ϕ f (t k ).
Proof. Statement 2 : Statement 1 of Theorem 4 still holds by following the same steps as before, but this time the vector ∆(t k ) in (17) will also have non-zero elements in its bottom m + 1 entries. Namely, the i-th component of ∆(t k ), with i = n + 1, . . . , n + m + 1, is now A similar procedure to (9) and (10) could be proposed for the computation of y f (t k ) by exploiting the fact that the noiseless output also corresponds to a multisine (thus, a more adequate reconstruction scheme could be designed). However, Remark 5 of [20] suggests that, as the number of iterations tends to infinity, the GEE at the converging point does not depend on the intersample behavior of the output. Thus, if the iterations converge, a more precise filtering of the output is not needed.

Extension to arbitrary input signals
The previous method is naturally suited for multisine inputs due to the simplicity of the filtered outputs at stationary regime, since they are also multisine signals. By introducing a Delta transform description, the computations in (9) and (10) can be generalised for an arbitrary input signal with arbitrary accuracy. For this procedure, only regular sampling is considered, although extensions to irregular sampling are also possible.
The proposed algorithm for computing ϕ f (t k ) and ϕ f (t k ) at the j-th iteration of the SRIVC-type method goes as follows: 1. Given the sampling period h of y(t), (over)sample u(t) with sampling period δh, where δ 1. 2. From θ j , form the prefilters of interest, namely p i /A j (p) and B j (p)p l /A 2 j (p) for i = 0, . . . , m; l = 1, . . . , n. 3. Compute the Delta equivalent [18] of the prefilters, and calculate the response at instants t k of each filter to the fast-sampled version of u(t) in the Delta domain.
The Delta domain description leads to an exact simulation of the underlying CT system when δh → 0. Thus, the SRIVC-type method with prefilters computed as in steps 1, 2 and 3 above calculates the filtered regressor and instrument vectors accurately if the oversampling period δh is chosen small as compared to the sampling period. Due to potentially high sampling rates, the use of the Delta operator is needed for ameliorating rounding errors and ill-conditioning problems regarding the sensitivity of the coefficients of the prefilters. Note that step 3 of the algorithm can be easily performed by using the Delta domain Toolbox in MATLAB.
Remark 7. The precision of this procedure will depend on the over-sampling period δh. Via extensive simulations, we have found that sampling at least 100 times faster is usually enough to provide reliable estimates.

Simulation examples
Via numerical simulations we compare the performance of the standard SRIVC method with the proposed SRIVC-type method, which is labeled SRIVC-c. For a multisine input, we examine the consistency of both methods for different regular sampling period and also for irregular sampling. We also study the consistency of each method under an arbitrary input excitation for different regular sampling periods. For the following tests, we consider the system G * (p) = 1.25 0.25p 2 + 0.7p + 1 , where the parameters of interest are a * 1 = 0.25, a * 2 = 0.7, and b * 0 = 1.25. Regarding the implementation of the standard SRIVC method, we have used the srivc command from the CONTSID toolbox version 7.3 for MAT-LAB [7], under default initialisation and tolerance settings. It was set to estimate the best model among the correct model structure with a FOH as the intersample behaviour.

Multisine input: Regular sampling
We first test if the algorithms provide consistent estimates of the parameter vector [a 1 , a 2 , b 0 ] . The system in (22) is excited with the CT input u(t) = sin(0.714t) + sin(1.428t) + sin(2.142t).
The noiseless output is computed analytically by assuming that it corresponds to the output of the system at the stationary regime, i.e.,  SRIVC-c estimator accurately identifies all parameters while the standard SRIVC method fails to recover the true parameter vector as N increases. The MSEs for the SRIVC-c estimator decrease to zero as expected for consistency, while at least two out of the three estimated parameters given by the SRIVC method are biased, which empirically indicates that the SRIVC estimator is not consistent in this example.

Multisine input: Different sampling periods
We now study the effect of the intersample behaviour on the SRIVC-type estimates. Under the same input and noise variance as the previous simulation, we test the performance of each algorithm for a fixed number of output measurements (N = 2000) with different regular sampling periods. Since the rise time of the system is approximately 2 seconds, a good choice for the sampling period should be between 0.2 and 0.5 seconds according to the criterion suggested in [2]. In order to cover fast, normal and slow sampling, we test with sampling periods h = 0.06, 0.2 and 0.6[s].
The sample mean and mean square error of each parameter over 300 Monte Carlo runs for each sampling period are shown in Table 1. On average, the SRIVC-c estimator delivers the true values of every parameter for all sampling periods in this study, whereas the SRIVC estimator only performs well when the sampling period is small. For h = 0.6[s], the large sampling period exaggerates the interpolation error of the input signal in the standard SRIVC estimator, which severely degrades its performance. This is confirmed by the order of magnitude of difference in MSE of the parameters given by the two estimators.

Multisine input: Irregular sampling
We consider the same system described before, with the same input and noise variance. In this simulation study, 2000 irregularly sampled output measurements are obtained. The sampling interval is distributed uniformly between h lb and h hb , where the lower bound is fixed at h lb = 0.05, while the upper bound is varied from 0.1 to 0.6. A total of 6 Monte Carlo simulations are performed with each simulation containing 300 runs. Figure 3 shows the mean value of each parameter, with their standard deviation around this value. As expected, the SRIVC-c estimator provides accurate estimates for all sampling period ranges in this study. On the other hand, the SRIVC estimator has a degrading performance as the sampling range increases, which could be attributed to the approximation errors in the prefilter calculations due to incorrect assumptions on the intersampling behaviour.

Arbitrary input application: Chirp signals
The next goal is to check whether the algorithms can provide accurate estimates for arbitrary input signals. Now, the system in (22) is excited with an up-chirp signal, which is a CT signal that increases in frequency with time. These signals are widely used in signal processing applications such as radar systems and seismology and have been used for system identification [27,19].
The chirp signal used in this example is where f 0 = 0.1[Hz] and f 1 = 0.6 [Hz], and T f = 500[s] is the length of one period of the chirp signal. In this case, we determine the true system output using the explicit Runge-Kutta formulae RK5(4) [5]. The output is sampled every h = 0.5[s] with δ = 0.001, and the measurement noise has variance 0.05, which corresponds to approximately 10% of the energy of the noiseless output. For the computation of the SRIVC-c estimate, we follow the algorithm described in Section 4.2. The number of periods of the input signal vary from 1 to 10, which is equivalent to sample sizes ranging from 1000 to 10000, and 300 Monte Carlo runs are performed for each case. The empirical evidence in Figure 4 suggests that the SRIVC-c estimator can also lead to accurate estimates for signals that are not described exactly by hold reconstructions nor multisines.

Conclusions
In this paper, we have derived an algorithm for continuous-time system identification that is consistent for a wide class of input signals that have a known intersample behaviour. This estimator extends the applicability of the standard SRIVC method to inputs that are not exactly described by hold devices. We put forward a comprehensive analysis of the generic consistency of this estimator for multisine inputs, and extensive simulations have confirmed the advantages of this estimator over the widely popular SRIVC method. Further research on this topic concerns variance analyses of these estimators, and theoretical guarantees for irregular sampling schemes. Using this result, and the identity cos(α) cos(β) = [cos(α + β) + cos(α − β)]/2, the second term in the sum in (24a) is zero. Moreover, in (24b) the term for j = l is a sum of sinusoids whose sum tends to zero as N tends to infinity, while for j = l constants appear. Thus, which leads to cos(φ j r − φ l r ) = Re{Ā(iω r )} |Ā(iω r )| cos π 2 (l − j) + ∠A * (iω r ) + Im{Ā(iω r )} |Ā(iω r )| sin π 2 (l − j) + ∠A * (iω r ) .