1. Introduction

Dynamical systems with delays are important class of models describing phenomena appearing in many areas, for example in industry or biology. One of the practical problems related to such models is a need to estimate their parameters based on measurements carried out in the real system (process).

There are many works dealing with this problem [1], [4], [12], [14], [15], [16], [17], [18]. Unfortunately, most of proposed approaches assume that the analyzed system is linear and both input and output signals for delaying elements can be measured.

More general and universal approaches, for non-linear systems with delays have been proposed in [13] and [8]. Both approaches depends on gradient-based minimization of an appropriately defined objective function. The latter approach uses so-called structural adjoint sensitivity analysis, which decreases significantly computational effort when many parameters are estimated. Moreover, this approach is more general and can be applied to any dynamical system presented in structural form as block diagram containing any number of delay elements.

All above mentioned methods are focused only on estimation of time delays and eventually other parameters of the mathematical model. But general task of identification of dynamical systems requires also estimation of initial conditions in the situation when they are unknown.

In case of one discrete delay element the initial condition (its "state" for time $t=0$) is a function of time specified for an interval $[-\tau, 0]$, where $\tau$ is a delay time of this element.

There are relatively little works related to the problem of estimation of initial functions for systems with delays. In paper [2] a gradient based approach to estimation of initial functions for non-linear systems described by retarded type delay differential equations (RDDE). Another paper [3] deals with systems described by neutral type delay differential equations (NDDE). In both works a gradient-based estimation of initial functions is done by using adjoint sensitivity analysis.

In this work a more general structural adjoint sensitivity analysis is utilized. It is especially useful for models described by block diagrams and was originally developed for neural networks [5] and afterwards was used for different models described by ordinary differential equations [6,7,11], systems with delays [8] and age-structured models [10]. Recently it has been used for spatiotemporal models of tumor growth [9].

The structural sensitivity approach can be applied to any non-linear dynamical system presented as a block diagram and containing many discrete delay elements. Therefore, it may be used for wider class than analysed in [2] and [3]. For example it may be applied for systems containing delays in input (control) channel, which is not allowed in RDDE and NDDE models. The proposed approach can be used for systems which output signals are measurable continuously and for sampled systems where the information output signals is available only at discrete time moments.

---

2. Problem formulation

Let us consider a model of dynamical system with one isolated delay element presented in Fig. 1. We do not assume any particular structure of the model $M$, for example RDDE or NDDE, but we assume that it is given in structural form － as a block diagram containing basic elements such as:

1. Linear static element represented by a gain matrix $A$.

2. Linear continuous-time dynamical element represented by a transfer function $K(s)$.

3. Linear discrete-time dynamical element represented by a transfer function $K(z)$.

4. Non-linear static element described by a function $f(\cdot)$.

5. Summing junction.

6. Branching node.

7. Ideal d-c pulser, placed between discrete-time part of the system and the continuous-time part, which output signal contains Dirac pulses multiplied by instantaneous value of its discrete-time input.

8. Ideal c-d pulser, placed between continuous-time part of the system and the discrete-time part, which output signal contains Kronecker pulses multiplied by instantaneous value of its continuous-time input.

Using such a set of elements one can present any non-linear hybrid continuous-discrete dynamical system with delays of arbitrary structure as a block diagram.

For the sake of simplicity the system presented in Fig. 1 contains only one delay element but in general case there can be more delay elements with different delay times.

The delay element is described by the input-output relation

$r(t)=q(t-\tau)$

(1)

with the initial condition

$q(t)=\varphi(t) \;\;\;\;\textrm{for} \;\;\;\; t \in [-\tau,0]$

(2)

The function $\varphi(t)$ is called the initial function of the delay element.

The delay element can be also mathematically described in Laplace operator domain by its transfer function which is frequently used for example in control systems theory. The transfer function $K(s)$ is defined as a ratio of the output of a system to the input of a system, in the Laplace domain, under zero initial conditions. Taking into account the properties of the Laplace transform, it can be shown that the transfer function of the delay element has the form

