Parameter identification of nonlinear multirate time-delay system with uncertain output delays

Abstract

The joint parameter and time-delay estimation problems for a class of nonlinear multirate time-delay system with uncertain output delays are addressed in this paper. The practical process typically has time-delay properties and the process data are often multirate, sampled with output data inevitably corrupted by uncertain delays. The linear parameter varying (LPV) finite impulse response (FIR) multirate time-delay model is initially built to describe the considered system. The problems of over-parameterization and the existence of both continuous model parameters and discrete time-delays have made the conventional maximum likelihood difficult to solve the considered problems. In order to handle these problems, the joint parameter and time-delay estimation for the LPV FIR multirate time-delay model are formulated in the expectation-maximization scheme, and the algorithm to estimate the model parameters and time-delays is derived, simultaneously based on multirate process data. The efficacy of the proposed method is verified through a numerical simulation and a practical chemical plant.

Keywords

System identification nonlinear time-delay model multirate data the continuous stirred tank reactor system expectation-maximization algorithm

Introduction

Modern industrial processes typically perform complex production tasks and exhibit complicated process dynamics, leading to a nonlinear system (Sun et al., 2013; Yao et al., 2014a; Yin et al., 2016; Zhao et al., 2015). The uncertainty of the nonlinear system model can greatly limit the control performance, and therefore the system model should be derived first (Yao et al., 2014b). The nonlinear system identification has been a longstanding and extremely challenging theoretical and engineering problem, and many approaches have been developed to solve it. Among the reported nonlinear system identification approaches, linear parameter varying (LPV) modelling has become one of the most attractive methods due to its relatively simple model structure and close connection with LPV control analysis and synthesis (Toth et al., 2012; Zhong et al., 2010, 2015).

The LPV model has a parametrized linear model structure and the model parameters are typically described as a linear combination of certain basis functions with respect to time-varying scheduling variables. The LPV identification methods can be classified into state-space-based approaches and transfer function-based approaches, according to the adopted model structure. The identification of the LPV state-space model was considered in Wills and Ninness (2011) and an expectation-maximization (EM) algorithm-based method was developed. The assumption of a known state sequence was omitted and the unknown states were estimated by using Kalman smooth in the expectation step. The subspace identification of the LPV model was handled in van Wingerden and Verhaegen (2009). A predictor containing the LPV equivalent of the Markov parameters was constructed based on an introduced factorization, and the state sequence was estimated by using this predictor. The system matrices were then recovered from the estimated state sequence. However, the state-space-based methods are normally computationally complex due to the estimation of an unknown state sequence and the estimated system matrices not being unique. The prediction error method was extended to LPV identification in Zhao et al. (2012). The cost function was built based on one-step-ahead prediction error and the numerical optimization algorithm was used to optimize the cost function, to find the optimal parameter estimates. In Laurain et al. (2009), the instrumental variable method was proposed for LPV Box–Jenkins model parameter estimation. The robust identification of the LPV system was presented in Yang et al. (2015). The robust LPV model was constructed with the $t$ -distribution, and the formulas to estimate the model parameters iteratively were derived.

The practical process data set is often not ideal, and the multirate sampled data and time-delays are common (Gopaluni, 2010; Sun et al., 2011; Yan et al., 2016). The non-ideal data have made the process modelling and control very challenging (Sun et al., 2015; Yin et al., 2017). The quality variables and process variables are usually sampled in slow-rate and fast-rate, respectively, due to sensor performance limitations, practical system requirements, etc., resulting in multirate data. In Xie et al. (2013), the multirate process identification with random delays was considered. An initial finite impulse response (FIR) model was identified and the output error (OE) model was then identified based on complete data, with missing output data compensated with output estimates of the FIR model. The multirate system identification can be seen as a special case of the missing data treatment problem. In Yang and Yin (2017), the variational Bayesian approach to identify the linear systems with missing data was proposed, and the missing data were compensated with its distribution. This process often has time-delay properties and the output data are inevitably corrupted by uncertain delays due to transmission, communication limitations (Zhu et al., 2016), etc. The mode-dependent time-varying delays were considered in Zhang et al. (2015) when dealing with the asynchronous state estimation issue for a class of Markov jump neural networks in the $l_{2}$ - $l_{\infty}$ sense. A non-iterative method was proposed in Masmoudi et al. (2011) to estimate the time-delay and other model parameters, subsequently in the maximum likelihood estimation framework. The grid searching method was adopted in Weyer (2001) to estimate the time-delay of the open water channel.

This paper aims to develop a joint parameter and time-delay estimation algorithm for a nonlinear multirate time-delay system. The LPV modelling approach is adopted and the LPV FIR multirate time-delay model is considered. The multirate data are common and the output data are inevitably corrupted by uncertain delays. To our best knowledge, the LPV time-delay identification with multirate data and uncertain output delays has not been reported in the literature. The conventional maximum likelihood method cannot be used to handle the estimation problems. In order to solve the considered problems, the probability model is constructed with Gaussian distribution. Mathematical formulation is performed in the EM algorithm scheme, in which the time-delay process property and uncertain output delays are dealt with simultaneously.

