# Identification of LTI Time-Delay Systems with Missing Output Data Using GEM Algorithm

Authors: Xianqiang Yang; Hamid Reza Karimi
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper considers the parameter estimation for linear time-invariant (LTI) systems in an input-output setting with output error (OE) time-delay model structure. The problem of missing data is commonly experienced in industry due to irregular sampling, sensor failure, data deletion in data preprocessing, network transmission fault, and so forth; to deal with the identification of LTI systems with time-delay in incomplete-data problem, the generalized expectation-maximization (GEM) algorithm is adopted to estimate the model parameters and the time-delay simultaneously. Numerical examples are provided to demonstrate the effectiveness of the proposed method.


---

## Body

## 1. Introduction

The advanced process control theories have enjoyed rapid development in the past several decades to meet the growing demands of closed-loop system performances, such as improved process safety and efficiency of plant operation, consistent product quality, and economic optimization [1]. These control strategies have improved the process automation and stability through providing control solutions for the process operated under the abnormal working conditions, such as process fault [2–5], network transmission delay [6], data packet dropouts [7], and modeling error. Generally, the implementation of these control strategies relies on the understanding of the process dynamics and the availability of an accurate mathematical model of the process. In view of the difficulties and complexities imposed by modeling using first principle method, the data-driven modeling method, in which the process model is retrieved from the process data, has become a main modeling method.Typically, the process data used in process modeling are generated by performing an identification experiment, in which a testing signal is designed and utilized to excite the process. Most of the conventional parameter estimation methods, such as prediction error method (PEM), instrumental variable (IV) method, and subspace method, assume that the identification data are sampled regularly and recorded properly. However, this is not always true in practical industry. For example, in the development of an inferential model for the sulfur content in the gas oil product, the sulfur concentration cannot be measured directly and the lab analysis is required which takes a long time. The process variable can be sampled in every minute, but the sulfur concentration is only available in every twelve hours. Another example is the industrial process with data transmission through the network. The recorded process data are corrupted by many network-induced problems, such as transmission delay and packet dropout or missing. Therefore, parameter estimation with irregular data has not been extensively investigated in the literature.Time-delays are commonly encountered in various engineering systems, such as chemical processes, mechanical systems, network control systems, transmission line, and economic systems [8]. Since the existence of time-delay usually causes performance degradation of the inferential model and is frequently a source of instability of the closed-loop system, it should be handled carefully in the modeling process. Common methods to estimate the time-delays are nonparametric methods (e.g., step test or correlation analysis) and grid searching method. For example, Wang and Zhang [9] considered the robust identification problems of linear continuous time-delay systems from step responses. A linear regression equation was derived from the solution of the output time response and its various-order integrals and solved by using IV-least squares (LS) method. The parameters of the transfer function were then recovered from the LS solution. Weyer [10] considered to build a model for the open water channel. The model parameters were estimated using the grid searching method in which one model was established for each time-delay in a range. The final model was selected as the model which gave the best prediction performance for a validation data set. The time-delay estimation methods mentioned above are to estimate the time-delay and model parameters in a separate way.Missing data problem is very common in process industry. A special example is the irregularly sampling system. Many critical parameters, such as the product concentration, steam quality, and90
% boiling point, cannot be measured directly by using the sensors. These parameters are measured through lab analysis, so only the slow rate data are available. However, the process variables, such as the temperature, pressure, and flow rate, can be measured on-line in fast rate by using the sensors. Therefore, we can treat the data samples between the slow rate data as missing data. Another example is the network control system in which data transmit via the wireless network or internet. Data missing occurs due to data packet dropout or missing. Other reasons for missing data are sensor fault, data recording system malfunction, and so forth. Several methods have been reported in the literature to handle missing data problem in system identification. For example, F. Ding and J. Ding [11] proposed an auxiliary model-based approach to cope with the problems of parameter estimation and output estimation with irregularly missing output data using the PEM method. The outputs of the auxiliary model were used in the identification process. Zhu et al. [12] considered the identification of systems with slowly and irregularly sampled output data. The output error method was employed to estimate the fast rate model based on the fast input and slow output data. However, the methods mentioned above just used part of the process data, which may lead to information missing. Moreover, the statistical properties of the model parameters and the process noise cannot be given in these methods.The work introduced in this paper aims at handling the identification problem of the LTI systems with missing output data in the presence of time-delay. The identification problem is formulated under the scheme of the generalized expectation-maximization (GEM) algorithm and the time-delay and missing output data are handled simultaneously. The GEM algorithm consists of expectation step (E-step) and maximization step (M-step). In the M-step, the maximization problem is transformed into an equivalent minimization problem and this problem is solved by using a general numerical optimization algorithm.The rest of this paper is organized as follows. The problem statement is presented in Section2. A brief revisit of the GEM algorithm and the mathematical formulation of the identification of LTI time-delay systems with incomplete data set are given in Section 3. Numerical examples are presented in Section 4 to show the effectiveness of the proposed method. The conclusions are given in Section 5.


## 2. Problem Statement

