A Recursive System Identification With Non-Uniform Temporal Feedback Under Coprime Collaborative Sensing

Ouyang, Jinhua; Chen, Xu

doi:10.1115/1.4063481

Abstract

We present a system identification method based on recursive least-squares (RLS) and coprime collaborative sensing, which can recover system dynamics from non-uniform temporal data. Focusing on systems with fast input sampling and slow output sampling, we use a polynomial transformation to reparameterize the system model and create an auxiliary model that can be identified from the non-uniform data. We show the identifiability of the auxiliary model using a Diophantine equation approach. Numerical examples demonstrate successful system reconstruction and the ability to capture fast system response with limited temporal feedback.

1 Introduction

A common assumption for real-time control systems design is that the sampling of input and output signals is uniform, periodic, and synchronous [1]. In the information-rich world, however, data streams are often non-uniform and asynchronous. (In fact, real-time control system implementations often have to adjust the sampling rate to deal with irregular data [2].) While non-uniformly sampled data intuitively contain more temporal information for system analysis and controls [3,4], they violate the classical real-time control framework, and most existing methods for non-uniformly sampled systems are heuristic and specific [5]. It remains not well understood how to systematically leverage non-uniform data streams for real-time dynamic systems. In particular, as the first critical step in real-time controls, classic system identification requires synchronous input and output data when building the model of a dynamical system [6].

From a signal processing point of view, non-uniform data are naturally dense in certain temporal regions where more information about the system dynamics can be revealed [4]. The non-uniformly sampled data can be collected by triggering the sensor with events, by randomized sampling, or by fusing measurements from multiple sensors. On the one hand, the temporal resolution is increased due to the data irregularity [3]. On the other hand, it challenges conventional system identification algorithms.

One approach to identifying a system under non-uniform data is based on the approximation theory [7]. Briefly, the non-uniformly collected data are approximated or reconstructed by a sequence of uniform samples, and then, the conventional system identification algorithms can be applied to the resulting uniform data [8]. Several techniques have been proposed for the data reconstruction, including linear [9], polynomial [10], and spline interpolations [11]. Other works on system identification subject to non-uniformly sampled data have also been conducted using the expectation maximization approach [12–14], the maximum likelihood estimation [14,15], the lifting operator [16,17], and the output error method [18].

Stepping further beyond the existing approaches, this paper contributes to a novel system identification that leverages the temporal advantage of non-uniform sampling but overcomes the obstacle imposed by non-uniform data collection for general input–output models. We first propose a coprime collaborative sensing scheme, which generates one set of data that appears non-uniform with respect to time while, in the meantime, having systematic underlying sampling patterns. Next, we implement a model reparameterization tailored for the selected sensing scheme based on polynomial transformation to construct an auxiliary model that can be directly identified with the available observations. Then, a recursive least-squares (RLS)-based algorithm is designed to identify the auxiliary model and to illustrate the feasibility of working with the mechanism of collaborative sampling and model reparameterization. Lastly, the parameters of the original fast system model are recovered by removing the highest common factors between the denominator and numerator polynomials.

The remainder of this paper is organized as follows. In Sec. 2, technical preliminaries regarding the model reparameterization are reviewed and introduced. The proposed coprime collaborative sensing and system modeling are formally defined in Sec. 3. In Sec. 4, we derive recursive system identification algorithms based on the proposed sensing scheme and model reparameterization strategies. Section 5 contains multiple classes of numerical examples. Section 6 concludes the paper.

2 Preliminaries

In this section, we review the transfer operator and model parameterization in standard single-rate and multi-rate system identification. Consider the deterministic autoregressive-moving-average model for a linear time-invariant system:

y (k) = \frac{q^{- d} B (q^{- 1})}{A (q^{- 1})} u (k)

(1)

where d is an integer number of sampling periods contained in the time delay of the systems, and q is the time-domain shift operator defined as qy(k) = y(k + 1) and q⁻¹y(k) = y(k − 1), and

A (q^{- 1}) = 1 + a_{1} q^{- 1} + \dots + a_{n_{a}} q^{- n_{a}}

and

B (q^{- 1}) = b_{1} q^{- 1} + b_{2} q^{- 2} + \dots + b_{n_{b}} q^{- n_{b}}

are polynomials of q⁻¹. Equation can be rewritten as

A (q^{- 1}) y (k) = q^{- d} B (q^{- 1}) u (k)

(2)

or alternatively

(1 + \dots + a_{n_{a}} q^{- n_{a}}) y (k) = (b_{1} q^{- 1} + \dots + b_{n_{b}} q^{- n_{b}}) u (k - d)

⁠. Expanding and rearranging yield

\begin{aligned} y (k) & = - a_{1} y (k - 1) - \dots - a_{n_{a}} y (k - n_{a}) \\ + b_{1} u (k - d - 1) + \dots + b_{n_{b}} u (k - d - n_{b}) \end{aligned}

(3)

or y(k) = θ^Tϕ(k), where