The paper is organized as follows. The problem statement is given in the following section. The EM algorithm is briefly introduced in section ‘Mathematical formation of the nonlinear multirate time-delay system identification problems in the EM framework’. The performance verification of the proposed method is presented in section ‘Verification’ and the final section concludes the work.

Problem statement

Consider the nonlinear multirate system described by the following linear parameter varying (LPV) finite impulse response (FIR) multirate time-delay model

\begin{matrix} x_{k} = G (z_{k}, q^{- 1}) u_{k - τ} \\ y_{t_{i}} = x_{t_{i} - λ_{i}} + ε_{t_{i}} \end{matrix}

(1)

where $z_{k}$ and $u_{k}$ are the fast-rate scheduling variable and input variable, respectively, with sampling period $Δ t$ , $τ$ is the unknown system time-delay, $x_{k}$ is the unmeasurable noise-free output variable, $q^{- 1}$ is the back shift operator, $ε_{t_{i}}$ is the Gaussian white noise, $y_{t_{i}}$ is the slow-rate output variable that has a sampling period that is the integral multiple of $Δ t$ , $λ_{i}$ is the uncertain output delay, and $t_{i}$ is the time instant that the output data is available. $G (z_{k}, q^{- 1})$ is the fast-rate transfer function that has a FIR model structure with time-varying parameter

G (z_{k}, q^{- 1}) = c_{1} (z_{k}) q^{- 1} + c_{2} (z_{k}) q^{- 2} + \dots + c_{n} (z_{k}) q^{- n}

(2)

where the parameters are a linear combination of meromorphic functions of $z_{k}$ at time instant $k Δ t$

c_{j} (z_{k}) = c_{j, 0} + \sum_{m = 1}^{M} c_{j, m} f_{m} (z_{k}), j = 1, \dots, n

(3)

${f_{m} (z_{k})}_{m = 1, \dots, M}$ are meromorphic functions, and $n$ and $M$ are the orders of the LPV FIR model and the parameter functions, respectively.

Based on equation (1), it can be further derived that

x_{t_{i} - λ_{i}} = q^{- λ_{i}} x_{t_{i}} = G (z_{t_{i} - λ_{i}}, q^{- 1}) u_{t_{i} - λ_{i} - τ}

(4)

Then, the model (1) is rewritten into a compact form

y_{t_{i}} = G (z_{t_{i} - λ_{i}}, q^{- 1}) u_{t_{i} - λ_{i} - τ} + ε_{t_{i}}

(5)

The measurement noise $ε_{t_{i}}$ is assumed to follow the Gaussian distribution

ε_{t_{i}} ~ N (ε_{t_{i}} | 0, δ)

(6)

where $δ$ is the variance. According to equation (5), the output is also Gaussian distributed

y_{t_{i}} ~ N (y_{t_{i}} | G (z_{t_{i} - λ_{i}}, q^{- 1}) u_{t_{i} - λ_{i} - τ}, δ)

(7)

Denote the recorded slow-rate output data ${y_{t_{i}}}_{i = 1, \dots, N}$ , fast-rate input data ${u_{k}}_{k = 1, \dots, L}$ , and scheduling variable data ${z_{k}}_{k = 1, \dots, L}$ as $y_{t_{1} : t_{N}}$ , $u_{1 : L}$ , and $z_{1; L}$ , respectively. Denote the output delays ${λ_{i}}_{i = 1, \dots, N}$ as $λ_{1 : N}$ . The considered parameter estimation problem is to estimate the model parameters $θ$ , system time-delay $τ$ , and noise variance $δ$ based on the observed multirate process data $Y = {y_{t_{1}, t_{N}}, u_{1 : L}, z_{1 : L}}$ with output data corrupted by the uncertain output delays $λ_{1 : N}$ .

Mathematical formulation of the nonlinear multirate time-delay system identification problems in the EM framework

The conventional maximum likelihood method is commonly used to solve the formulated estimation problem through maximizing the log likelihood function of the observed data set. The log likelihood function of $Y$ is

\begin{matrix} L (y_{t_{1} : t_{N}}, u_{1 : L}, z_{1 : L} | θ, δ) \\ = - \frac{N}{2} \log δ - \sum_{i = 1}^{N} \frac{{(y_{t_{i}} - G (z_{t_{i} - λ_{i}}, q^{- 1}) u_{t_{i} - λ_{i} - τ})}^{2}}{2 δ} \\ - \frac{N}{2} \log 2 π \end{matrix}

(8)

Note that $L (y_{t_{1} : t_{N}}, u_{1 : L}, z_{1 : L} | θ, δ)$ contains both the continuous variables $θ$ and $δ$ and the discrete variables $τ$ and $λ_{1 : N}$ , and has as many unknown parameters as the observed data, so it cannot be optimized directly to find the parameter estimates. Here, the maximum likelihood estimation problem is formulated and solved in the EM algorithm framework.

Brief introduction to the EM algorithm

The EM algorithm simplifies the maximum likelihood estimation problem by introducing the latent variables and missing variables, and conducts the parameter optimization by alternating the expectation step (E-step) and the maximization step (M-step). The EM algorithm begins with the construction of the observed data set $Y$ and the missing data set $Z$ , leading to the complete data set ${Z, Y}$ . The moments of functions of latent variables and missing variables are calculated in the E-step, and the marginal likelihood of the complete data is maximized with respect to unknown parameters in the M-step to find new parameter estimates. The overflows of the EM algorithm are expressed as follows (McLachlan and Krishnan, 2007).