Consider the LTI system described by the following output error (OE) time-delay model:(1)
y
t
=
G
(
z
-
1
)
u
t
-
τ
+
e
t
,
where τ is the time-delay which is assumed to be integer multiples of the sampling period, e
t is the Gaussian white noise with zero mean and variance σ
2, and y
t and u
t are the output and input, respectively. The transfer function G
(
z
-
1
) has the following form:
(2)
G
(
z
-
1
)
=
∑
i
=
1
n
b
b
i
z
-
i
1
+
∑
i
=
1
n
a
a
i
z
-
i
.
Here, we assume that the model orders n
a and n
b are known a priori and the time-delay τ is uniformly distributed in a known range of [
d
1
,
d
2
].The identification data{
y
t
,
u
t
}
t
=
1,2
,
…
,
N are collected. We denote {
y
t
}
t
=
1,2
,
…
,
N as Y and {
u
t
}
t
=
1,2
,
…
,
N as U. Since part of the output data are missing completely at random (MCAR), the output data set Y can be divided into Y
obs
=
{
y
t
}
t
=
t
1
,
…
,
t
α and Y
mis
=
{
y
t
}
t
=
m
1
,
…
,
m
β. Therefore, the identification problem is to estimate the parameters θ
=
{
a
1
,
…
,
a
n
a
,
b
1
,
…
,
b
n
b
}, the noise variance σ
2, and the time-delay τ based on the identification data Y
obs and U.


## 3. Parameter Estimation Using the GEM Algorithm

### 3.1. GEM Algorithm Revisit

The GEM algorithm is a general-purpose iterative optimization algorithm to derive the maximum likelihood (ML) estimate and it has attracted great attentions of the researcher due to its flexibility in handling the missing data or hidden state [13]. Denote the missing data set by C
mis and the observed data set by C
obs. The main idea of the GEM algorithm is that, instead of optimizing the likelihood of the observed Y
obs, the conditional expectation of the complete data likelihood function with respect to the missing data set is calculated in the E-step and the maximization problem is solved in the M-step. The procedures of the GEM algorithm to calculate the ML estimate can be described as follows [13]:E-step: given theC
obs and the parameter estimate Θ
s in previous iteration, the Q-function can be calculated by
(3)
Q
(
Θ
∣
Θ
s
)
=
E
C
mis
∣
C
obs
,
Θ
s
{
log
⁡
p
(
C
mis
,
C
obs
∣
Θ
)
}
,M-step: find theΘ
s
+
1 to increase Q
(
Θ
∣
Θ
s
) over its value at Θ
s; that is,
(4)
Q
(
Θ
s
+
1
∣
Θ
s
)
⩾
Q
(
Θ
s
∣
Θ
s
)
.The E-step and M-step alternate until the relative change of the parameter estimate between neighboring iterations is smaller than a prespecified arbitrary small constant or the maximal iteration number is achieved.


### 3.2. LTI Time-Delay System Identification with Missing Output Using GEM Algorithm