\begin{aligned} ϕ (k) & = [- y (k - 1), - y (k - 2), \dots, - y (k - n_{a}), u (k - d - 1), \\ u (k - d - 2), \dots, u (k - d - n_{b})]^{T} \\ θ & = [a_{1}, a_{2}, \dots, a_{n_{a}}, b_{1}, b_{2}, \dots, b_{n_{b}}]^{T} \end{aligned}

When the system input and output are sampled at different rates (slower output sampling in this paper), the available data become y(Jk) = {y(k − J), y(k − 2J), …} and u(k) = {u(k − 1), u(k − 2), …}, where J is a positive integer representing the ratio between the input and output sampling rates. The original single-rate model described in Eq. can be transformed into a dual-rate version that can be identified directly from the available measurements [19]. The solution approach is first to recognize the factorization:

1 - x^{J} = (1 - x) (1 + x + x^{2} + \dots + x^{J - 1})

(4)

Next, consider the characteristic equation A(q⁻¹) in the multiplication form:

A (q^{- 1}) ≜ \prod_{i = 1}^{n_{a}} [1 - (λ_{i} q)^{- 1}]

(5)

where λ_i’s are the reciprocals of the poles of the system, and n_a is the order of the characteristic equation (i.e., the number of poles). Observing the structure of Eq. , we notice that by designing a polynomial:

\begin{aligned} F_{J} (q^{- 1}) & = \prod_{i = 1}^{n_{a}} [1 + (λ_{i} q)^{- 1} + (λ_{i} q)^{- 2} + \dots + (λ_{i} q)^{- J + 1}] \\ = 1 + f_{1} q^{- 1} + \dots + f_{n_{a} J - n_{a}} q^{- n_{a} J + n_{a}} \end{aligned}

(6)

the original characteristic equation described in Eq. can be transferred into

\begin{aligned} A_{J} (q^{- J}) & = A (q^{- 1}) F_{J} (q^{- 1}) \\ = [1 - (λ_{1} q)^{- J}] [1 - (λ_{2} q)^{- J}] \dots [1 - (λ_{n_{a}} q)^{- J}] \\ = 1 + a_{J, 1} q^{- J} + a_{J, 2} q^{- 2 J} + \dots + a_{J, n_{a}} q^{- n_{a} J} \end{aligned}

with a down-sampled observation space. Applying the same transformation polynomial shown in Eq. to the numerator of Eq. yields a multi-rate system model:

\begin{aligned} y (k) & = \frac{q^{- d} B (q^{- 1}) F_{J} (q^{- 1})}{A (q^{- 1}) F_{J} (q^{- 1})} u (k) \\ = \frac{q^{- d} B_{J} (q^{- 1})}{A_{J} (q^{- J})} u (k) \end{aligned}

(7)

or in a form similar to Eq. :

A_{J} (q^{- J}) y (k) = q^{- d} B_{J} (q^{- 1}) u (k)

(8)

where

B_{J} (q^{- 1}) = b_{J, 1} q^{- 1} + \dots + b_{J, n_{a} (J - 1) + n_{b}} q^{- n_{a} (J - 1) - n_{b}}

⁠. Rewriting Eq. in a predictor form similar to Eq. , we have

\begin{aligned} y (k) & = - a_{J, 1} y (k - J) - \dots - a_{J, n_{a}} y (k - n_{a} J) \\ + b_{J, 1} u (k - d - 1) + \dots + b_{J, n_{b}} u (k - d - n_{a} (J - 1) - n_{b}) \end{aligned}

where the output prediction is precisely a linear combination of u(k) and y(Jk). Therefore, directly identifying the multi-rate model parameters becomes possible after the aforementioned model reparameterization.

The key insight of the introduced model reparameterization is to recognize that the historical measurements required for system identification depend solely on the order of system polynomials (i.e., A(q⁻¹) and B(q⁻¹)). By designing a transformation polynomial, we can freely adjust the order of system polynomials. Consequently, the challenge posed by input and output asynchronism in identifying system dynamics is effectively circumvented.

3 Proposed Coprime Collaborative Sensing and Model Reparameterization

Figure 1 illustrates the proposed coprime collaborative sensing scheme, where multiple sensors with coprime sampling rates collaboratively sense the system output. Assuming the fundamental sampling period is T, and S represents the set of sensors sampling rate, we define the coprime sampling rate as S = {aT, bT, cT, …}, where a, b, c, … are coprime integers. The data collected from these sensors are then combined chronologically, assuming that all sensors begin sampling simultaneously. The coprime sampling rates result in fewer measurements overlapping when multiple sensor measurements are fused, providing the highest temporal resolution as more details of the system response become available. This enables the parameter estimation to be updated with the maximum information entropy precisely when all sensor measurements overlap.

Fig. 1

View large Download slide

The proposed collaborative sensing scheme of multiple sensors with coprime sampling rates. The illustration depicts the case when three coprime sensors’ data are merged for use (boxed in dashed lines). The instants enclosed by solid lines represent valid measurements for updating parameter estimation (i.e., when all sensors’ measurements overlap).

When n sensors of different rates are used, the available output measurements become