1. Initialize parameters and set $s = 1$ .

2. E-step: Calculate the following $Q$ -function based on current parameter estimates $Ω^{s}$

Q (Ω | Ω^{s}) = E_{Z | Y, Ω^{s}} {\log p (Z, Y | Ω)}

(9)

3. M-step: Maximize the $Q$ -function to find new parameter estimates

Ω^{s + 1} = \arg max_{Ω} Q (Ω | Ω^{s})

(10)

4. Once convergence of the algorithm is achieved, stop (else, set $s = s + 1$ and return to step 2).

Mathematical derivations

The observed data set collects all the available data and is then built as $Y = {y_{t_{1} : t_{N}}, u_{1 : L}, z_{1 : L}}$ . The system delay $τ$ is treated as a discrete random variable and is assumed to be uniformly distributed in the integer range $[d_{1}, d_{2}]$ . The uncertain output delay $λ_{i}$ is also treated as a discrete random variable and assumed to follow a uniform distribution in the integer range $[γ_{1}, γ_{2}]$ . In order to make the formulated parameter and time-delay estimation problems tractable, the system delay $τ$ and uncertain output delays $λ_{1 : N}$ are treated as latent variables and are incorporated into the missing data set $Z = {τ, λ_{1 : N}}$ . The unknown model parameters are $Ω = {θ, δ}$ .

It is assumed that the delay ranges $[d_{1}, d_{2}]$ and $[γ_{1}, γ_{2}]$ can be roughly pre-specified. In practice, the delay range can be easily set according to physical laws. For example, if the transmission distance and speed range are known, the delay range can be calculated. The delay range can also be determined by statistical tests or set to an arbitrary range that is large enough to contain the true time-delay.

The log likelihood function of ${Z, Y}$ is derived based on the probability chain rule

\begin{matrix} \log p (Z, Y | Ω) = \log p (y_{t_{1} : t_{N}}, u_{1 : L}, z_{1 : L}, τ, λ_{1 : N} | Ω) \\ = \log p (y_{t_{1} : t_{N}} | u_{1 : L}, z_{1 : L}, τ, λ_{1 : N}, Ω) \\ + \log p (τ | u_{1 : L}, z_{1 : L}, λ_{1 : N}, Ω) \\ + \log p (λ_{1 : N} | u_{1 : L}, z_{1 : L}, Ω) \\ + \log p (u_{1 : L}, z_{1 : L} | Ω) \end{matrix}

(11)

The term $\log p (y_{t_{1} : t_{N}} | u_{1 : L}, z_{1 : L}, τ, λ_{1 : N}, Ω)$ in (11) is further derived as

\begin{matrix} \log p (y_{t_{1} : t_{N}} | u_{1 : L}, z_{1 : L}, τ, λ_{1 : N}, Ω) \\ = \log p (y_{t_{N}} | y_{1 : t_{N} - 1}, u_{1 : L}, z_{1 : L}, τ, λ_{1 : N}, Ω) \\ + \log p (y_{1 : t_{N} - 1} | u_{1 : L}, z_{1 : L}, τ, λ_{1 : N}, Ω) \\ = \sum_{i = 1}^{N} \log p (y_{t_{i}} | y_{1 : t_{i} - 1}, u_{1 : L}, z_{1 : L}, τ, λ_{1 : N}, Ω) \\ = \sum_{i = 1}^{N} \log p (y_{t_{i}} | u_{1 : t_{i} - 1}, z_{t_{i}}, τ, λ_{i}, Ω) \end{matrix}

(12)

The derivation is based on the fact that output $y_{t_{i}}$ depends only on the input data $u_{1 : t_{i} - 1}$ , current scheduling variable $z_{t_{i}}$ , the system time-delay $τ$ , output delay $λ_{i}$ , and the model parameters $Ω$ .

The term $\log p (τ | u_{1 : L}, z_{1 : L}, λ_{1 : N}, Ω)$ in (11) is simplified to

\log p (τ | u_{1 : L}, z_{1 : L}, λ_{1 : N}, Ω) = \log p (τ),

(13)

because the system time-delay is independent of all the data and other parameters.

The term $p (λ_{1 : N} | u_{1 : L}, z_{1 : L}, Ω)$ in (11) is expanded to

\begin{matrix} \log p (λ_{1 : N} | u_{1 : L}, z_{1 : L}, Ω) \\ = \sum_{i = 1}^{N} \log p (λ_{i} | λ_{1 : i - 1}, u_{1 : L}, z_{1 : L}, Ω) \\ = \sum_{i = 1}^{N} \log p (λ_{i}) \end{matrix}

(14)

The above formula is derived using the fact that the output delays are independent of each other, and are also independent of all the data and other model parameters.

The final term in (11) is a constant and is denoted as

\log p (u_{1 : L}, z_{1 : L} | Ω) = C_{1}

(15)

The log likelihood function of ${Z, Y}$ in (11) is written as