$K(s)=\frac{R(s)}{Q(s)}=\frac{ \mathcal{L} \{r(t)\} }{\mathcal{L} \{q(t)\}}=e^{-s\tau}$

(3)

where $\mathcal{L} \{ \cdot \}$ stands for the Laplace transform.

We also assume that the output signal $d(t)$ of the real identified system, also referred to as plant, can be measured only at discrete time moments $t_1, t_2, \dots ,t_N \in [0,t_f]$ where $N$ is a number of measurements, and $t_\textrm{f}$ is a final time. These measurements will be denoted by $d(1),d(2),\dots,d(N)$, and corresponding instantaneous values of the output signal of the model by $y(1),y(2),\dots,y(N)$.

Let us define an objective function which is a measure of discrepancy between the measurements and the output of the model

$J=\frac{1}{2}\sum\limits_{n=1}^{N}(y(n)-d(n))^2$

(4)

Problem 1. Find the initial function of the delay element $\varphi(t)$ minimizing the objective function (4).

The above task will be solved iteratively using the gradient-based approach. Hence, we need to solve the folowing sub-problem

Problem 2. Find the gradient of the objective function (4) in the space of the initial function $\varphi(t)$.

To solve the Problem 2 we will use the adjoint sensitivity analysis. In works [5], [7] rules for construction on the sensitivity model and the adjoint model have been presented. In addition in [8] such rules has been extended to systems with delays and it has been shown how to perform the sensitivity analysis with respect to delay times. Now, we are going to show how to calculate the gradient of the objective function in the space of the initial function of the delay element

---

3. Model of the delay element, its sensitivity model and the adjoint model

Before we start to solve problems formulated in previous section let us present one delay element in the form which will be more suitable for further analysis. This form comes from the observation that the delay element with non-zero initial condition can be replaced by a delay element with zero initial conditions and with additional signal $\psi(t)$ additively introduced as presented in Fig. 2.

A function $\psi(t)$ is related to the initial function of the delay element $\varphi(t)$ by the following relation