\begin{aligned} y (J_{1} k, \dots, J_{n} k) & = {y (k - J_{1}), y (k - 2 J_{1}), \dots \\ ⋮ \\ y (k - J_{n}), y (k - 2 J_{n}), \dots} \end{aligned}

Based on the aforementioned multi-rate model reparameterization, we know that there will be n unique transformation polynomials

\begin{aligned} F_{J_{1}} (q^{- 1}) & = \prod_{i = 1}^{n_{a}} [1 + (λ_{i} q)^{- 1} + (λ_{i} q)^{- 2} + \dots + (λ_{i} q)^{- J_{1} + 1}] \\ ⋮ \\ F_{J_{n}} (q^{- 1}) & = \prod_{i = 1}^{n_{a}} [1 + (λ_{i} q)^{- 1} + (λ_{i} q)^{- 2} + \dots + (λ_{i} q)^{- J_{n} + 1}] \end{aligned}

if we want to identify the system model with the given output observation space. Let

F_{J_{*}} (q^{- 1}) = \frac{1}{n} \sum_{i = 1}^{n} F_{j_{i}} (q^{- 1})

(9)

Multiplying Eq. to both sides of Eq. , we obtain the auxiliary system model:

A_{J_{*}} (q^{- J_{1}}, \dots, q^{- J_{n}}) y (k) = q^{- d} B_{J_{*}} (q^{- 1}) u (k)

(10)

where

\begin{aligned} A_{J_{*}} (q^{- J_{1}}, \dots, q^{- J_{n}}) & = 1 + \frac{a_{J_{1}, 1}}{n} q^{- J_{1}} + \dots + \frac{a_{J_{1}, n_{a}}}{n} q^{- n_{a} J_{1}} \\ ⋮ \\ + \frac{a_{J_{n}, 1}}{n} q^{- J_{n}} + \dots + \frac{a_{J_{n}, n_{a}}}{n} q^{- n_{a} J_{n}} \end{aligned}

Therefore, the auxiliary model can be directly identified since the required measurements now match the available measurements after rewriting Eq. in a predictor form. We will focus on the case where two output sensors of different rates are used in the remaining content for notational simplicity.

4 Recursive System Identification Under Collaborative Sensing

4.1 Recursive Least-Squares Formulation.

We present the recursive least-squares formulation for the case where two coprime sensors are deployed collaboratively in the scheme above for the output measurement. First, we design transformation polynomials for the characteristic equation of the original system model as follows:

\begin{aligned} F_{J_{1}} (q^{- 1}) & = \prod_{i = 1}^{n_{a}} (1 + (λ_{i} q)^{- 1} + \dots + (λ_{i} q)^{1 - J_{1}}) \\ F_{J_{2}} (q^{- 1}) & = \prod_{i = 1}^{n_{a}} (1 + (λ_{i} q)^{- 1} + \dots + (λ_{i} q)^{1 - J_{2}}) \end{aligned}

where J₁ and J₂ are coprime integers. Without loss of generality, we assume that a smaller index denotes the sensor with a faster sampling rate (i.e., J₁ < J₂). Summing up the two polynomials above and implementing the normalization in Eq. yield

F_{J_{*}} (q^{- 1}) = \frac{1}{2} [F_{J_{1}} (q^{- 1}) + F_{J_{2}} (q^{- 1})]

(11)

Next, multiplying the polynomial shown in Eq. to both sides of the original system model in Eq. yields

F_{J_{*}} (q^{- 1}) A (q^{- 1}) y (k) = q^{- d} F_{J_{*}} (q^{- 1}) B (q^{- 1}) u (k)

or

A_{J_{*}} (q^{- J_{1}}, q^{- J_{2}}) y (k) = q^{- d} B_{J_{*}} (q^{- 1}) u (k)

(12)

where

\begin{aligned} A_{J_{*}} (q^{- J_{1}}, q^{- J_{2}}) & = 1 + \frac{a_{J_{1}, 1}}{2} q^{- J_{1}} + \dots + \frac{a_{J_{1}, n_{a}}}{2} q^{- n_{a} J_{1}} \\ + \frac{a_{J_{2}, 1}}{2} q^{- J_{2}} + \dots + \frac{a_{J_{2}, n_{a}}}{2} q^{- n_{a} J_{2}} \\ B_{J_{*}} (q^{- 1}) & = b_{J_{*}, 1} q^{- 1} + b_{J_{*}, 2} q^{- 2} + \dots + b_{J_{*}, n_{a} (J_{2} - 1) + n_{b}} q^{- n_{a} (J_{2} - 1) - n_{b}} \end{aligned}

The order of

B_{J_{*}} (q^{- 1})

here comes from the sum of the order of B(q⁻¹), i.e., n_b, and that of

F_{J_{*}} (q^{- 1})

⁠, i.e., n_a(J₂ − 1). For simplicity, let χ = n_a(J₂ − 1) + n_b. Equation can be rearranged as