\begin{matrix} \log p (Z, Y | Ω) \\ = \sum_{i = 1}^{N} \log p (y_{t_{i}} | u_{1 : t_{i} - 1}, z_{t_{i}}, τ, λ_{i}, Ω) \\ + \log p (τ) + \sum_{i = 1}^{N} \log p (λ_{i}) + C_{1} \end{matrix}

(16)

The $Q$ -function in (9) is then expressed as

\begin{matrix} Q (Ω | Ω^{s}) = E_{τ, λ_{1 : N} | Y, Ω^{s}} {- \frac{N}{2} \log 2 π - \frac{N}{2} \log δ \\ - \sum_{i = 1}^{N} \frac{{(y_{t_{i}} - G (z_{t_{i} - λ_{i}}, q^{- 1}) u_{t_{i} - λ_{i} - τ})}^{2}}{2 δ} \\ + \log p (τ) + \sum_{i = 1}^{N} \log p (λ_{i}) + C_{1}} \end{matrix}

(17)

Taking the conditional expectations of $Q (Ω | Ω^{s})$ with respect to the discrete random variables $τ$ and $λ_{1 : N}$ , the $Q$ -function is written as

Q (Ω | Ω^{s}) = J (θ, δ) + C

(18)

where

\begin{matrix} J (θ, δ) = - \frac{N}{2} \log δ - \frac{1}{2 δ} \sum_{d = d_{1}}^{d_{2}} p (τ = d | Y, Ω^{s}) \sum_{i = 1}^{N} \sum_{γ = γ_{1}}^{γ_{2}} \\ p (λ_{i} = γ | τ = d, Y, Ω^{s}) (y_{t_{i}} - G (z_{t_{i} - γ}, q^{- 1}) u_{t_{i} - γ - d})^{2} \end{matrix}

(19)

and

\begin{matrix} C = \sum_{d = d_{1}}^{d_{2}} p (τ = d | Y, Ω^{s}) \log p (τ = d) + \sum_{d = d_{1}}^{d_{2}} \\ p (τ = d | Y, Ω^{s}) \sum_{i = 1}^{N} \sum_{γ = γ_{1}}^{γ_{2}} p (λ_{i} = γ | τ = d, Y, Ω^{s}) \\ \log p (λ_{i} = γ) - \frac{N}{2} \log 2 π + C_{1} \end{matrix}

(20)

In order to obtain the complete form of the $Q$ -function, the following posterior distributions should be calculated

1 . p (τ = d | Y, Ω^{s})

2 . p (λ_{i} = γ | τ = d, Y, Ω^{s})

The posterior probability of $τ = d$ is

\begin{array}{l} p (τ = d | Y, Ω^{s}) \\ = \frac{\prod_{i = 1}^{N} p (y_{t_{i}} | τ = d, z_{t_{i}}, u_{1 : t_{i} - 1}, Ω^{s}) p (τ = d)}{\sum_{d = d_{1}}^{d_{2}} \prod_{i = 1}^{N} p (y_{t_{i}} | τ = d, z_{t_{i}}, u_{1 : t_{i} - 1}, Ω^{s}) p (τ = d)} \\ = {\prod_{i = 1}^{N} \sum_{γ = γ_{1}}^{γ_{2}} p (y_{t_{i}} | τ = d, λ_{i} = γ, z_{t_{i}}, u_{1 : t_{i} - 1}, Ω^{s}) \times p (λ_{i} = γ) p (τ = d)} / {\sum_{d = d_{1}}^{d_{2}} \prod_{i = 1}^{N} \sum_{γ = γ_{1}}^{γ_{2}} p (y_{t_{i}} | τ = d, \\ λ_{i} = γ, z_{t_{i}}, u_{1 : t_{i} - 1}, Ω^{s}) p (λ_{i} = γ) p (τ = d)} \end{array}

(21)

The posterior probability of $λ_{i} = γ$ is

\begin{matrix} p (λ_{i} = γ | τ = d, Y, Ω^{s}) \\ = {p (y_{t_{i}} | τ = d, λ_{i} = γ, u_{1 : t_{i} - 1}, z_{t_{i}}, Ω^{s}) p (λ_{i} = γ)} \\ / {\sum_{γ = γ_{1}}^{γ_{2}} p (y_{t_{i}} | τ = d, λ_{i} = γ, u_{1 : t_{i} - 1}, z_{t_{i}}, Ω^{s}) \\ \times p (λ_{i} = γ)} \end{matrix}

(22)

Based on equations (1) and (5), the system model can be written in linear regression form

y_{t_{i}} = ϕ_{τ, t_{i} - λ_{i}}^{T} θ + ε_{t_{i}}

(23)

where

\begin{matrix} ϕ_{τ, t_{i} - λ_{i}} = [u_{t_{i} - λ_{i} - τ - 1} f_{1} (z_{t_{i} - λ_{i}}) u_{t_{i} - λ_{i} - τ - 1} \dots \\ f_{M} (z_{t_{i} - λ_{i}}) u_{t_{i} - λ_{i} - τ - 1} u_{t_{i} - λ_{i} - τ - 2} f_{1} (z_{t_{i} - λ_{i}}) \\ \times u_{t_{i} - λ_{i} - τ - 2} \dots f_{M} (z_{t_{i} - λ_{i}}) u_{t_{i} - λ_{i} - τ - 2} \dots \\ u_{t_{i} - λ_{i} - τ - n} \dots f_{M} (z_{t_{i} - λ_{i}}) u_{t_{i} - λ_{i} - τ - n}]^{T} \end{matrix}