$\psi (t) = \left\{ {\begin{array}{*{20}c} {\varphi (t - \tau )}&{{\rm{for}}}&{t \le \tau } \\ 0&{{\rm{for}}}&{t > \tau } \\ \end{array}} \right.$

(5)

Therefore the task of finding the gradient of the objective function in the initial function space $\varphi(t)$ for time interval $[-\tau,0]$ can be replaced by the following problem:

Problem 3. Find the gradient of the objective function (4) in the space of the input signal $\psi(t)$ for time interval $t \in [0,\tau]$.

The delay time $\tau$ has also been presented as an input "signal" of the delay element presented in Fig. 2. This can be utilized in the case when one looks also for the gradient (partial derivative) of the objective function with respect to the delay time.

The sensitivity model of the delay element presented in Fig. 2, which describes relationship between variations of all signals $\bar q(t)$, $\bar \psi(t)$, $\bar \tau$ and $\bar r(t)$ is presented in Fig. 3a. Since the input signal $\psi(t)$ enters additively the the model from Fig. 2, its variation $\bar \psi(t)$ enters in the same way the sensitivity model from Fig. 3a. The rest part of the sensitivity model has been developed and justified in previous work [8].

Rules for construction of the adjoint system presented in works [5] and [7] specify, among others, that the directions of all signals should be reversed and all summing junctions should be replaced by branching nodes. As a result we obtain the adjoint system of one delay element presented in Fig. 3b. The output signal $\beta(t)$ corresponds to the input signal $\psi(t)$ in the original model. It will be used (after reversing in time) as a solution to the Problem 3.

---

4. Problem solution

To solve the Problem 3 (and consecutively Problem 2 and Problem 2) let us extend the general model presented in Fig. 1. The extended model, presented in Fig. 4, takes into account that we minimize the objective function (4). It is obtained by using a non-linear element calculating the quadratic function in (4) and the discrete transfer function $z \over {z-1}$ realizing summing over time. Thanks to these extensions the additional output signal $\tilde J(n)$ has such a property that its final value is equal to the objective function: $\tilde J(N)=J$.

Moreover, the extended model has an additional input signal $\psi(t)$, which has been discussed in previous section. We will calculate the sensitivity of $\tilde J(N)$ with respect to this input signal. Since in this article we are not interested in finding the sensitivity of $J$ w.r.t. the delay time, the additional input signal $\tau$, presented previously in Fig. 2, is now omitted.

The extended model presented in Fig. 4 is an example of a hybrid continuous-discrete-time system. It contains both, continuous-time part (for time $t$) and discrete-time part (for discrete time moments $n$) and the interfacing c-d sampler. Rules for construction of the adjoint for such system were presented in our previous works: [5], [7]. Using them it is easy to construct the adjoint system, which is presented in Fig. 5. The non-stationary gain $e(N-n+1)$ resulted as a reversed in time derivative of the previous non-linear quadratic function in the extended model. The block denoted by $\widehat M$ is a system adjoint to the part $M$ of the original model and can be constructed based on its structure using the same rules.

The adjoint system stimulated by the Kronecker pulse $\delta(n)$ generates as an output the signal $\beta(t)$, which, after reversing in time, is the searched gradient of the objective function in the space of the input signal $\psi(t)$:

$\beta(t_f-t)=\nabla_{\psi(t)}J$

(6)

This signal in the time interval $[0,\tau]$ is a solution to the Problem 3. The same signal, shifted in time according to (5), is a solution to the Problem 2 and can be used during gradient-based optimization procedure, which gives an estimated solution to the Problem 1.

---

5. Numerical examples

To illustrate how the proposed approach works, we provide results of six numerical examples. They were performed under different conditions which are shown in Table 1. First of all, three different models, with different number of delays and their location, were used. Structures (block diagrams) of models A, B and C are presented in figures 6, 12 and 14 respectively. Moreover, in all examples times of discrete measurements $t_1, t_2, \dots ,t_N$ are equidistant but sampling time is different. Finally, we show cases where delay times are estimated in addition to estimation of initial functions.

Each model, A, B and C, is a first-order system, described by a first-order delay differential equation and hence contains one integrating element － transfer function $\frac{1}{s}$. In each numerical example it is assumed zero initial condition for the integrator and non-zero initial condition(s) for delay element(s) i.e. initial function(s).

Each model has one scalar external input signal $u(t)$ stimulating the system and one scalar output signal $y(t)$. In all numerical examples $u(t)$ is assumed to be the a step function i.e. it is constant $u(t)=1$ for $t\geq 0$.

In every numerical example measurements $d(1),d(2),\dots,d(N)$ are obtained by simulation the virtual plant, which has the same structure as the model¹. In all cases the final time of simulation $t_f=3~[\textrm{s}]$ and delay time(s) in the plant $\tau = \tau_1 = \tau_2 = 1~[\textrm{s}]$.

¹In fact, in presented numerical examples both plant and the model are "models" but we consistently use two different names to emphasize that the plant generates measurements and the model is fitted to measurements.

The gradient obtained by simulation of the adjoint model is used in the simplest iterative gradient-based optimization procedure:

$\psi^{k+1}(t)=\psi^{k}(t)-c \nabla_{\psi(t)}J$

(7)

where $k$ is an index if current iteration ad $c$ is a positive constant assuring convergence of the procedure. The $c$ parameter has been chosen separately for each example to speed up the estimation procedure and preserve its convergence. In examples where time delay $\tau$ is unknown, a similar updating rule is used:

$\tau^{k+1}=\tau^{k}-c \nabla_{\tau}J$

(8)

with the same value of $c$ that applied in (7).

Example 1. In the first example the model A is used. It is presented in Fig. 6.

It is a simple first-order system with one delay described by the following delay differential equation:

$\dot y(t)=-y(t)-y(t-\tau)+u(t)$

(9)

An adjont system for the Model A, created by using rules described in [5], [7], is presented in Fig. 7.

This is a part of the overall adjont system from Fig. 5 and generates a function $\beta(t)$ which is a reversed in time searched gradient according to (6).

The unknown (estimated) initial function $\psi(t)$ of the delay element applied in the plant is a stepwise function presented in Fig. 8a by a dashed line. In all examples we consequently present only secondary initial functions $\psi(t)$ associated with the primary initial function $\varphi(t)$ by the relation (5). Of course the original initial function $\varphi(t)$ has the same shape but is specified for shifted time interval $[-\tau,0]$.

In this example, and in the next one, it is assumed that measurements are quasi-continuous i.e. sampling time is infinitesimally small² $t_s \rightarrow 0^+$ and there is no effect of sampling.

²For the computer simulation it is the same as the variable step size (with assumed upper limit) used by ODE solver.

The results of the estimation of the initial function obtained after nearly 500 iterations of the gradient-descent optimization procedure (7) are presented in Fig. 8. The starting initial function $\psi^0(t)$ for the optimization procedure, in this and in the rest of examples, was chosen as a constant zero function. The estimate of the initial function $\psi(t)$ is presented in Fig. 8a － solid line. It can be observed that it differs from the true initial function in the plant － dashed line, especially around jumps of the true initial function. Nevertheless, the objective function reached a very low value, approx. $10^{-4}$ － see Fig. 8c － and can be less for longer optimization process. In general, one can see that the estimation process is convergent. The output of the model is very close to the output of the plant, see Fig. 8b where dashed line for the plant is nearly invisible. The absolute value of the prediction error is small as well － Fig. 8d.

Example 2. The only difference between Example 2 and the previous Example 1 is that delay time $\tau$ is estimated together with the initial function. Here is applied the similar gradient-based approach (and the same adjoint model) described in our previous work [8]. In order to obtain the gradient (scalar partial derivative) of the objective function w.r.t. delay time, the second output signal $\gamma(t)$ of the adjoint model presented in Fig. 7 has to be used. The reader interested in further details concerned with $\tau$ estimation is referred to our previous work [8]. The initial value of the delay time $\tau^0$ for the estimation procedure is 1.2 [s].

Once again, it can bee seen that the estimation process is convergent. The estimated value of $\tau$ reached the true value 1 [s] used in the plant － Fig. 9b. One can see that the value of the objective function is not strictly decreasing function of the iteration number and there is visible "bump". This is because we used the simplest gradient descent optimization procedure with constant $c$ parameter for which such bumps may appear. Of course they can be eliminated by reducing the parameter $c$ but at the cost of slowing down the process of estimation. Another possible approach is to apply more more sophisticated gradient-based optimization algorithms.

Example 3. In this example measurements are no longer quasi-countinuous. The sampling time $t_s=0.1~[\textrm{s}]$. The rest of of conditions are the same as in the Example 2, see Table 1. Results of this numerical example are presented in Fig. 10.

The estimated initial functions of the delay element differs from the true initial function used in the plant, see Fig. 10a. There are characteristic "jumps" caused by sampling. As previously, the delay time is estimated correctly － Fig. 10b. The difference between $y(t)$ and $d(t)$, i.e. the prediction error $e(t)$, presented in Fig. 10d, is significant and this discrepancy is caused by the the difference between the true and estimated initial function.

Nevertheless, the performance index, which takes into account only discrete-time measurements, reached a very small value. Furthermore, one can see that the prediction error $e(t)$ is significant only between sampling times. The prediction error taken at sampling times $e(n)$, see Fig. 10d, is close to zero. The conclusion coming from this example is that for a sampled-data system:

● the gradient of the objective function in the initial function space is calculated correctly,

● the output the model fits the discrete-time measurements and estimation procedure is convergent in the sense of the objective function value,

● the initial function of the delay, in general, is not convergent to the true initial function.

The last conclusion is more general. It is impossible to reconstruct perfectly a continuous function based only on discrete-time data, without further assumptions, like for example assumption about a frequency band limits in Nyquist-Shannon theorem. However, from the practical point of view, one can see that the initial function is estimated pretty well and it is close to the true function.

Example 4. In this example we show results of estimation of the delay time $\tau$ only. We also assume that there is no information about the true initial function and it is set to constant zero function in the model. We used the same stepwise initial function in the plant like in the previous examples. The initial value of the delay time $\tau$ for the estimation procedure is 1.2 [s] like in previous examples. Let us look at results of the gradient-based estimation process presented in Fig. 11.

One can see that the objective function is decreased. It means that the gradient of the objective function w.r.t the delay time is calculated correctly. However, the delay time is not estimated correctly. Even when $\tau$ reached the true value 1, see Fig. 11a, about 50-th iteration, the optimization procedure does not stop and continues to decrease $\tau$ until is reaches the lower constraint which is set to 0.

The conclusion coming from this example is that it is still worthwhile to estimate the initial functions of delays, even when we know that this estimate is not accurate (like in Example 3) or when we are not interested in information about the initial function at all. Simultaneous estimation of model's parameters and initial functions of delays improves estimation results of these parameters.

Example 5. In the next two examples we show results of initial functions estimation when there are more delays in the system. In Example 5 Model B with two delays, presented in Fig. 12, is used. The structure of the system is similar to the model A, except for the second delay acting in the upper branch.

The initial function for the second delay used in the plant is a sine wave presented in Fig. 13b － dashed line.

One can see from Fig. 13 that both initial function are estimated correctly like in Examples 1 and 2. Once again output signals of the model and the plant presented in Fig. 13d are nearly indistinguishable due to very small prediction error. Now, let us go to the next example where we will see a non-identifiable case.

Example 6. The structure of the model C used in this example is presented in Fig. 14. Like the model B, it also contains a second delay but placed in input channel (the input signal $u(t)$ is delayed by $\tau_2$).

Let us look at results of the gradient-based estimation process presented in Fig. 15. One can see that both initial function are estimated incorrectly, see Fig. 15a and 15b. Nevertheless, the output of the model is close to the output of the plant Fig. 15c. It means the solution is not unique and initial functions of these two delays are not identifiable. Besides the functions used in the plant, there are also other (at least two found in this example) optimal functions minimizing the performance index $J$. This effect can be explained when we analyze carefully how these two initial functions act in the system. Both functions influence additively the system ($\psi_1(t)$ with sign "$-$"and $\psi_2(t)$ with sign "$+$") through the same summing junction － see the structure of the model C from Fig. 14. It means that any change in function $\psi_1(t)$ can be compensated by change in function $\psi_2(t)$ and vice versa. The optimal solution have to preserve the difference between two initial functions used in the plant. Indeed, Fig. 15d, where these two differences (for the plant and the fitted model) are shown, confirms this observation.

---

6. Conclusions

In this work a gradient-based approach to estimation of initial functions for sampled non-linear systems with delays has been presented. The gradient of the appropriately defined quadratic objective function in the space of the initial function is obtained by using so-called structural sensitivity analysis.

Six numerical examples: for different sampling times, different number of delays and where delay time has been also estimated together with initial functions, have been presented. All these examples have shown that it is possible to efficiently calculate the gradient of the objective function in the space of initial functions.

Nevertheless, for some cases we have encountered the problem of non-identifiability － the objective function has been minimized, but the estimated initial function has differed from the reference initial function. It has been shown that discrete (non-continuous) measurements cause non-identifiability of continuous initial functions. The observed differences between outputs of the plant and the model comes from the nature of measurements － the output of the plant is measured only at (relatively rare) discrete moments, for which the prediction error is negligible but between them it stays significant. Two, or more initial functions can also be non-identifiable, even for (quasi-) continuous measurements.

Results obtained on this work suggest further investigation the of non-identifiability problem of initial functions. There are also some possibilities to decrease (not to eliminate) the prediction error between sampling times and they will be investigated in the future works.

---

Acknowledgments

This work was funded by the Polish National Science Centre under grant DEC-2013/11/B/ST7/01713. Part of calculations were performed on the Ziemowit computational cluster (http://www.ziemowit.hpc.polsl.pl) created in the POIG.02.01.00-00-166/08 project (BIO-FARMA) and expanded in the POIG.02.03.01-00-040/13 project (Syscancer). Ronald Hancock (Laval University, Laval, QC, Canada) is greatly acknowledged for valuable comments and English editing of the manuscript.

---