\begin{aligned} y (k) & = - \frac{a_{J_{1}, 1}}{2} y (k - J_{1}) - \dots - \frac{a_{J_{1}, n_{a}}}{2} y (k - n_{a} J_{1}) \\ - \frac{a_{J_{2}, 1}}{2} y (k - J_{2}) - \dots - \frac{a_{J_{2}, n_{a}}}{2} y (k - n_{a} J_{2}) \\ + b_{J_{*}, 1} u (k - d - 1) + \dots + b_{J_{*}, χ} u (k - d - χ) \end{aligned}

or in a vector form:

y (k) = θ^{T} ϕ (k)

(13)

where

(14)

From the vector form of the predictor function in Eq. , we see that the required historical output measurements are integer multiples of J₁ or J₂ steps behind the current instant k. At time instant iJ₁J₂, we have

\hat{y} (i J_{1} J_{2}) = ϕ^{T} (i J_{1} J_{2}) \hat{θ} (i J_{1} J_{2}), i = 0, 1, 2, \dots

Consider the performance index:

J_{k} = \sum_{i = 1}^{k} e (i ϰ)^{2} = \sum_{i = 1}^{k} {[y (i ϰ) - ϕ^{T} (i ϰ) \hat{θ} (k ϰ)]}^{2}

where

ϰ = J_{1} J_{2}

for brevity of notation. The solution

\hat{θ} (k ϰ)

can then be obtained using techniques from single-rate recursive least-squares, and the parameter adaptation algorithm (PAA) is as follows:

\hat{θ} ((k + 1) ϰ) = \hat{θ} (k ϰ) + F (k + 1) ϕ ((k + 1) ϰ) ϵ^{o} ((k + 1) ϰ)

(15)

F (k + 1) = F (k) - \frac{F (k) ϕ ((k + 1) ϰ) ϕ^{T} ((k + 1) ϰ) F (k)}{1 + ϕ^{T} ((k + 1) ϰ) F (k) ϕ ((k + 1) ϰ)}

(16)

where the a priori output estimation

{\hat{y}}^{o}

and a priori output estimation error

ϵ^{o}

are defined as

\begin{aligned} {\hat{y}}^{o} ((k + 1) ϰ) & = ϕ^{T} ((k + 1) ϰ) \hat{θ} (k ϰ) \\ ϵ^{o} ((k + 1) ϰ) & = y ((k + 1) ϰ) - {\hat{y}}^{o} ((k + 1) ϰ) \end{aligned}

The stability of the PAA follows from standard hyperstability analysis for system identification and adaptive control [6]. Between

k ϰ

and

(k + 1) ϰ

⁠, we keep the data asynchronous and hold the parameter estimate:

\hat{θ} (k ϰ + j) = \hat{θ} (k ϰ), j = 1, 2, \dots, ϰ - 1

⁠.

4.2 PAA Convergence and Identifiability Analysis.

Parameter convergence in standard system identification requires the model to be irreducible, meaning that the polynomial orders cannot be further reduced and there are no common factors between B(q⁻¹) and A(q⁻¹).

When the input signal persistently excites the system dynamics, the convergence condition reduces to the existence of the solution to the following Diophantine equation associated with the polynomial parameters [6]

\tilde{A} (q^{- 1}) B (q^{- 1}) - A (q^{- 1}) \tilde{B} (q^{- 1}) = 0

where

\tilde{A} (q^{- 1})

and

\tilde{B} (q^{- 1})

represent the difference between the ground truth and the estimated system polynomials. In deterministic cases, parameters converge to the true values when the only solution to the Diophantine equation is

\tilde{A} (q^{- 1}) = 0

and

\tilde{B} (q^{- 1}) = 0

⁠. For the proposed collaborative sensing, it can be shown that parameter convergence still holds due to the following lemma.

Diophantine multiplicative equations

Lemma 1

Given a polynomial of the following form:

F (z^{- 1}) = f_{0} + f_{1} z^{- 1} + f_{2} z^{- 2} + \dots + f_{m} z^{- m}

where not all f_i’s are zero, and

\begin{aligned} α (z^{- 1}) & = 1 + α_{1} z^{- 1} + \dots + α_{n} z^{- n} \\ β (z^{- 1}) & = β_{1} z^{- 1} + β_{2} z^{- 2} + \dots + β_{n} z^{- n} \end{aligned}

Then, the Diophantine equation

F (z^{- 1}) [α (z^{- 1}) \underset{unknown}{\underset{⏟}{σ (z^{- 1})}} + β (z^{- 1}) \underset{unknown}{\underset{⏟}{γ (z^{- 1})}}] = 0

(17)

has the unique zero solution for σ(z⁻¹) and γ(z⁻¹) (i.e., σ(z⁻¹) = 0 and γ(z⁻¹) = 0), if the numerators of α(z⁻¹) and β(z⁻¹) are coprime, and the orders of σ and γ are restricted to be less than n as follows:

\begin{aligned} σ (z^{- 1}) & = σ_{0} + σ_{1} z^{- 1} + \dots + σ_{n - 1} z^{- (n - 1)} \\ γ (z^{- 1}) & = γ_{0} + γ_{1} z^{- 1} + \dots + γ_{n - 1} z^{- (n - 1)} \end{aligned}