Here, we treat the time-delayτ as a hidden state variable. The observed data set C
obs is constructed as C
obs
=
{
Y
obs
,
U
} and the missing data set C
mis is constructed as C
mis
=
{
Y
mis
,
τ
}. The parameter vector Θ is constructed as Θ
=
{
θ
,
σ
2
}.Based on the Bayesian property, the likelihood function of the complete data set can be decomposed into(5)
p
(
C
mis
,
C
obs
∣
Θ
)
=
p
(
Y
,
U
,
τ
∣
Θ
)
=
p
(
Y
∣
U
,
τ
,
Θ
)
p
(
τ
∣
U
,
Θ
)
p
(
U
∣
Θ
)
.The termp
(
Y
∣
U
,
τ
,
Θ
) can be further decomposed into
(6)
p
(
Y
∣
U
,
τ
,
Θ
)
=
p
(
y
N
,
…
,
y
1
∣
u
N
,
…
,
u
1
,
τ
,
Θ
)
=
p
(
y
N
∣
y
N
-
1
,
…
,
y
1
,
u
N
,
…
,
u
1
,
τ
,
Θ
)
×
p
(
y
N
-
1
,
…
,
y
1
∣
u
N
,
…
,
u
1
,
τ
,
Θ
)
=
∏
t
=
1
N
‍
p
(
y
t
∣
y
t
-
1
,
…
,
y
1
,
u
N
,
…
,
u
1
,
τ
,
Θ
)
.
Based on (1) and (2), y
t depends only on the previous input sequence u
t
-
1
:
1
=
{
u
t
-
1
,
…
,
u
1
}, the time-delay τ, and the parameter vector Θ. Therefore, (6) can be rewritten as
(7)
p
(
Y
∣
U
,
τ
,
Θ
)
=
∏
t
=
1
N
‍
p
(
y
t
∣
u
t
-
1
:
1
,
τ
,
Θ
)
.Since the time-delayτ is uniformly distributed in the range [
d
1
,
d
2
], the probability of τ taking any value in this range is a constant. Since the input U is measurable data and it is independent of the parameter vector Θ, the term p
(
U
∣
Θ
) is a constant. Therefore, the last two terms of (5) will not play a role in the following derivations. The complete data likelihood function can be further written as
(8)
p
(
Y
,
U
,
τ
∣
Θ
)
=
∏
t
=
1
N
‍
p
(
y
t
∣
u
t
-
1
:
1
,
τ
,
Θ
)
C
1
,
where C
1
=
p
(
τ
∣
U
,
Θ
)
p
(
U
∣
Θ
).Therefore, the conditional expectation of the log complete data densityQ
(
Θ
∣
Θ
s
) in (3) can be written as
(9)
Q
(
Θ
∣
Θ
s
)
=
E
Y
mis
,
τ
∣
C
obs
,
Θ
s
{
log
⁡
p
(
Y
,
U
,
τ
∣
Θ
)
}
=
E
Y
mis
,
τ
∣
C
obs
,
Θ
s
{
∑
t
=
1
N
‍
log
⁡
p
(
y
t
∣
u
t
-
1
:
1
,
τ
,
Θ
)
+
log
⁡
C
1
}
=
E
Y
mis
,
τ
∣
C
obs
,
Θ
s
{
∑
t
=
m
1
m
β
‍
log
⁡
p
(
y
t
∣
u
t
-
1
:
1
,
τ
,
Θ
)
1111111111111111
+
∑
t
=
t
1
t
α
‍
log
⁡
p
(
y
t
∣
u
t
-
1
:
1
,
τ
,
Θ
)
+
log
⁡
C
1
∑
t
=
m
1
m
β
‍
}
.The expectation is firstly taken with respect to the discrete variableτ; then we have
(10)
Q
(
Θ
∣
Θ
s
)
=
E
Y
mis
∣
C
obs
,
τ
,
Θ
s
{
∑
d
=
d
1
d
2
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
11111111111111111
×
∑
t
=
m
1
m
β
log
⁡
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
)
11111111111111111
+
∑
d
=
d
1
d
2
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
11111111111111111
×
∑
t
=
t
1
t
α
log
⁡
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
)
11111111111111111
+
log
⁡
C
1
∑
d
=
d
1
d
2
p
(
τ
=
d
C
o
b
s
,
Θ
s
)
}
.The expectation is then taken with respect to the continuous variableY
miss, so we have
(11)
Q
(
Θ
∣
Θ
s
)
=
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
∑
t
=
m
1
m
β
‍
∫
‍
p
(
y
t
∣
C
obs
,
τ
=
d
,
Θ
s
)
×
log
⁡
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
)
d
y
t
+
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
∑
t
=
t
1
t
α
‍
log
⁡
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
)
+
log
⁡
C
1
.In order to calculateQ
(
Θ
∣
Θ
s
), the unknown terms should be calculated firstly. Consider
(12)
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
=
p
(
τ
=
d
∣
Y
obs
,
U
,
Θ
s
)
=
p
(
Y
obs
∣
U
,
τ
=
d
,
Θ
s
)
p
(
τ
=
d
∣
U
,
Θ
s
)
∑
d
=
d
1
d
2
p
(
Y
obs
∣
U
,
τ
=
d
,
Θ
s
)
p
(
τ
=
d
∣
U
,
Θ
s
)
=
∏
t
=
t
1
t
α
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
s
)
p
(
τ
=
d
)
∑
d
=
d
1
d
2
∏
t
=
t
1
t
α
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
s
)
p
(
τ
=
d
)
=
∏
t
=
t
1
t
α
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
s
)
∑
d
=
d
1
d
2
∏
t
=
t
1
t
α
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
s
)
,
∫
‍
p
(
y
t
∣
C
obs
,
τ
=
d
,
Θ
s
)
log
⁡
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
)
d
y
t
=
∫
‍
p
(
y
t
∣
C
obs
,
τ
=
d
,
Θ
s
)
log
⁡
1
2
π
σ
2
×
exp
⁡
{
-
(
y
t
-
G
(
z
-
1
)
u
t
-
d
)
2
2
σ
2
}
d
y
t
=
-
1
2
log
⁡
(
2
π
σ
2
)
-
1
2
σ
2
∫
‍
p
(
y
t
∣
C
obs
,
τ
=
d
,
Θ
s
)
111111111111111111111
×
(
y
t
-
G
(
z
-
1
)
u
t
-
d
)
2
d
y
t
=
-
1
2
log
⁡
(
2
π
σ
2
)
-
1
2
σ
2
(
(
σ
s
)
2
+
(
G
s
(
z
-
1
)
u
t
-
d
)
2
)
+
1
σ
2
(
G
(
z
-
1
)
u
t
-
d
)
(
G
s
(
z
-
1
)
u
t
-
d
)
-
1
2
σ
2
(
G
(
z
-
1
)
u
t
-
d
)
2
=
-
1
2
log
⁡
(
2
π
σ
2
)
-
1
2
σ
2
(
(
σ
s
)
2
+
(
G
(
z
-
1
)
u
t
-
d
-
G
s
(
z
-
1
)
u
t
-
d
)
2
)
,
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
)
=
1
2
π
σ
2
exp
⁡
{
-
(
y
t
-
G
(
z
-
1
)
u
t
-
d
)
2
2
σ
2
}
.Therefore, theQ-function can be rewritten as
(13)
Q
(
Θ
∣
Θ
s
)
=
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
∑
t
=
m
1
m
β
‍
{
-
1
2
log
⁡
(
2
π
σ
2
)
-
1
2
σ
2
111111111
×
(
(
σ
s
)
2
+
(
G
(
z
-
1
)
u
t
-
d
-
G
s
(
z
-
1
)
u
t
-
d
)
2
)
-
1
2
log
⁡
(
2
π
σ
2
)
-
1
2
σ
2
(
(
σ
s
)
2
}
+
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
∑
t
=
t
1
t
α
‍
{
-
1
2
log
⁡
(
2
π
σ
2
)
-
1
2
σ
2
(
y
t
-
G
(
z
-
1
)
u
t
-
d
)
2
}
+
log
⁡
C
1
.In the M-step of the GEM algorithm, the unknown parameters should be estimated to increase theQ-function by solving an optimization problem. Taking the gradient of the Q-function (13) with respect to the σ
2 and setting it to zeros, we have
(14)
σ
^
2
=
1
N
{
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
∑
t
=
m
1
m
β
‍
(
(
G
(
z
-
1
)
u
t
-
d
-
G
s
(
z
-
1
)
u
t
-
d
)
2
+
(
σ
s
)
2
)
+
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
∑
t
=
t
1
t
α
‍
(
y
t
-
G
(
z
-
1
)
u
t
-
d
)
2
∑
d
=
d
1
d
2
‍
}
.
Substituting σ
^
2 into the Q-function (13), we get
(15)
Q
(
Θ
∣
Θ
s
)
=
-
N
2
log
⁡
(
2
π
)
-
N
2
-
N
2
log
⁡
{
1
N
[
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
1111111111111111
×
(
∑
t
=
m
1
m
β
‍
(
(
G
(
z
-
1
)
u
t
-
d
-
G
s
(
z
-
1
)
u
t
-
d
)
2
1111111111111111111111111
+
(
σ
s
)
2
)
11111111111111111111
+
∑
t
=
t
1
t
α
‍
(
y
t
-
G
(
z
-
1
)
u
t
-
d
)
2
∑
t
=
m
1
m
β
‍
(
(
G
(
z
-
1
)
u
t
-
d
-
G
s
(
z
-
1
)
u
t
-
d
)
2
)
∑
d
=
d
1
d
2
‍
]
1
N
[
∑
d
=
d
1
d
2
‍
p
}
.
Based on the monotonicity of the log function, the problem is transformed into optimizing the following cost function:
(16)
J
(
θ
∣
θ
s
)
=
1
N
[
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
(
∑
t
=
m
1
m
β
(
(
G
(
z
-
1
)
u
t
-
d
-
G
s
(
z
-
1
)
u
t
-
d
)
2
+
(
σ
s
)
2
)
‍
+
∑
t
=
t
1
t
α
‍
(
y
t
-
G
(
z
-
1
)
u
t
-
d
)
2
)

∑
d
=
d
1
d
2
‍
]
.Here, we introduce the variablex
t
τ denoting the noise-free output with time-delay τ. Based on (1) and (2), we have
(17)
x
t
τ
=
∑
i
=
1
n
b
b
i
z
-
i
1
+
∑
i
=
1
n
a
a
i
z
-
i
u
t
-
τ
=
(
ϕ
t
-
τ
τ
)
T
θ
,
where ϕ
t
-
τ
τ
=
[
-
x
t
-
1
τ
,
…
,
-
x
t
-
n
a
τ
,
u
t
-
τ
-
1
,
…
,
u
t
-
τ
-
n
b
]
T. Therefore, the cost function can be rewritten as
(18)
J
(
θ
∣
θ
s
)
=
1
N
[
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
(
∑
t
=
m
1
m
β
‍
(
(
(
ϕ
t
-
d
d
)
T
θ
-
(
x
t
d
)
s
)
2
+
(
σ
s
)
2
)
+
∑
t
=
t
1
t
α
‍
(
y
t
-
(
ϕ
t
-
d
d
)
T
θ
)
2
)

∑
d
=
d
1
d
2
‍
]
,
where (
x
t
d
)
s
=
(
ϕ
t
-
d
d
)
T
θ
s. However, the cost function (18) cannot be optimized directly due to the unmeasurable {
x
t
d
}
t
=
1
,
…
,
N
,
d
=
d
1
,
…
,
d
2. Here, we adopt the auxiliary model principle and the auxiliary model can be constructed based on the estimates obtained in the previous iteration. That is,
(19)
x
^
t
τ
=
∑
i
=
1
n
b
b
^
i
z
-
i
1
+
∑
i
=
1
n
a
a
^
i
z
-
i
u
t
-
τ
=
(
ϕ
^
t
-
τ
τ
)
T
θ
^
,
where ϕ
^
t
-
τ
τ
=
[
-
x
^
t
-
1
τ
,
…
,
-
x
^
t
-
n
a
τ
,
u
t
-
τ
-
1
,
…
,
u
t
-
τ
-
n
b
]
T. Therefore, the cost function (18) with ϕ
t
-
τ
τ substituted by ϕ
^
t
-
τ
τ can be optimized by using the damped Newton algorithm,
(20)
θ
s
+
1
=
θ
s
-
[
R
s
]
-
1
J
′
(
θ
s
∣
θ
s
)
,
where
(21)
J
′
(
θ
s
∣
θ
s
)
=
2
N
[
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
(
∑
t
=
m
1
m
β
‍
ϕ
^
t
-
d
d
(
(
ϕ
^
t
-
d
d
)
T
θ
s
-
(
x
t
d
)
s
)
1111111111
+
∑
t
=
t
1
t
α
‍
ϕ
^
t
-
d
d
(
(
ϕ
^
t
-
d
d
)
T
θ
s
-
y
t
)
)
]
,
R
s
=
2
N
[
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
∑
t
=
1
N
‍
ϕ
^
t
-
d
d
(
ϕ
^
t
-
d
d
)
T
]
+
λ
I
,
where λ is a constant.The time-delayτ can be selected as the delay in the range with maximal posterior probability. That is,
(22)
τ
^
=
argmax
d
⁡
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
.
The E-step and M-step alternate until the convergence condition of the GEM algorithm is met.


## 3.1. GEM Algorithm Revisit

The GEM algorithm is a general-purpose iterative optimization algorithm to derive the maximum likelihood (ML) estimate and it has attracted great attentions of the researcher due to its flexibility in handling the missing data or hidden state [13]. Denote the missing data set by C
mis and the observed data set by C
obs. The main idea of the GEM algorithm is that, instead of optimizing the likelihood of the observed Y
obs, the conditional expectation of the complete data likelihood function with respect to the missing data set is calculated in the E-step and the maximization problem is solved in the M-step. The procedures of the GEM algorithm to calculate the ML estimate can be described as follows [13]:E-step: given theC
obs and the parameter estimate Θ
s in previous iteration, the Q-function can be calculated by
(3)
Q
(
Θ
∣
Θ
s
)
=
E
C
mis
∣
C
obs
,
Θ
s
{
log
⁡
p
(
C
mis
,
C
obs
∣
Θ
)
}
,M-step: find theΘ
s
+
1 to increase Q
(
Θ
∣
Θ
s
) over its value at Θ
s; that is,
(4)
Q
(
Θ
s
+
1
∣
Θ
s
)
⩾
Q
(
Θ
s
∣
Θ
s
)
.The E-step and M-step alternate until the relative change of the parameter estimate between neighboring iterations is smaller than a prespecified arbitrary small constant or the maximal iteration number is achieved.


## 3.2. LTI Time-Delay System Identification with Missing Output Using GEM Algorithm

Here, we treat the time-delayτ as a hidden state variable. The observed data set C
obs is constructed as C
obs
=
{
Y
obs
,
U
} and the missing data set C
mis is constructed as C
mis
=
{
Y
mis
,
τ
}. The parameter vector Θ is constructed as Θ
=
{
θ
,
σ
2
}.Based on the Bayesian property, the likelihood function of the complete data set can be decomposed into(5)
p
(
C
mis
,
C
obs
∣
Θ
)
=
p
(
Y
,
U
,
τ
∣
Θ
)
=
p
(
Y
∣
U
,
τ
,
Θ
)
p
(
τ
∣
U
,
Θ
)
p
(
U
∣
Θ
)
.The termp
(
Y
∣
U
,
τ
,
Θ
) can be further decomposed into
(6)
p
(
Y
∣
U
,
τ
,
Θ
)
=
p
(
y
N
,
…
,
y
1
∣
u
N
,
…
,
u
1
,
τ
,
Θ
)
=
p
(
y
N
∣
y
N
-
1
,
…
,
y
1
,
u
N
,
…
,
u
1
,
τ
,
Θ
)
×
p
(
y
N
-
1
,
…
,
y
1
∣
u
N
,
…
,
u
1
,
τ
,
Θ
)
=
∏
t
=
1
N
‍
p
(
y
t
∣
y
t
-
1
,
…
,
y
1
,
u
N
,
…
,
u
1
,
τ
,
Θ
)
.
Based on (1) and (2), y
t depends only on the previous input sequence u
t
-
1
:
1
=
{
u
t
-
1
,
…
,
u
1
}, the time-delay τ, and the parameter vector Θ. Therefore, (6) can be rewritten as
(7)
p
(
Y
∣
U
,
τ
,
Θ
)
=
∏
t
=
1
N
‍
p
(
y
t
∣
u
t
-
1
:
1
,
τ
,
Θ
)
.Since the time-delayτ is uniformly distributed in the range [
d
1
,
d
2
], the probability of τ taking any value in this range is a constant. Since the input U is measurable data and it is independent of the parameter vector Θ, the term p
(
U
∣
Θ
) is a constant. Therefore, the last two terms of (5) will not play a role in the following derivations. The complete data likelihood function can be further written as
(8)
p
(
Y
,
U
,
τ
∣
Θ
)
=
∏
t
=
1
N
‍
p
(
y
t
∣
u
t
-
1
:
1
,
τ
,
Θ
)
C
1
,
where C
1
=
p
(
τ
∣
U
,
Θ
)
p
(
U
∣
Θ
).Therefore, the conditional expectation of the log complete data densityQ
(
Θ
∣
Θ
s
) in (3) can be written as
(9)
Q
(
Θ
∣
Θ
s
)
=
E
Y
mis
,
τ
∣
C
obs
,
Θ
s
{
log
⁡
p
(
Y
,
U
,
τ
∣
Θ
)
}
=
E
Y
mis
,
τ
∣
C
obs
,
Θ
s
{
∑
t
=
1
N
‍
log
⁡
p
(
y
t
∣
u
t
-
1
:
1
,
τ
,
Θ
)
+
log
⁡
C
1
}
=
E
Y
mis
,
τ
∣
C
obs
,
Θ
s
{
∑
t
=
m
1
m
β
‍
log
⁡
p
(
y
t
∣
u
t
-
1
:
1
,
τ
,
Θ
)
1111111111111111
+
∑
t
=
t
1
t
α
‍
log
⁡
p
(
y
t
∣
u
t
-
1
:
1
,
τ
,
Θ
)
+
log
⁡
C
1
∑
t
=
m
1
m
β
‍
}
.The expectation is firstly taken with respect to the discrete variableτ; then we have
(10)
Q
(
Θ
∣
Θ
s
)
=
E
Y
mis
∣
C
obs
,
τ
,
Θ
s
{
∑
d
=
d
1
d
2
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
11111111111111111
×
∑
t
=
m
1
m
β
log
⁡
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
)
11111111111111111
+
∑
d
=
d
1
d
2
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
11111111111111111
×
∑
t
=
t
1
t
α
log
⁡
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
)
11111111111111111
+
log
⁡
C
1
∑
d
=
d
1
d
2
p
(
τ
=
d
C
o
b
s
,
Θ
s
)
}
.The expectation is then taken with respect to the continuous variableY
miss, so we have
(11)
Q
(
Θ
∣
Θ
s
)
=
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
∑
t
=
m
1
m
β
‍
∫
‍
p
(
y
t
∣
C
obs
,
τ
=
d
,
Θ
s
)
×
log
⁡
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
)
d
y
t
+
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
∑
t
=
t
1
t
α
‍
log
⁡
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
)
+
log
⁡
C
1
.In order to calculateQ
(
Θ
∣
Θ
s
), the unknown terms should be calculated firstly. Consider
(12)
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
=
p
(
τ
=
d
∣
Y
obs
,
U
,
Θ
s
)
=
p
(
Y
obs
∣
U
,
τ
=
d
,
Θ
s
)
p
(
τ
=
d
∣
U
,
Θ
s
)
∑
d
=
d
1
d
2
p
(
Y
obs
∣
U
,
τ
=
d
,
Θ
s
)
p
(
τ
=
d
∣
U
,
Θ
s
)
=
∏
t
=
t
1
t
α
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
s
)
p
(
τ
=
d
)
∑
d
=
d
1
d
2
∏
t
=
t
1
t
α
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
s
)
p
(
τ
=
d
)
=
∏
t
=
t
1
t
α
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
s
)
∑
d
=
d
1
d
2
∏
t
=
t
1
t
α
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
s
)
,
∫
‍
p
(
y
t
∣
C
obs
,
τ
=
d
,
Θ
s
)
log
⁡
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
)
d
y
t
=
∫
‍
p
(
y
t
∣
C
obs
,
τ
=
d
,
Θ
s
)
log
⁡
1
2
π
σ
2
×
exp
⁡
{
-
(
y
t
-
G
(
z
-
1
)
u
t
-
d
)
2
2
σ
2
}
d
y
t
=
-
1
2
log
⁡
(
2
π
σ
2
)
-
1
2
σ
2
∫
‍
p
(
y
t
∣
C
obs
,
τ
=
d
,
Θ
s
)
111111111111111111111
×
(
y
t
-
G
(
z
-
1
)
u
t
-
d
)
2
d
y
t
=
-
1
2
log
⁡
(
2
π
σ
2
)
-
1
2
σ
2
(
(
σ
s
)
2
+
(
G
s
(
z
-
1
)
u
t
-
d
)
2
)
+
1
σ
2
(
G
(
z
-
1
)
u
t
-
d
)
(
G
s
(
z
-
1
)
u
t
-
d
)
-
1
2
σ
2
(
G
(
z
-
1
)
u
t
-
d
)
2
=
-
1
2
log
⁡
(
2
π
σ
2
)
-
1
2
σ
2
(
(
σ
s
)
2
+
(
G
(
z
-
1
)
u
t
-
d
-
G
s
(
z
-
1
)
u
t
-
d
)
2
)
,
p
(
y
t
∣
u
t
-
1
:
1
,
τ
=
d
,
Θ
)
=
1
2
π
σ
2
exp
⁡
{
-
(
y
t
-
G
(
z
-
1
)
u
t
-
d
)
2
2
σ
2
}
.Therefore, theQ-function can be rewritten as
(13)
Q
(
Θ
∣
Θ
s
)
=
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
∑
t
=
m
1
m
β
‍
{
-
1
2
log
⁡
(
2
π
σ
2
)
-
1
2
σ
2
111111111
×
(
(
σ
s
)
2
+
(
G
(
z
-
1
)
u
t
-
d
-
G
s
(
z
-
1
)
u
t
-
d
)
2
)
-
1
2
log
⁡
(
2
π
σ
2
)
-
1
2
σ
2
(
(
σ
s
)
2
}
+
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
∑
t
=
t
1
t
α
‍
{
-
1
2
log
⁡
(
2
π
σ
2
)
-
1
2
σ
2
(
y
t
-
G
(
z
-
1
)
u
t
-
d
)
2
}
+
log
⁡
C
1
.In the M-step of the GEM algorithm, the unknown parameters should be estimated to increase theQ-function by solving an optimization problem. Taking the gradient of the Q-function (13) with respect to the σ
2 and setting it to zeros, we have
(14)
σ
^
2
=
1
N
{
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
∑
t
=
m
1
m
β
‍
(
(
G
(
z
-
1
)
u
t
-
d
-
G
s
(
z
-
1
)
u
t
-
d
)
2
+
(
σ
s
)
2
)
+
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
∑
t
=
t
1
t
α
‍
(
y
t
-
G
(
z
-
1
)
u
t
-
d
)
2
∑
d
=
d
1
d
2
‍
}
.
Substituting σ
^
2 into the Q-function (13), we get
(15)
Q
(
Θ
∣
Θ
s
)
=
-
N
2
log
⁡
(
2
π
)
-
N
2
-
N
2
log
⁡
{
1
N
[
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
1111111111111111
×
(
∑
t
=
m
1
m
β
‍
(
(
G
(
z
-
1
)
u
t
-
d
-
G
s
(
z
-
1
)
u
t
-
d
)
2
1111111111111111111111111
+
(
σ
s
)
2
)
11111111111111111111
+
∑
t
=
t
1
t
α
‍
(
y
t
-
G
(
z
-
1
)
u
t
-
d
)
2
∑
t
=
m
1
m
β
‍
(
(
G
(
z
-
1
)
u
t
-
d
-
G
s
(
z
-
1
)
u
t
-
d
)
2
)
∑
d
=
d
1
d
2
‍
]
1
N
[
∑
d
=
d
1
d
2
‍
p
}
.
Based on the monotonicity of the log function, the problem is transformed into optimizing the following cost function:
(16)
J
(
θ
∣
θ
s
)
=
1
N
[
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
(
∑
t
=
m
1
m
β
(
(
G
(
z
-
1
)
u
t
-
d
-
G
s
(
z
-
1
)
u
t
-
d
)
2
+
(
σ
s
)
2
)
‍
+
∑
t
=
t
1
t
α
‍
(
y
t
-
G
(
z
-
1
)
u
t
-
d
)
2
)

∑
d
=
d
1
d
2
‍
]
.Here, we introduce the variablex
t
τ denoting the noise-free output with time-delay τ. Based on (1) and (2), we have
(17)
x
t
τ
=
∑
i
=
1
n
b
b
i
z
-
i
1
+
∑
i
=
1
n
a
a
i
z
-
i
u
t
-
τ
=
(
ϕ
t
-
τ
τ
)
T
θ
,
where ϕ
t
-
τ
τ
=
[
-
x
t
-
1
τ
,
…
,
-
x
t
-
n
a
τ
,
u
t
-
τ
-
1
,
…
,
u
t
-
τ
-
n
b
]
T. Therefore, the cost function can be rewritten as
(18)
J
(
θ
∣
θ
s
)
=
1
N
[
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
(
∑
t
=
m
1
m
β
‍
(
(
(
ϕ
t
-
d
d
)
T
θ
-
(
x
t
d
)
s
)
2
+
(
σ
s
)
2
)
+
∑
t
=
t
1
t
α
‍
(
y
t
-
(
ϕ
t
-
d
d
)
T
θ
)
2
)

∑
d
=
d
1
d
2
‍
]
,
where (
x
t
d
)
s
=
(
ϕ
t
-
d
d
)
T
θ
s. However, the cost function (18) cannot be optimized directly due to the unmeasurable {
x
t
d
}
t
=
1
,
…
,
N
,
d
=
d
1
,
…
,
d
2. Here, we adopt the auxiliary model principle and the auxiliary model can be constructed based on the estimates obtained in the previous iteration. That is,
(19)
x
^
t
τ
=
∑
i
=
1
n
b
b
^
i
z
-
i
1
+
∑
i
=
1
n
a
a
^
i
z
-
i
u
t
-
τ
=
(
ϕ
^
t
-
τ
τ
)
T
θ
^
,
where ϕ
^
t
-
τ
τ
=
[
-
x
^
t
-
1
τ
,
…
,
-
x
^
t
-
n
a
τ
,
u
t
-
τ
-
1
,
…
,
u
t
-
τ
-
n
b
]
T. Therefore, the cost function (18) with ϕ
t
-
τ
τ substituted by ϕ
^
t
-
τ
τ can be optimized by using the damped Newton algorithm,
(20)
θ
s
+
1
=
θ
s
-
[
R
s
]
-
1
J
′
(
θ
s
∣
θ
s
)
,
where
(21)
J
′
(
θ
s
∣
θ
s
)
=
2
N
[
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
×
(
∑
t
=
m
1
m
β
‍
ϕ
^
t
-
d
d
(
(
ϕ
^
t
-
d
d
)
T
θ
s
-
(
x
t
d
)
s
)
1111111111
+
∑
t
=
t
1
t
α
‍
ϕ
^
t
-
d
d
(
(
ϕ
^
t
-
d
d
)
T
θ
s
-
y
t
)
)
]
,
R
s
=
2
N
[
∑
d
=
d
1
d
2
‍
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
∑
t
=
1
N
‍
ϕ
^
t
-
d
d
(
ϕ
^
t
-
d
d
)
T
]
+
λ
I
,
where λ is a constant.The time-delayτ can be selected as the delay in the range with maximal posterior probability. That is,
(22)
τ
^
=
argmax
d
⁡
p
(
τ
=
d
∣
C
obs
,
Θ
s
)
.
The E-step and M-step alternate until the convergence condition of the GEM algorithm is met.


## 4. Simulation Examples

### 4.1. A Numerical Simulation Example

Consider the following LTI time-delay system described by the OE time-delay model:(23)
y
t
=
0.3
z
-
1
1
-
0.7
z
-
1
u
t
-
3
+
e
t
.
The input data and output data are generated by simulation and the noise e
t with zero mean and variance 0.01 is added to the output. The input and output data are shown in Figure 1. In the simulation, 12.5
% output data are randomly missing. The parameter range of the time-delay is set to [
1,5
]. The method proposed in this paper is used to estimate the parameters and the time-delay. The parameter estimate trajectories of the model parameters and the noise variance are shown in Figures 2 and 3, respectively. The estimated time-delay is 3 which is consistent with the true time-delay. To further verify the effectiveness of the proposed method, the simulations are also performed with 25
% output data missing and 50
% output data missing. The estimated parameters after 13 iterations are summarized in Table 1. It can be seen from these figures and the table that the proposed GEM algorithm has a good identification performance.Table 1
Estimated parameters after 13 iterations.
True value
a
=
-
0.7
b
=
0.3
τ
=
3
σ
2
=
0.01
Proportion of missing output
a
b
τ
σ
2
Full data set
−0.695
0.3026
3
0.0114
12.5
%
−0.694
0.3042
3
0.0116
25
%
−0.7012
0.2994
3
0.0121
50
%
−0.6979
0.3017
3
0.0124The input and output data.
(a)
(b)Figure 2
The estimated model parameters in each iteration.Figure 3
The estimated noise variance in each iteration.


### 4.2. The Continuous Stirred Tank Reactor

The Continuous Stirred Tank Reactor (CSTR) is a benchmark example used to test the performances of different modeling and control algorithms and the first principle model of the CSTR is described as [14]
(24)
d
C
A
d
t
=
q
V
(
C
A
0
-
C
A
)
-
k
0
C
A
exp
⁡
(
-
E
R
T
)
,
d
T
d
t
=
q
V
(
T
0
-
T
)
-
(
Δ
H
)
k
0
C
A
ρ
C
p
exp
⁡
(
-
E
R
T
)
+
ρ
c
C
p
c
ρ
C
p
V
q
c
{
1
-
exp
⁡
(
-
h
A
q
c
ρ
C
p
)
}
(
T
c
0
-
T
)
,
where the product concentration C
A and the temperature T are output variables and the coolant flow rate q
c is the input variable. The steady state values of the process variables can be found in Gopaluni [14]. In this simulation, the CSTR is operated at a steady state working point which is at q
c
=
97 L/min, C
A
=
0.0795mol/L, and T
=
443.4566K. The task here is to build a first-order model between C
A and q
c. The input and output data are generated through simulation and the noise with zero mean and variance 1
×
1
0
-
7 is added to the output data. Since time is needed to measure the concentration C
A, so the measurement delay with 1.5 minutes is also added to the output data. The input and output data is shown in Figure 4. In this simulation, 25
% output data are randomly missing and the parameter range of the time-delay is set to [
0.9,2.1
]. The proposed GEM algorithm is used to estimate the unknown parameters. The estimated parameters are a
^
=
-
0.468, b
^
=
0.0015, σ
^
2
=
1.5
×
1
0
-
7, and τ
^
=
1.5. The self-validation and the cross-validation results are shown in Figures 5 and 6. It can be seen from these results that the proposed method has a good identification performance and the estimated model can capture the dynamic behavior of the CSTR.The input and output data.
(a)
(b)Figure 5
The self-validation results. The blue line is the real process data and the red line is the simulated output of the estimated model.Figure 6
The cross-validation results. The blue line is the real process data and the red line is the simulated output of the estimated model.


## 4.1. A Numerical Simulation Example

Consider the following LTI time-delay system described by the OE time-delay model:(23)
y
t
=
0.3
z
-
1
1
-
0.7
z
-
1
u
t
-
3
+
e
t
.
The input data and output data are generated by simulation and the noise e
t with zero mean and variance 0.01 is added to the output. The input and output data are shown in Figure 1. In the simulation, 12.5
% output data are randomly missing. The parameter range of the time-delay is set to [
1,5
]. The method proposed in this paper is used to estimate the parameters and the time-delay. The parameter estimate trajectories of the model parameters and the noise variance are shown in Figures 2 and 3, respectively. The estimated time-delay is 3 which is consistent with the true time-delay. To further verify the effectiveness of the proposed method, the simulations are also performed with 25
% output data missing and 50
% output data missing. The estimated parameters after 13 iterations are summarized in Table 1. It can be seen from these figures and the table that the proposed GEM algorithm has a good identification performance.Table 1
Estimated parameters after 13 iterations.
True value
a
=
-
0.7
b
=
0.3
τ
=
3
σ
2
=
0.01
Proportion of missing output
a
b
τ
σ
2
Full data set
−0.695
0.3026
3
0.0114
12.5
%
−0.694
0.3042
3
0.0116
25
%
−0.7012
0.2994
3
0.0121
50
%
−0.6979
0.3017
3
0.0124The input and output data.
(a)
(b)Figure 2
The estimated model parameters in each iteration.Figure 3
The estimated noise variance in each iteration.


## 4.2. The Continuous Stirred Tank Reactor

The Continuous Stirred Tank Reactor (CSTR) is a benchmark example used to test the performances of different modeling and control algorithms and the first principle model of the CSTR is described as [14]
(24)
d
C
A
d
t
=
q
V
(
C
A
0
-
C
A
)
-
k
0
C
A
exp
⁡
(
-
E
R
T
)
,
d
T
d
t
=
q
V
(
T
0
-
T
)
-
(
Δ
H
)
k
0
C
A
ρ
C
p
exp
⁡
(
-
E
R
T
)
+
ρ
c
C
p
c
ρ
C
p
V
q
c
{
1
-
exp
⁡
(
-
h
A
q
c
ρ
C
p
)
}
(
T
c
0
-
T
)
,
where the product concentration C
A and the temperature T are output variables and the coolant flow rate q
c is the input variable. The steady state values of the process variables can be found in Gopaluni [14]. In this simulation, the CSTR is operated at a steady state working point which is at q
c
=
97 L/min, C
A
=
0.0795mol/L, and T
=
443.4566K. The task here is to build a first-order model between C
A and q
c. The input and output data are generated through simulation and the noise with zero mean and variance 1
×
1
0
-
7 is added to the output data. Since time is needed to measure the concentration C
A, so the measurement delay with 1.5 minutes is also added to the output data. The input and output data is shown in Figure 4. In this simulation, 25
% output data are randomly missing and the parameter range of the time-delay is set to [
0.9,2.1
]. The proposed GEM algorithm is used to estimate the unknown parameters. The estimated parameters are a
^
=
-
0.468, b
^
=
0.0015, σ
^
2
=
1.5
×
1
0
-
7, and τ
^
=
1.5. The self-validation and the cross-validation results are shown in Figures 5 and 6. It can be seen from these results that the proposed method has a good identification performance and the estimated model can capture the dynamic behavior of the CSTR.The input and output data.
(a)
(b)Figure 5
The self-validation results. The blue line is the real process data and the red line is the simulated output of the estimated model.Figure 6
The cross-validation results. The blue line is the real process data and the red line is the simulated output of the estimated model.


## 5. Conclusion

This paper considers the identification problem of LTI systems with irregular data set. The time-delay and the missing data are commonly encountered problems in process industry and the existence of these problems makes the process modeling a challenging task. The identification problem with incomplete data set in the presence of time-delay is formulated under the scheme of the GEM algorithm and the model parameters and the time-delay are estimated simultaneously in this algorithm. Numerical examples are presented to demonstrate the efficacy of the proposed method.



---

*Source: 242145-2014-03-17.xml*