(24)

and

θ = {[c_{1, 0} c_{1, 1} \dots c_{1, M} c_{2, 0} c_{2, 1} \dots c_{2, M} \dots c_{n, 0} \dots c_{n, M}]}^{T}

(25)

Therefore, the cost function $J (θ, δ)$ is transformed into

\begin{matrix} J (θ, δ) = - \frac{N}{2} \log δ - \frac{1}{2 δ} \sum_{d = d_{1}}^{d_{2}} p (τ = d | Y, Ω^{s}) \sum_{i = 1}^{N} \sum_{γ = γ_{1}}^{γ_{2}} \\ p (λ_{i} = γ | τ = d, Y, Ω^{s}) {(y_{t_{i}} - ϕ_{d, t_{i} - γ}^{T} θ)}^{2} \end{matrix}

(26)

To obtain the formula to update the variance estimate, the derivative is taken over $J (θ, δ)$ with respect to $δ$ and set to zero

\begin{matrix} \frac{\partial J (θ, δ)}{\partial δ} = - \frac{N}{2 δ} + \frac{1}{2 δ^{2}} \sum_{d = d_{1}}^{d_{2}} p (τ = d | Y, Ω^{s}) \sum_{i = 1}^{N} \sum_{γ = γ_{1}}^{γ_{2}} \\ p (λ_{i} = γ | τ = d, Y, Ω^{s}) {(y_{t_{i}} - ϕ_{d, t_{i} - γ}^{T} θ)}^{2} \\ = 0 \end{matrix}

(27)

Solving the above equation, it is derived that

\begin{matrix} δ^{s + 1} = \frac{1}{N} \sum_{d = d_{1}}^{d_{2}} p (τ = d | Y, Ω^{s}) \sum_{i = 1}^{N} \sum_{γ = γ_{1}}^{γ_{2}} p (λ_{i} = γ | τ = d, Y, Ω^{s}) {(y_{t_{i}} - ϕ_{d, t_{i} - γ}^{T} θ)}^{2} \end{matrix}

(28)

Similarly, the estimate of the model parameters $θ$ is derived by solving $\partial J (θ, δ) / \partial θ = 0$ and given as

\begin{matrix} θ^{s + 1} = {\sum_{d = d_{1}}^{d_{2}} p (τ = d | Y, Ω^{s}) \sum_{i = 1}^{N} \sum_{γ = γ_{1}}^{γ_{2}} p (λ_{i} = γ | τ = d, \\ {Y, Ω^{s}) ϕ_{d, t_{i} - γ} ϕ_{d, t_{i} - γ}^{T}}}^{- 1} {\sum_{d = d_{1}}^{d_{2}} p (τ = d | Y, Ω^{s}) \\ \sum_{i = 1}^{N} \sum_{γ = γ_{1}}^{γ_{2}} p (λ_{i} = γ | τ = d, Y, Ω^{s}) ϕ_{d, t_{i} - γ} y_{t_{i}}} . \end{matrix}

(29)

The estimate of time-delay $τ$ is obtained by solving the discrete optimization problem

τ^{s + 1} = \arg max_{d \in [d_{1}, d_{2}]} p (τ = d | Y, Ω^{s})

(30)

The proposed parameter estimation algorithm for a nonlinear multirate time-delay system described by the LPV FIR multirate time-delay model is summarized in Algorithm 1.

Algorithm 1 Parameter estimation for a nonlinear multirate time-delay system with uncertain output delays.
Require: The collected data $Y = {y_{t_{1} : t_{N}}, u_{1 : L}, z_{1 : L}}$ Ensure: The optimal parameter estimates $Ω^{} = {θ^{}, δ^{}}$ and delay estimate $τ^{}$ 1: Perform parameter initialization $Ω^{1}$ and let $s = 1$ . 2: repeat 3: Calculate the posterior probability $p (τ = d \| Y, Ω^{s})$ for $d \in [d_{1}, d_{2}]$ according to equation (21); Calculate the posterior probability $p (λ_{i} = γ \| τ = d, Y, Ω^{s})$ for $γ \in [γ_{1}, γ_{2}]$ and $d \in [d_{1}, d_{2}]$ according to equation (22). 4: Update the parameter estimates $θ$ and variance $δ$ according to equations (29) and (28), respectively. Update the delay estimate $τ^{s + 1}$ by solving discrete optimization problem (30). 5: Set $s = s + 1$ 6: until convergence 7: The obtained optimal estimates are denoted as $θ^{}$ , $δ^{}$ , and $τ^{*}$ , respectively.

Algorithm 1 Parameter estimation for a nonlinear multirate time-delay system with uncertain output delays.

Require: The collected data

Y = {y_{t_{1} : t_{N}}, u_{1 : L}, z_{1 : L}}

Ensure: The optimal parameter estimates

Ω^{*} = {θ^{*}, δ^{*}}

and delay estimate

τ^{*}