The proof is provided in the Appendix.

4.3 Parameter Recovery.

Recall the reparameterized system model:

y (k) = \frac{q^{- d} B (q^{- 1}) F_{J_{*}} (q^{- 1})}{A (q^{- 1}) F_{J_{*}} (q^{- 1})} u (k) = \frac{q^{- d} B_{J_{*}} (q^{- 1})}{A_{J_{*}} (q^{- J_{1}}, q^{- J_{2}})} u (k)

By applying the aforementioned RLS-system identification algorithm, the intermediate parameter vector, i.e., the coefficients of

B_{J_{*}} (q^{- 1})

and

A_{J_{*}} (q^{- J_{1}}, q^{- J_{2}})

⁠, can be identified directly. By removing the highest-order common factor from

B_{J_{*}} (q^{- 1})

and

A_{J_{*}} (q^{- J_{1}}, q^{- J_{2}})

⁠, the original fast model polynomials B(q⁻¹) and A(q⁻¹) can then be obtained.

5 Case Study

We present three cases with different system setups, including a practical example in motion controls. We assume that two sensors are deployed for the output data collection. For the first two simulation cases, J₁ and J₂ are 2 and 3 times slower than the input sampling rate, respectively, and an input pseudo-random binary sequence (PRBS) signal is generated at 1024 Hz. A sufficiently long time horizon is selected to ensure parameter convergence (within ten iterations). For the third motion control example implemented on a hard drive drive (HDD) benchmark, we assume that J₁ and J₂ are 9 and 13 times slower. The PRBS signal is generated at 50,400 Hz. Algorithm 1 outlines the implementation steps for the proposed algorithm.

Collaborative sensing RLS system identification

Algorithm 1

Input: $u (k)$ ⁠, $y (k)$ ⁠, $F$ ⁠, $J_{1}$ ⁠, $J_{2}$ ⁠, $n_{a}$ ⁠, $n_{b}$ ⁠, $d$

$θ$ ⁠, $ϕ \leftarrow$ $n_{a}$ ⁠, $n_{b}$ ⁠, $d$ while $t \leq t_{operation}$ do

ift mod $J_{1} \cdot J_{2} = 0$ then

Update $θ$ ⁠, $F$ ⁠; // refer to Eqs. (15) and (16)

Update $ϕ$ ⁠; // refer to Eq. (14)

else

ift mod $J_{1} = 0$ then

Update $ϕ_{J_{1}}$ ⁠, $ϕ_{u}$ // refer to Eq. (14)

end

ift mod $J_{2} = 0$ then

Update $ϕ_{J_{2}}$ ⁠, $ϕ_{u}$ // refer to Eq. (14)

end

$B_{J_{*}} (q^{- 1})$ ⁠, $A_{J_{*}} (q^{- J_{1}}, q^{- J_{2}}) \leftarrow θ$

$B (q^{- 1})$ ⁠, $A (q^{- 1}) \leftarrow B_{J_{*}} (q^{- 1})$ ⁠, $A_{J_{*}} (q^{- J_{1}}, q^{- J_{2}})$

Return: $B (q^{- 1})$ ⁠, $A (q^{- 1})$

5.1 Third-Order System.

Consider

y (k) = \frac{q^{- 2} + 0.5 q^{- 3}}{1 + 0.9 q^{- 1} + 0.26 q^{- 2} + 0.024 q^{- 3}} u (k)

where the poles are at −0.2, −0.3, and −0.4, and the zero is at −0.5. Rewrite the transfer function into the general form for system identification as follows:

G (q^{- 1}) ≜ \frac{q^{- d} B (q^{- 1})}{A (q^{- 1})} = \frac{q^{- 1} (q^{- 1} + 0.5 q^{- 2})}{1 + 0.9 q^{- 1} + 0.26 q^{- 2} + 0.024 q^{- 3}}

From the general form, we record the hyperparameters for the algorithm, which are d = 1, n_a = 3, and n_b = 2. The model parameters needed to be identified are

\tilde{B} (q^{- 1}) : [1.0, 0.5]

⁠,

\tilde{A} (q^{- 1}) : [1.0, 0.9, 0.26, 0.024]

⁠. The identified system response is shown in Fig. 2 and the parameters convergence of the auxiliary model is shown in Fig. 3. We also plotted the Nyquist frequencies of the individual sensors and observed the accurate model identification beyond the limitations of the individual sensors.

Fig. 2

View large Download slide

The frequency response comparison of the third-order system, and the identification beyond the Nyquist frequency

Fig. 3

View large Download slide

Illustration of parameter convergence of the auxiliary model for 30 iterations

5.2 Higher-Order System.

Consider a fourth-order system

G (q^{- 1}) = \frac{q^{- 1} (q^{- 1} + 1.5 q^{- 2} + 0.56 q^{- 3})}{1 + 1.4 q^{- 1} + 0.71 q^{- 2} + 0.154 q^{- 3} + 0.012 q^{- 4}}

