JOURNAL OF INDUSTRIAL AND doi:10.3934/jimo.2013.9.471
MANAGEMENT OPTIMIZATION
Volume 9, Number 2, April 2013 pp. 471–486
A UNIFIED PARAMETER IDENTIFICATION METHOD FOR
NONLINEAR TIME-DELAY SYSTEMS
Qinqin Chai
a,b
Ryan Loxton
b
and Kok Lay Teo
b
Chunhua Yang
a
a
School of Information Science & Engineering, Central South University, Changsha, China
b
Department of Mathematics & Statistics, Curtin University, Perth, Australia
(Communicated by Cheng-Chew Lim)
Abstract. This paper deals with the problem of identifying unknown time-
delays and model parameters in a general nonlinear time-delay system. We
propose a unified computational approach that involves solving a dynamic op-
timization problem, whose cost function measures the discrepancy between
predicted and observed system output, to determine optimal values for the un-
known quantities. Our main contribution is to show that the partial derivatives
of this cost function can be computed by solving a set of auxiliary time-delay
systems. On this basis, the parameter identification problem can be solved
using existing gradient-based optimization techniques. We conclude the paper
with two numerical simulations.
1. Introduction. Time-delay systems arise in many important applications, in-
cluding medicine [21], chromatography [24], aerospace engineering [5], and chemical
reactor control [9]. Time-delays are often the cause of unpredictable and unusual
system behavior. For example, it is known that the introduction of time-delays can
destabilize systems that would otherwise be uniformly asymptotically stable [4, 19].
If the time-delays in a system are fixed and known, then existing powerful optimal
control algorithms (e.g. the control parameterization method; see [3, 25]) can be
applied to compute an optimal control for the system. In many systems, however,
the time-delays are unknown, which renders most of the existing optimal control
algorithms unusable. In this case, the time-delays must first be estimated before
an optimal control algorithm is applied. The problem of estimating the time-delays
(and possibly other unknown system parameters) from a given set of experimental
data is one of the key problems in the study of time-delay systems [20]. Such
problems are known as parameter identification problems.
Parameter identification for time-delay systems has been the subject of vigorous
research activity over the last two decades. An exact least squares algorithm for the
estimation of a single input delay is reported in [8]. Algebraic techniques [2] and
the steepest descent algorithm [6] have also been proposed for the identification of
input delays. In [26], information theory is used to identify time-delays for systems
in which each nonlinear term contains at most one unknown delay. Furthermore,
in [18] a genetic algorithm is developed for identifying a single time-delay in a
2010 Mathematics Subject Classification. Primary: 93B30; Secondary: 90C30, 90C90, 37N40.
Key words and phrases. Time-delay system, parameter identification, nonlinear optimization.
471
472 QINQIN CHAI, RYAN LOXTON, KOK LAY TEO AND CHUNHUA YANG
linear discrete-time system. Lyapunov function methods have also been used to
design delay estimators in [7]. One of the limitations of the existing identification
methods in [2, 6, 7, 8, 18, 26] is that they are mainly designed for systems with a
single input delay and no unknown parameters. Nonlinear systems with multiple
unknown delays are rarely considered in the literature.
In this paper, we consider a general nonlinear delay-differential system with un-
known time-delays and unknown system parameters. We formulate the problem of
identifying these unknown quantities as a nonlinear optimization problem in which
the cost function measures the least-squares error between predicted and observed
system output. This type of parameter identification problem was previously con-
sidered in [14] for systems in which each nonlinear component contains at most one
unknown delay and no unknown system parameters. However, in many real-world
systems, such as the purification process of zinc sulphate solution [22], the nonlinear
terms contain both delays and parameters that need to be identified. Our goal in
this paper is to extend the approach pioneered in [14] to these more complicated
systems. The key idea is to introduce a set of auxiliary delay-differential systems,
and then express the gradient of the least-squares cost function in terms of the so-
lution of these auxiliary systems. On this basis, numerical integration can be used
to solve the auxiliary systems, and thereby obtain the gradient of the cost function,
which is the main information needed to solve the parameter identification prob-
lem via numerical optimization techniques [16]. Based on this idea, we propose a
computational algorithm for identifying the unknown time-delays and system pa-
rameters in a general nonlinear system. We then demonstrate the effectiveness of
this algorithm on two nonlinear parameter identification problems, one being the
parameter identification problem for the zinc sulphate purification process.
2. Problem formulation. Consider the following nonlinear time-delay system:
˙
x(t) = f(t, x(t),
˜
x(t), ζ), t ∈ [0, T ], (1)
x(t) = φ(t), t ≤ 0, (2)
where T > 0 is a given terminal time; x(t) = [x
1
(t), . . . , x
n
(t)]
>
∈ R
n
is the state
vector ;
˜
x(t) = [x(t − τ
1
)
>
, . . . , x(t − τ
m
)
>
]
>
∈ R
nm
is the delayed state vector; and
ζ = [ζ
1
, . . . , ζ
r
]
>
∈ R
r
is a vector of unknown system parameters. Furthermore,
f : R × R
n
× R
nm
× R
r
→ R
n
and φ : R → R
n
are given functions.
The time-delays in (1)-(2) are unknown quantities that need to be determined.
We assume that the ith time-delay belongs to the interval [a
i
, b
i
], where a
i
and b
i
are given constants such that 0 ≤ a
i
< b
i
. Hence, the unknown time-delays satisfy
the following bound constraints:
a
i
≤ τ
i
≤ b
i
, i = 1, . . . , m. (3)
Any vector τ = [τ
1
, . . . , τ
m
]
>
∈ R
m
that satisfies (3) is called a candidate time-delay
vector. Let T denote the set of all such candidate time-delay vectors.
In addition to the time-delays, the system parameters in (1)-(2) are also unknown
quantities that need to be determined. We suppose that
c
j
≤ ζ
j
≤ d
j
, j = 1, . . . , r, (4)
where c
j
and d
j
are given real numbers such that 0 ≤ c
j
< d
j
. Note that there is
no loss of generality in assuming that c
j
≥ 0; if c
j
< 0, then we may replace ζ
j
with
ζ
j
+ c
j
. Any vector ζ = [ζ
1
, . . . , ζ
r
]
>
∈ R
r
that satisfies (4) is called a candidate
parameter vector. Let Z denote the set of all such candidate parameter vectors.
PARAMETER IDENTIFICATION FOR NONLINEAR TIME-DELAY SYSTEMS 473
The output of system (1)-(2) is given by
y(t) = g(x(t), ζ), t ∈ [0, T ], (5)
where g : R
n
× R
r
→ R
p
is a given function.
We assume that the following conditions are satisfied.
Assumption 1. The given functions f and g are continuously differentiable, and
φ is twice continuously differentiable.
Assumption 2. There exists a real number L
1
> 0 such that
|f(t, x,
˜
x, ζ)| ≤ L
1
(1 + |x| + |
˜
x| + |ζ|), (t, x,
˜
x, ζ) ∈ R × R
n
× R
nm
× R
r
,
where | · | denotes the Euclidean norm.
On the basis of Assumptions 1 and 2, the dynamic system (1)-(2) admits a unique
solution corresponding to each pair (τ , ζ) ∈ T × Z [1]. We denote this solution by
x(·|τ , ζ). Substituting x(·|τ , ζ) into (5) gives y(·|τ , ζ), the predicted system output
corresponding to (τ , ζ) ∈ T × Z. More formally,
y(t|τ , ζ) = g(x(t|τ , ζ), ζ), t ∈ [0, T]. (6)
Suppose that the output from system (1)-(2) has been measured experimentally at
times t = t
l
, l = 1, . . . , q, where each t
l
∈ [0, T ]. Let
ˆ
y
l
∈ R
p
denote the measured
output at time t = t
l
. Then the problem of identifying the unknown time-delays
and system parameters can be formulated mathematically as follows.
Problem (P). Choose τ ∈ T and ζ ∈ Z to minimize the following cost function:
J(τ , ζ) =
q
X
l=1
y(t
l
|τ , ζ) −
ˆ
y
l
2
. (7)
Problem (P) is a nonlinear dynamic optimization problem whose decision vari-
ables are the delays and model parameters in system (1)-(2). Our aim is to select
optimal values for these delays and parameters so that the predicted system output
best fits the experimental data. Almost all of the existing optimization techniques
for time-delay systems are based on the assumption that the delays are fixed and
known (see, for example, [3, 12, 25]). Problem (P) is unique in that the delays
are not fixed, but are instead decision variables to be chosen optimally. The cost
function in Problem (P) is also highly non-standard, as it depends on the system’s
state at a set of discrete time points, not just at the terminal time. Such cost
functions have been considered in [13, 17] for non-delay systems, and in [22] for
systems with fixed delays. However, the computational techniques developed in
references [13, 17, 22] are not applicable to Problem (P) because the time-delays in
system (1)-(2) are variable.
3. Gradient computation. Problem (P) involves choosing a finite number of
decision variables to minimize the cost function (7). Thus, in principle, Problem (P)
can be viewed as a nonlinear programming problem. Standard algorithms for solving
nonlinear programming problems—for example, sequential quadratic programming
or interior-point methods [16]—typically require the gradient of the cost function,
which is difficult to determine in Problem (P) because the delays and parameters
influence (7) implicitly through the dynamic system (1)-(2). The aim of this section
is to develop an efficient computational method for computing the gradient of the
cost function in Problem (P). This method, which is inspired by our earlier work in
474 QINQIN CHAI, RYAN LOXTON, KOK LAY TEO AND CHUNHUA YANG
[10, 11, 14, 15], can be integrated with a standard nonlinear programming algorithm
to solve Problem (P).
3.1. Preliminaries. Throughout this subsection, we assume that k ∈ {1, . . . , m}
and (τ , ζ) ∈ T × Z are arbitrary but fixed. For simplicity, we write x(t) instead of
x(t|τ , ζ), and x
(t) instead of x(t|τ + e
k
, ζ), where e
k
denotes the kth unit basis
vector in R
m
.
Define
I = [a
k
− τ
k
, b
k
− τ
k
].
Note that I 6= ∅ and 0 ∈ I. Clearly,
∈ I ⇐⇒ τ + e
k
∈ T .
For each ∈ I, define
ϕ
(t) = x
(t) − x(t), t ≤ T,
and
θ
,i
(t) = x
(t − τ
i
− δ
ki
) − x(t − τ
i
), t ≤ T, i = 1, . . . , m,
where δ
ki
denotes the Kronecker delta function. Furthermore, let
θ
(t) = [(θ
,1
(t))
>
, . . . , (θ
,m
(t))
>
]
>
∈ R
nm
, t ≤ T.
Clearly,
θ
,i
(t) = ϕ
(t − τ
i
), t ≤ T, i 6= k, (8)
ϕ
(t) = 0, t ≤ 0. (9)
In the sequel, we will use the notation
∂
∂
˜
x
i
to denote differentiation with respect to
the ith delayed state in
˜
x(t) (i.e. differentiation with respect to x(t − τ
i
)).
Now, define
χ(t) =
(
˙
φ(t), if t ≤ 0,
f(t, x(t),
˜
x(t), ζ), if t ∈ (0, T ].
We immediately see that, for almost all t ∈ (−∞, T],
˙
x(t) = χ(t).
Let
¯
b > 0 be a fixed constant such that
¯
b ≥ b
i
for each i = 1, . . . , m. By following
similar arguments to those used in [14], it is possible to show that there exists a
positive real number L
2
> 0 such that for each ∈ I,
|x
(t)|, |χ(t)| ≤ L
2
, t ∈ [−
¯
b, T ], (10)
and
|ϕ
(t)|, max
i=1,...,m
|θ
,i
(t)| ≤ L
2
||, t ∈ [0, T ]. (11)
Furthermore, one can also show that for almost all t ∈ [0, T ],
lim
→0
θ
,k
(t) − ϕ
(t − τ
k
)
= −χ(t − τ
k
). (12)
See Appendix B of [14] for more details.
PARAMETER IDENTIFICATION FOR NONLINEAR TIME-DELAY SYSTEMS 475
3.2. State variation with respect to time-delays. The solution of system
(1)-(2) is normally viewed as a function of time, with τ and ζ being fixed vectors.
Now, by instead fixing t ∈ (−∞, T ], while allowing the vectors τ and ζ to vary,
we obtain a new function x(t|·, ·) : T × Z → R
n
whose value at (τ , ζ) ∈ T × Z is
x(t|τ , ζ). In the following result, we show that x(t|·, ·) is differentiable with respect
to the time-delays. This result is central to the development of a computational
procedure for solving Problem (P).
Theorem 3.1. Let t ∈ (0, T ] be a fixed time point. Then x(t|·, ·) is differentiable
with respect to τ
k
on T × Z. In fact, for each (τ , ζ) ∈ T × Z,
∂x(t|τ , ζ)
∂τ
k
= Λ
k
(t|τ , ζ), k = 1, . . . , m, (13)
where Λ
k
(·|τ , ζ) satisfies the auxiliary time-delay system
˙
Λ
k
(t) =
∂f (t, x(t),
˜
x(t), ζ)
∂x
Λ
k
(t) +
m
X
i=1
∂f (t, x(t),
˜
x(t), ζ)
∂
˜
x
i
Λ
k
(t − τ
i
)
−
∂f (t, x(t),
˜
x(t), ζ)
∂
˜
x
k
χ(t − τ
k
)
(14)
with initial condition
Λ
k
(t) = 0, t ≤ 0. (15)
Proof. Let k ∈ {1, . . . , m} and (τ , ζ) ∈ T × Z be arbitrary but fixed. As in
Subsection 3.1, we write x
(t) instead of x(t|τ +e
k
, ζ), and x(t) instead of x(t|τ , ζ).
For each ∈ I \ {0}, define
ρ() =
Z
T
0
−1
θ
,k
(s) −
−1
ϕ
(s − τ
k
) + χ(s − τ
k
)
ds. (16)
It follows from (9), (10), and (11) that for each ∈ I \ {0},
−1
θ
,k
(s) −
−1
ϕ
(s − τ
k
) + χ(s − τ
k
)
≤ 3L
2
, s ∈ [0, T ].
Hence, the integrand in (16) is uniformly bounded with respect to ∈ I \ {0}.
Furthermore, it follows from (12) that the integrand in (16) converges to zero almost
everywhere on [0, T ] as → 0. Thus, from the Lebesgue dominated convergence
theorem,
lim
→0
ρ() = lim
→0
Z
T
0
−1
θ
,k
(s) −
−1
ϕ
(s − τ
k
) + χ(s − τ
k
)
ds = 0.
Now, keeping ∈ I \ {0} fixed for the time being, we define
¯
f(s, α) = f(s, x(s) + αϕ
(s),
˜
x(s) + αθ
(s), ζ), (s, α) ∈ [0, T ] × [0, 1].
Then by the chain rule,
∂
¯
f(s, α)
∂α
=
∂
¯
f(s, α)
∂x
ϕ
(s) +
m
X
i=1
∂
¯
f(s, α)
∂
˜
x
i
θ
,i
(s), (17)
where
∂
¯
f(s, α)
∂x
=
∂f (s, x(s) + αϕ
(s),
˜
x(s) + αθ
(s), ζ)
∂x
, (18)
∂
¯
f(s, α)
∂
˜
x
i
=
∂f (s, x(s) + αϕ
(s),
˜
x(s) + αθ
(s), ζ)
∂
˜
x
i
. (19)