1: Perform parameter initialization

Ω^{1}

and let

s = 1

.
2: repeat
3: Calculate the posterior probability

p (τ = d | Y, Ω^{s})

for

d \in [d_{1}, d_{2}]

according to equation (21); Calculate the posterior probability

p (λ_{i} = γ | τ = d, Y, Ω^{s})

for

γ \in [γ_{1}, γ_{2}]

and

d \in [d_{1}, d_{2}]

according to equation (22).
4: Update the parameter estimates

θ

and variance

δ

according to equations (29) and (28), respectively. Update the delay estimate

τ^{s + 1}

by solving discrete optimization problem (30).
5: Set

s = s + 1

6: until convergence
7: The obtained optimal estimates are denoted as

θ^{*}

δ^{*}

, and

τ^{*}

, respectively.

Verification

Numerical example

Consider a nonlinear multirate time-delay system described by the following LPV FIR multirate time-delay model

\begin{matrix} x_{k} = c_{1} (z_{k}) u_{t - τ - 1} + c_{2} (z_{k}) u_{t - τ - 2} + c_{3} (z_{k}) u_{t - τ - 3}, \\ y_{t_{i}} = x_{t_{i} - λ_{i}} + ε_{t_{i}} \end{matrix}

(31)

where

\begin{matrix} c_{1} (z_{k}) = 1.5328 - 1.3438 z_{k} + 0.1562 z_{k}^{2} \\ c_{2} (z_{k}) = 1.3813 - 1.8750 z_{k} + 0.6250 z_{k}^{2} \\ c_{3} (z_{k}) = 1.1969 - 2.0625 z_{k} + 0.9375 z_{k}^{2} \end{matrix}

(32)

The input signal is chosen as the uniformly distributed white noise $u_{k} = - 1 + 2 U (0, 1)$ and the scheduling variable $z_{t}$ varies as a periodic signal $z_{k} = 0.5 + 0.5 \sin (0.15 π k)$ . The sampling period of the input and scheduling variables is $Δ t$ and the sampling period of the output variable is $3 Δ t$ . The input delay $τ = 3$ and the output delay $λ_{i}$ follow a uniform distribution in the integer range $[0, 2]$ . The input and slow-rate output data for parameter estimation are shown in Figure 1.

Figure 1.

The input data and output data denoted as *.

The proposed joint parameter and time-delay estimation algorithm is applied to the training data set and the convergence condition is set to $| | θ^{s + 1} - θ^{s} | |_{2} / | | θ^{s} | |_{2} < 1 e - 5$ . In order to evaluate the accuracy of estimated parameters in each iteration, the relative parameter estimation error (RE) $(| | θ_{t r u e} - θ^{s} | |_{2} / | | θ_{t r u e} | |_{2}) * 100 %$ is calculated. The relative parameter estimation error versus iteration is shown in Figure 2. In this simulation, the time-delay $τ$ is set as the integer range $[0, 10]$ and the estimated time-delay versus iteration is presented in Figure 3. As shown in Figures 2 and 3, the estimated parameters converge to true parameters after the fourteenth iteration.

Figure 2.

The relative parameter estimation error versus iteration.

Figure 3.

The estimated time-delay versus iteration.

The comparison of the proposed algorithm with the prediction error method (PEM) in Zhao et al. (2012) is conducted. In order to apply the PEM method, the considered LPV FIR model is transformed into linear regression form, one cost function is constructed based on one-step-ahead prediction error, and the model parameter estimates are finally obtained by optimizing this cost function. Note that the PEM method cannot estimate unknown time-delays. Therefore, the proposed algorithm is compared with the PEM method under two conditions: the unknown time-delays in the output data are accurately known (denoted as PEM-Accurate) and the unknown time-delays in the output data are unknown and ignored in the estimation process (denoted as PEM-Ignore). The estimated parameters are given in Table 1. The relative parameter estimation error between the true parameters and the estimated parameters is also calculated. It can be seen from this table, that when the time-delays in the output data are accurately known, the performance of the PEM method is comparable with the proposed method. However, when the time-delays in the output data are unknown and ignored, the PEM method cannot provide satisfactory parameter estimates.

Table 1.

The estimated parameters of the proposed method and PEM method.

True	Proposed	PEM-Accurate	PEM-Ignore
1.5328	1.5376	1.5310	0.6489
−1.3438	−1.3524	−1.2844	−0.8516
0.1562	0.1467	0.0826	0.2868
1.3813	1.3487	1.3496	0.9095
−1.8750	−1.6870	−1.7062	−0.4707
0.6250	0.4746	0.4974	−0.3422
1.1969	1.1813	1.1797	1.1022
−2.0625	−1.9881	−1.9871	−1.0269
0.9375	0.8820	0.8756	0.2476
RE	6.4110%	6.2420%	58.8156%