where the poles are at −0.2, −0.3, −0.4, and −0.5, and the zeros are at −0.6 and −0.7. The hyperparameters are d = 1, n_a = 4, and n_b = 3. The identified parameters are

\tilde{B} (q^{- 1}) : [0.9999, 1.4999, 0.5602]

⁠,

\tilde{A} (q^{- 1}) : [1.0000, 1.3999, 0.7102, 0.1541, 0.0120]

⁠.

Figure 4 compares the original and identified system responses.

Fig. 4

View large Download slide

Higher-order system frequency response comparison, and the identification beyond the Nyquist frequency

5.3 Hard Drive Drive Benchmark System.

Consider the major first two modes in the voice coil motor of an HDD benchmark system [20]:

G (q^{- 1}) = \frac{0.00033 q^{- 1} + 0.003 q^{- 2} + 0.00264 q^{- 3} + 0.00019 q^{- 4}}{1 - 3.559 q^{- 1} + 5.091 q^{- 2} - 3.506 q^{- 3} + 0.9739 q^{- 4}}

where the poles are at 0.7792 + 0.6055i, 0.7792 − 0.6055i, 1, and 0.9999, and the zeros are at −8.2069, −0.0831, and −0.8835. The plant has common characteristics that relate torque/force to position in motion control. The hyperparameters are d = 0, n_a = 4, n_b = 4, M = 9, and N = 13. The identified parameters are

\tilde{B} (q^{- 1}) : [0.00033, 0.00303, 0.00265, 0.00020], \tilde{A} (q^{- 1}) :

[1, - 3.55851, 5.09094, - 3.50634, 0.97391]

⁠. Figure 5 compares the original and identified system responses.

Fig. 5

View large Download slide

HDD benchmark system frequency response comparison, and the identification beyond the Nyquist criterion

In all cases, the proposed algorithm was observed to have accurately identified the underlying system dynamics beyond the individual sensor’s Nyquist sampling limit.

6 Conclusion and Future Work

This paper presented a novel framework for non-uniformly sampled system identification based on the proposed coprime collaborative sensing and the RLS-based algorithm. Leveraging a polynomial transformation and characteristics of coprime numbers, we showed how the algorithm can recover fast system models beyond the Nyquist frequencies of multiple slow sensors. Example applications in motion control illustrate the effectiveness of the process. Future work includes optimal sensor rate selection, minimum data requirements, and addressing noise in stochastic environments.

Footnote

1

Paper presented at the 2023 Modeling, Estimation, and Control Conference (MECC 2023), Lake Tahoe, NV, Oct. 2–5, Paper No. MECC2023-88.

Acknowledgment

This material is based upon work supported by the National Science Foundation under Grants Nos. CMMI-1953155 and CMMI-2141293. The opinions, findings, and conclusions, or recommendations expressed are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Appendix: Diophantine Multiplicative Equations Proof

Proof

Given

F (z^{- 1}) = f_{0} + f_{1} z^{- 1} + f_{2} z^{- 2} + \dots + f_{m} z^{- m}

let

\tilde{η} (z^{- 1}) = α (z^{- 1}) σ (z^{- 1}) + β (z^{- 1}) γ (z^{- 1})

(A1)

where

\begin{aligned} α (z^{- 1}) & = 1 + α_{1} z^{- 1} + \dots + α_{n} z^{- n} \\ β (z^{- 1}) & = β_{1} z^{- 1} + β_{2} z^{- 2} + \dots + β_{n} z^{- n} \end{aligned}

and

\begin{aligned} σ (z^{- 1}) & = σ_{0} + σ_{1} z^{- 1} + \dots + σ_{n - 1} z^{- (n - 1)} \\ γ (z^{- 1}) & = γ_{0} + γ_{1} z^{- 1} + \dots + γ_{n - 1} z^{- (n - 1)} \end{aligned}

are unknown a priori.

We show that

F (z^{- 1}) \tilde{η} (z^{- 1}) = 0

holds only when

\tilde{η} (z^{- 1}) = 0

and subsequently σ(z⁻¹) and γ(z⁻¹) must all be zero. After forming the Sylvester matrix, Eq. is equivalent to the linear equality:

S [\begin{matrix} 1 \\ α_{1} \\ ⋮ \\ α_{n - 1} \\ γ_{0} \\ γ_{1} \\ ⋮ \\ γ_{n - 1} \end{matrix}] = [\begin{matrix} {\tilde{η}}_{0} \\ {\tilde{η}}_{1} \\ ⋮ \\ {\tilde{η}}_{n - 1} \\ {\tilde{η}}_{n} \\ {\tilde{η}}_{n + 1} \\ ⋮ \\ {\tilde{η}}_{2 n - 1} \end{matrix}]

where

S = {[\begin{matrix} 1 & 0 & \dots & 0 & 0 & \dots & \dots & 0 \\ α_{1} & ⋱ & ⋱ & ⋮ & β_{1} & ⋱ & ⋮ \\ ⋮ & ⋱ & ⋱ & 0 & ⋮ & ⋱ & ⋱ & ⋮ \\ ⋮ & ⋱ & 1 & ⋮ & ⋱ & 0 \\ α_{n - 1} & α_{1} & β_{n - 1} & β_{1} \\ α_{n} & ⋱ & ⋮ & β_{n} & ⋱ & ⋮ \\ 0 & ⋱ & ⋱ & ⋮ & 0 & ⋱ & ⋱ & ⋮ \\ ⋮ & ⋱ & ⋱ & α_{n - 1} & ⋮ & ⋱ & ⋱ & β_{n - 1} \\ 0 & \dots & 0 & α_{n} & 0 & \dots & 0 & β_{n} \end{matrix}]}_{2 n, 2 n}

Then, the coefficients of the filter product

F (z^{- 1}) \tilde{η} (z^{- 1}) = η_{0} + η_{1} z^{- 1} + \dots + η_{m + 2 n - 1} z^{- m - 2 n + 1}

satisfies

\underset{η_{*}}{\underset{⏟}{[\begin{matrix} η_{0} \\ η_{1} \\ ⋮ \\ η_{m + 2 n - 1} \end{matrix}]}} = F_{*} \underset{{\tilde{η}}_{*}}{\underset{⏟}{[\begin{matrix} {\tilde{η}}_{0} \\ {\tilde{η}}_{1} \\ ⋮ \\ {\tilde{η}}_{2 n - 1} \end{matrix}]}} = F_{*} S \underset{ξ_{*}}{\underset{⏟}{[\begin{matrix} σ_{0} \\ σ_{1} \\ ⋮ \\ σ_{n - 1} \\ γ_{0} \\ ⋮ \\ γ_{n - 1} \end{matrix}]}}

where

F_{*} = {[\begin{matrix} f_{0} & 0 & \dots & 0 \\ f_{1} & ⋱ & ⋱ & ⋮ \\ ⋮ & ⋱ & ⋱ & 0 \\ ⋮ & ⋱ & f_{0} \\ f_{m - 1} & f_{1} \\ f_{m} & ⋱ & ⋮ \\ 0 & ⋱ & ⋱ & ⋮ \\ ⋮ & ⋱ & ⋱ & f_{m - 1} \\ 0 & \dots & 0 & f_{m} \end{matrix}]}_{m + 2 n - 1, 2 n}

and all columns of

F_{*}

are linearly independent. Thus, if

η_{*} = 0

⁠,

{\tilde{η}}_{*}

must be a zero vector. If the numerators of α(z⁻¹) and β(z⁻¹) are copirme, S will be nonsingular and thus

{\tilde{η}}_{*}

and

ξ_{*}

form a one-to-one mapping. The unique solution to Eq. is thus

ξ_{*} = 0

⁠.

References

1.

Isermann

,

R.

,

1989

,

Digital Control Systems

,

Springer Science & Business Media

,

Heidelberg

.

Google Scholar

Crossref

2.

Gambier

,

A.

,

2004

, “

Real-Time Control Systems: A Tutorial

,”

5th Asian Control Conference (IEEE Cat. No. 04EX904)

,

Melbourne, Australia

,

July 20–23

, Vol. 2, IEEE, pp.

1024

–

1031

.

3.

Khan

,

S.

,

Goodall

,

R. M.

, and

Dixon

,

R.

,

2013

, “

Non-Uniform Sampling Strategies for Digital Control

,”

Int. J. Syst. Sci.

,

44

(

12

), pp.

2234

–

2254

.

Google Scholar

Crossref

4.

Ding

,

F.

,

Qiu

,

L.

, and

Chen

,

T.

,

2009

, “

Reconstruction of Continuous-Time Systems From Their Non-uniformly Sampled Discrete-Time Systems

,”

Automatica

,

45

(

2

), pp.

324

–

332

.

Google Scholar

Crossref

5.

Albertos

,

P.

, and

Crespo

,

A.

,

1999

, “

Real-Time Control of Non-uniformly Sampled Systems

,”

Control Eng. Pract.

,

7

(

4

), pp.

445

–

458

.

Google Scholar

Crossref

6.

Landau

,

I. D.

,

Lozano

,

R.

,

M’Saad

,

M.

, and

Karimi

,

A.

,

2011

,

Adaptive Control: Algorithms, Analysis and Applications

,

Springer Science & Business Media

,

London

.

Google Scholar

Crossref

7.

Powell

,

M. J. D.

,

1981

,

Approximation Theory and Methods

,

Cambridge University Press

,

Cambridge

.

Google Scholar

Crossref

8.

Raghavan

,

H.

,

Gopaluni

,

R. B.

,

Shah

,

S.

,

Pakpahan

,

J.

,

Patwardhan

,

R.

, and

Robson

,

C.

,

2005

, “

Gray-Box Identification of Dynamic Models for the Bleaching Operation in a Pulp Mill

,”

J. Process Control

,

15

(

4

), pp.

451

–

468

.

Google Scholar

Crossref

9.

Amirthalingam

,

R.

,

Sung

,

S. W.

, and

Lee

,

J. H.