The performance of the proposed algorithm is further verified through Monte Carlo simulations with 100 different noise sequences under signal-to-noise ratios (SNRs) 20 dB, 15 dB, 10 dB, 5 dB, 0 dB, and −5 dB. All the parameter estimates are recorded and the mean and standard deviation of these parameter estimates are calculated to evaluate their accuracy. The results are shown in Figure 4 and Table 2. It can be seen from these results, that the estimated parameters are close to the true parameters and the standard deviations of the estimated parameters increase with the decrease of SNR. The bias norm (BN) $| | θ_{true} - E (θ^{*}) | |_{2}$ and variance norm (VN) $| | E (θ^{*} - E (θ^{*})) | |_{2}$ are two performance indices used to quantitatively evaluate the estimated parameter accuracy. The time-delay estimation accuracy (TDA), which is actually the ratio of estimated time-delay equals the true time-delay, is also calculated. The results of BNs, VNs, and TDAs are summarized in Table 3. It can be seen from this table, at the SNRs of 5 dB and higher, the BNs and VNs are close to zero indicating the parameter estimates converge to the true parameters and the time-delay $τ$ is correctly estimated. The performance of the proposed algorithm deteriorates with the decline of the SNR.

Figure 4.

The mean and standard deviation of parameter estimates under different SNRs.

Table 2.

The mean and standard deviation (std) of parameter estimates under SNR=15 dB and SNR=5 dB.

	15 dB		5 dB
True	Mean	Std	Mean	Std
1.5328	1.5362	0.0258	1.5148	0.1043
−1.3438	−1.3575	0.1200	−1.3472	0.4822
0.1562	0.1611	0.1213	0.1644	0.4853
1.3813	1.3880	0.0282	1.3948	0.1364
−1.8750	−1.8883	0.1465	−2.0096	0.6582
0.6250	0.6314	0.1464	0.7530	0.6320
1.1969	1.2010	0.0293	1.2194	0.1153
−2.0625	−2.0344	0.1499	−2.1295	0.6454
0.9375	0.8918	0.1528	1.0073	0.6220

Table 3.

The BNs, VNs, and TDAs of estimated parameters from Monte Carlo simulations.

	BN	VN	TDA
20 dB	0.0367	0.0158	100%
15 dB	0.0581	0.0490	100%
10 dB	0.1054	0.2246	100%
5 dB	0.2120	0.8836	100%
0 dB	0.3116	3.7346	88%
−5 dB	0.8939	13.4474	48%

The continuous stirred tank reactor

The continuous stirred tank reactor (CSTR) is a commonly used reactor type in chemical and biological engineering. The reactants are pumped into the reactor, a zero-order exothermic reaction happens in this reactor, the temperature of the reactor is controlled by the coolant flow, and the products are finally pumped out of the reactor. The process dynamic of the CSTR is described by Gopaluni (2008) as

\begin{matrix} \frac{d C_{A} (t)}{dt} = \frac{q (t)}{V} (C_{A 0} (t) - C_{A} (t)) - k_{0} C_{A} (t) \exp (\frac{- E}{RT (t)}) \\ \frac{dT (t)}{dt} = - \frac{(Δ H) k_{0} C_{A} (t)}{ρ C_{p}} \exp (\frac{- E}{RT (t)}) + \frac{q (t)}{V} (T_{0} (t) \\ - T (t)) + \frac{ρ_{c} C_{pc}}{ρ C_{p} V} q_{c} (t) {1 - \exp (\frac{- hA}{q_{c} (t) ρ C_{p}})} \\ \times (T_{c_{0}} (t) - T (t)) \end{matrix}

(33)

where the physical explanations and steady values of these process variables are also given in Gopaluni (2008).

The $C_{A} (t)$ and the $q (t)$ are chosen as the output variable and input variable, respectively. Since the dynamics between the $C_{A} (t)$ and the $q (t)$ of the CSTR are affected by the coolant flow rate, the $q_{c} (t)$ is selected as the scheduling variable. The input signal is a random binary signal and the scheduling variable signal is $q_{c} (t) = \sin (0.02 π t) + 100$ . The $q_{c} (t)$ and the $q (t)$ are sampled in fast-rate and the $C_{A} (t)$ is sampled with four multiples of fast-rate sampling period with uncertain delays in the integer range $[0, 3]$ . The fast-rate input and slow-rate output data are shown in Figure 5 and the recorded output data are denoted as *. The proposed algorithm is applied to identify a LPV FIR model with orders $n = 20$ and $M = 1$ for the CSTR system, and the comparison of the real fast-rate output and the simulated output of the identified model is given in Figure 6. A new validation data set is used to verify the accuracy of the identified model and the result is presented in Figure 7. The fitting accuracy of the real output and model simulated output is evaluated by calculating the percentage variance accounted for (VAF) defined as $(1 - V a r (y_{r e a l} - y^{*}) / V a r (y_{r e a l})) * 100 %$ . The calculated VAFs for the identification data set and validation data set are $82.28 %$ and $75.68 %$ , respectively. It can be concluded that the proposed algorithm is able to handle the multirate data and uncertain output delays, and identify a LPV FIR model that captures the dynamics of the CSTR system well.

Figure 5.

The input and output data of the CSTR system.

Figure 6.

The comparison of the real output and model simulated output for the identification data set.

Figure 7.

The comparison of the real output and model simulated output for the validation data set.

Conclusion