,

2000

, “

Two-Step Procedure for Data-Based Modeling for Inferential Control Applications

,”

AIChE J.

,

46

(

10

), pp.

1974

–

1988

.

Google Scholar

Crossref

10.

Sircar

,

P.

, and

Sarkar

,

T. K.

,

1988

, “

System Identification From Nonuniformly Spaced Signal Measurements

,”

Signal Process.

,

14

(

3

), pp.

253

–

268

.

Google Scholar

Crossref

11.

Gillberg

,

J.

, and

Ljung

,

L.

,

2010

, “

Frequency Domain Identification of Continuous-Time Output Error Models, Part II: Non-uniformly Sampled Data and B-Spline Output Approximation

,”

Automatica

,

46

(

1

), pp.

11

–

18

.

Google Scholar

Crossref

12.

Dempster

,

A. P.

,

Laird

,

N. M.

, and

Rubin

,

D. B.

,

1977

, “

Maximum Likelihood From Incomplete Data Via the EM Algorithm

,”

J. R. Stat. Soc.: Ser. B (Methodol.)

,

39

(

1

), pp.

1

–

22

.

13.

Shumway

,

R. H.

, and

Stoffer

,

D. S.

,

1982

, “

An Approach to Time Series Smoothing and Forecasting Using the EM Algorithm

,”

J. Time Series Anal.

,

3

(

4

), pp.

253

–

264

.

Google Scholar

Crossref

14.

Isaksson

,

A. J.

,

1993

, “

Identification of ARX-Models Subject to Missing Data

,”

IEEE Trans. Automat. Contr.

,

38

(

5

), pp.

813

–

819

.

Google Scholar

Crossref

15.

Jones

,

R. H.

,

1980

, “

Maximum Likelihood Fitting of ARMA Models to Time Series With Missing Observations

,”

Technometrics

,

22

(

3

), pp.

389

–

395

.

Google Scholar

Crossref

16.

Li

,

D.

,

Shah

,

S. L.

, and

Chen

,

T.

,

2001

, “

Identification of Fast-Rate Models From Multirate Data

,”

Int. J. Control

,

74

(

7

), pp.

680

–

689

.

Google Scholar

Crossref

17.

Li

,

D.

,

Shah

,

S. L.

, and

Chen

,

T.

,

2002

, “

Analysis of Dual-Rate Inferential Control Systems

,”

Automatica

,

38

(

6

), pp.

1053

–

1059

.

Google Scholar

Crossref

18.

Zhu

,

Y.

,

Telkamp

,

H.

,

Wang

,

J.

, and

Fu

,

Q.

,

2009

, “

System Identification Using Slow and Irregular Output Samples

,”

J. Process. Control

,

19

(

1

), pp.

58

–

67

.

Google Scholar

Crossref

19.

Lu

,

W.

, and

Grant Fisher

,

D.

,

1988

, “

Output Estimation With Multi-Rate Sampling

,”

Int. J. Control

,

48

(

1

), pp.

149

–

160

.

Google Scholar

Crossref

20.

Atsumi

,

T.

,

2023

, “

Magnetic-Head Positioning Control System in HDDs

,” https://www.mathworks.com/matlabcentral/fileexchange/111515-magnetic-head-positioning-control-system-in-hdds, MATLAB Central File Exchange, Retrieved July 20, 2023.

A Recursive System Identification With Non-Uniform Temporal Feedback Under Coprime Collaborative Sensing^¹

Abstract

1 Introduction

2 Preliminaries

3 Proposed Coprime Collaborative Sensing and Model Reparameterization

4 Recursive System Identification Under Collaborative Sensing

4.1 Recursive Least-Squares Formulation.

4.2 PAA Convergence and Identifiability Analysis.

Diophantine multiplicative equations

4.3 Parameter Recovery.

5 Case Study

Collaborative sensing RLS system identification

5.1 Third-Order System.

5.2 Higher-Order System.

5.3 Hard Drive Drive Benchmark System.

6 Conclusion and Future Work

Footnote

Acknowledgment

Conflict of Interest

Data Availability Statement

Appendix: Diophantine Multiplicative Equations Proof

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

A Recursive System Identification With Non-Uniform Temporal Feedback Under Coprime Collaborative Sensing1

Abstract

1 Introduction

2 Preliminaries

3 Proposed Coprime Collaborative Sensing and Model Reparameterization

4 Recursive System Identification Under Collaborative Sensing

4.1 Recursive Least-Squares Formulation.

4.2 PAA Convergence and Identifiability Analysis.

Diophantine multiplicative equations

4.3 Parameter Recovery.

5 Case Study

Collaborative sensing RLS system identification

5.1 Third-Order System.

5.2 Higher-Order System.

5.3 Hard Drive Drive Benchmark System.

6 Conclusion and Future Work

Footnote

Acknowledgment

Conflict of Interest

Data Availability Statement

Appendix: Diophantine Multiplicative Equations Proof

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only

A Recursive System Identification With Non-Uniform Temporal Feedback Under Coprime Collaborative Sensing^¹