In this paper, the joint parameter and time-delay estimation problems for a nonlinear multirate time-delay system based on multirate data with output data corrupted by uncertain delays are considered. The LPV FIR multirate time-delay model is introduced and the iterative algorithm to simultaneously estimate the unknown model parameters and time-delays is proposed. The proposed algorithm can effectively handle the multirate data, the time-delay system property, and the uncertain output delays in the estimation process, providing satisfactory parameter estimates. Future research topics include the extensions of the proposed algorithm to an LPV system identification with random input delay and an LPV output error (OE) model identification.

Footnotes

Conflict of interest

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is supported in part by National Natural Science Foundation of China under Grant No. 61503097, the First-class General Financial Grant from the China Postdoctoral Science Foundation (Grant No. 2015M80264), and Financial Assistance under Heilongjiang Postdoctoral Fund (Grant No. LBH-Z15080).

References

Gopaluni

(2008) A particle filter approach to identification of nonlinear processes under missing observations. Canadian Journal of Chemical Engineering 86(6): 1081–1092.

Gopaluni

(2010) Nonlinear system identification under missing observations: The case of unknown model structure. Journal of Process Control 20(3): 314–324.

Laurain

Gilson

Toth

et al . (2009) Refined instrumental variable methods for identification of LPV Box–Jenkins models. Automatica 46(6): 959-967.

Masmoudi

Bellili

Affes

et al . (2011) A non-data-aided maximum likelihood time-delay estimator using importance sampling. IEEE Transactions on Signal Processing 59(10): 4505–4515.

McLachlan

Krishnan

(2007) The EM Algorithm and Extensions, 2nd edition.Hoboken: John Wiley & Sons.

Sun

Gao

Kaynak

(2011) Finite frequency H

\infty

control for vehicle active suspension systems. IEEE Transactions on Control Systems Technology 19(2): 416–422.

Sun

Gao

Kaynak

(2015) Vibration isolation for active suspensions with performance constraints and actuator saturation. IEEE Transactions on Mechatronics 20(2): 675–683.

Sun

Zhao

Gao

(2013) Saturated adaptive robust control for active suspension systems. IEEE Transactions on Industrial Electronics 60(9): 3889–3896.

Toth

Laurain

Gilson

et al . (2012) Instrumental variable scheme for closed-loop LPV model identification. Automatica 48(9): 2314–2320.

10.

van Wingerden

Verhaegen

(2009) Subspace identification of bilinear and LPV systems for open- and closed-loop data. Automatica 45(2): 372–381.

11.

Weyer (2001) System identification of an open water channel. Control Engineering Practice 9(12): 1289–1299.

12.

Wills

Ninness

(2011) System identification of linear parameter varying state-space models. In: Linear Parameter Varying System Identification: New Developments and Trends, World Scientific Publishing, Singapore, pp. 295–316.

13.

Xie

Yang

Huang

(2013) FIR model identification of multirate processes with random delays using EM algorithm. AIChE Journal 59(11): 4124–4132.

14.

Yan

Qian

Zhang

et al . (2016) H

\infty

fault detection for networked mechanical spring–mass systems with incomplete information. IEEE Transactions on Industrial Electronics 63(9): 5622–5631.

15.

Yang

Yan

(2015) Robust global identification of linear parameter varying systems with generalized expectation-maximization algorithm. IET Control Theory & Applications 9(7): 1103–1110.

16.

Yang

Yin

(2017) Variational Bayesian inference for FIR models with randomly missing measurements. IEEE Transactions on Industrial Electronics 64(5): 4217–4225.

17.

Yao

Jiao

(2014a) Extended-state-observer-based output feedback nonlinear robust control of hydraulic systems with backstepping. IEEE Transactions on Industrial Electronics 61(11): 6285–6293.

18.

Yao

Jiao

et al . (2014b) High-accuracy tracking control of hydraulic rotary actuators with modelling uncertainties. IEEE/ASME Transactions on Mechatronics 19(2): 633–641.

19.

Yin

Qiu

Gao

et al . (2017) Descriptor reduced-order sliding mode observers design for switched systems with sensor and actuator faults. Automatica 76: 282–292.

20.

Yin

Zhu

Qiu

et al . (2016) State estimation in nonlinear system using sequential evolutionary filter. IEEE Transactions on Industrial Electronics 63(6): 3786–3794.

21.

Zhang

Zhu

Zheng

(2015) Energy-to-peak state estimation for Markov jump RNNs with time-varying delays via nonsynchronous filter with nonstationary mode transitions. IEEE Transactions on Neural Networks and Learning Systems 26(10): 2346–2356.

22.

Zhao

Shmaliy

Huang

et al . (2015) Minimum variance unbiased FIR filter for discrete time-variant systems. Automatica 53(3): 355–361.

23.

Zhao

Huang

et al . (2012) Prediction error method for identification of LPV models. Journal of Process Control 22: 180–192.

24.

Zhong

Ding

(2010) Optimal fault detection for linear discrete time-varying systems. Automatica 46(8): 1395–1400.

25.

Zhong

Song

Ding

(2015) Parity space-based fault detection for linear discrete time-varying systems with unknown input. Automatica 59: 120–126.

26.

Zhu

Zhang

Zheng

(2016) Distributed

H_{\infty}

filtering for a class of discrete-time Markov jump Lur’e systems with redundant channels. IEEE Transactions on Industrial Electronics 63(3): 1876–1885.