# On the numerical approximation of some inverse problems governed by nonlinear delay differential equation

Mourad Nachaoui, Abdeljalil Nachaoui, Tamaz Tadumadze

<figure><img src="image_1.png" /></figure>

## To cite this version:

Mourad Nachaoui, Abdeljalil Nachaoui, Tamaz Tadumadze. On the numerical approximation of some inverse problems governed by nonlinear delay differential equation. RAIRO - Operations Research, 2022, 56 (3), pp.1553-1569. 10.1051/ro/2022080. hal-03912544

HAL Id: hal-03912544

https://hal.science/hal-03912544

Submitted on 24 May 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the numerical approximation of some inverse problems governed by nonlinear delay differential equation

Mourad Nachaoui, Abdeljalil Nachaoui, Tamaz Tadumadze

<figure><img src="image_1.png" /></figure>

## To cite this version:

Mourad Nachaoui, Abdeljalil Nachaoui, Tamaz Tadumadze. On the numerical approximation of some inverse problems governed by nonlinear delay differential equation. RAIRO - Operations Research, 2022, 56 (3), pp.1553-1569. 10.1051/ro/2022080. hal-03912544

HAL Id: hal-03912544

https://hal.science/hal-03912544v1

Submitted on 24 May 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ON THE NUMERICAL APPROXIMATION OF SOME INVERSE PROBLEMS GOVERNED BY NONLINEAR DELAY DIFFERENTIAL EQUATION

MOURAD NACHAoui<fnref n="1" />*, ABDELJALIL NACHAoui<fnref n="2" /> AND TAMAZ TADUMADZE<fnref n="3" />

**Abstract.** The paper deals with the approximate solving of an inverse problem for the nonlinear delay differential equation, which consists of finding the initial moment and delay parameter based on some observed data. The inverse problem is considered as a nonlinear optimal control problem for which the necessary conditions of optimality are formulated and proved. The obtained optimal control problem is solved by a method based on an improved parallel evolutionary algorithm. The efficiency of the proposed approach is demonstrated through various numerical experiments.

**Mathematics Subject Classification.** 34K05 (34A12), 34A55, 65Y05, 68W50.

Received January 2, 2022. Accepted May 18, 2022.

## 1. INTRODUCTION

The recent development of applied mathematics is characterized by increasing attempts to use mathematical modeling tools in different engineering, biological and medicinal fields. This mathematical modeling, which generally based on ordinary and partial differential equations, attempt to better our understanding of more and more complicated phenomena. However, it is becoming clear that the simplest models cannot capture the rich variety of dynamics observed in natural phenomena. One of the more significant approaches is the inclusion of lag terms in differential equations. The delays or lags can represent gestation times, incubation periods, transport delays, or can simply lump complicated biological processes together, accounting only for the time required for these processes to occur [11]. Such models have the advantage of combining a simple, intuitive derivation with a wide variety of possible behavior regimes for a single system. Thereby, involving the delay in ordinary equations equations have been a topic of much interest in the mathematical research literature for more than 50 years. Contributions range from classical applications and theoretical and computational methodologies [9] to modern applications in biology [10, 15]. A particular interest has been established mathematical models describing the immune response during infectious diseases which are formulated as systems of nonlinear delay-differential equations (DDEs) characterized by multiple constant delays, moderate size, and stiffness [4, 5]. This paper deals with a topic that has become increasingly relevant in current research: inverse problems [1–3, 20, 24].

<footnote>*Keywords*. Inverse problem, delay differential equation, necessary optimality conditions, numerical approximation, parallel evolutionary algorithms.</footnote>

<footnote n="1">¹ (EMI), Université Sultan Moulay Slimane, Béni-Mellal, Morocco.</footnote>

<footnote n="2">² Laboratoire de Mathématiques Jean Leray UMR6629 CNRS / Université de Nantes, Nantes, France.</footnote>

<footnote n="3">³ Tbilisi State University, Tbilisi, Georgia.</footnote>

<footnote>*Corresponding author: nachaoui@gmail.com</footnote>
# ON THE NUMERICAL APPROXIMATION OF SOME INVERSE PROBLEMS GOVERNED BY NONLINEAR DELAY DIFFERENTIAL EQUATION

MOURAD NACHAoui<fnref n="1" />*, ABDELJALIL NACHAoui<fnref n="2" /> AND TAMAZ TADUMADZE<fnref n="3" />

**Abstract.** The paper deals with the approximate solving of an inverse problem for the nonlinear delay differential equation, which consists of finding the initial moment and delay parameter based on some observed data. The inverse problem is considered as a nonlinear optimal control problem for which the necessary conditions of optimality are formulated and proved. The obtained optimal control problem is solved by a method based on an improved parallel evolutionary algorithm. The efficiency of the proposed approach is demonstrated through various numerical experiments.

**Mathematics Subject Classification.** 34K05 (34A12), 34A55, 65Y05, 68W50.

Received January 2, 2022. Accepted May 18, 2022.

## 1. INTRODUCTION

The recent development of applied mathematics is characterized by increasing attempts to use mathematical modeling tools in different engineering, biological and medicinal fields. This mathematical modeling, which generally based on ordinary and partial differential equations, attempt to better our understanding of more and more complicated phenomena. However, it is becoming clear that the simplest models cannot capture the rich variety of dynamics observed in natural phenomena. One of the more significant approaches is the inclusion of lag terms in differential equations. The delays or lags can represent gestation times, incubation periods, transport delays, or can simply lump complicated biological processes together, accounting only for the time required for these processes to occur [11]. Such models have the advantage of combining a simple, intuitive derivation with a wide variety of possible behavior regimes for a single system. Thereby, involving the delay in ordinary equations equations have been a topic of much interest in the mathematical research literature for more than 50 years. Contributions range from classical applications and theoretical and computational methodologies [9] to modern applications in biology [10, 15]. A particular interest has been established mathematical models describing the immune response during infectious diseases which are formulated as systems of nonlinear delay-differential equations (DDEs) characterized by multiple constant delays, moderate size, and stiffness [4, 5]. This paper deals with a topic that has become increasingly relevant in current research: inverse problems [1–3, 20, 24].

<footnote>*Keywords*. Inverse problem, delay differential equation, necessary optimality conditions, numerical approximation, parallel evolutionary algorithms.</footnote>

<footnote n="1">¹ (EMI), Université Sultan Moulay Slimane, Béni-Mellal, Morocco.</footnote>

<footnote n="2">² Laboratoire de Mathématiques Jean Leray UMR6629 CNRS / Université de Nantes, Nantes, France.</footnote>

<footnote n="3">³ Tbilisi State University, Tbilisi, Georgia.</footnote>

<footnote>*Corresponding author: nachaoui@gmail.com</footnote>
In particular these involving nonlinear delay systems [29]. It appears in many applications where parameter estimation is ubiquitous. Inverse problems arise when our objective is to recover information about a system from observations of the system. The system is usually closed which refers to the property that the information we wish to recover cannot be directly observed.

In this work, we consider an inverse problem governed by a nonlinear delay differential system. It consists in finding the initial time and the delay that appear in the considered system of nonlinear DDE. We propose an optimal control approach to solve this inverse problem.

The inverse problems are highly ill-posed in Hadamard's sense [13]. There is no standard way to solve these kinds of problems. Therefore, the ill-posedness of inverse problems means that it needs to be somewhat modified to be solved. In this context, the optimal control formulation presents an excellent alternative to handle many inverse problems for differential equations. Indeed, the optimal controls bring a powerful theoretical framework and allow deriving optimality conditions as well as convergence rates. This explains the extended studies of inverse problems using this formulation [7, 18, 19, 26].

In this work, we study the existence of the optimal solution and we derive the necessary optimal conditions. In the numerical case, we develop a method based on the evolutionary algorithm (EA) and parallel technique running on a multiprocessor machine with distributed memory [22]. Indeed, as it's well known that EA is readily lending itself to parallel processing because each individual in the population is independent of the others and therefore can be processed concurrently. Thus, through many concrete problems, we present some numerical experiments illustrated by the performance study in terms of the approximation quality and time-consuming.

We highlight that the proposed theoretical and numerical parts are complementary. Indeed, to face the ill-posedness of the considered inverse problem, in the theoretical part we have shown the existence of an optimal solution and derived the necessary optimal conditions. Obviously, the optimality conditions are given only the existence of a local solution. By complementary, in the numerical part, we propose a numerical approach based on (EA) which is a robust algorithm that converges to a global solution [6, 22]. Consequently, the proposed approach in its entirety is an efficient tool to solve the considered ill-posed problem.

The rest of the paper is organized as follows. The Section 2 is devoted to the statement of the inverse problem. In Section 3, we derive the necessary optimality conditions. The proof of the main theorem is presented in Section 4. The numerical approach, the developed parallel algorithm, the numerical experiments with discussions, and the analysis of the performance of the parallel algorithm are presented in Section 5. The conclusion and remarks are presented in Section 6.

## 2. STATEMENT OF THE INVERSE PROBLEM

Let $\mathbb{R}^n$ be the $n$-dimensional vector space of points $x = (x^1, \dots, x^n)^T$ with $|x|^2 = \sum_{i=1}^n (x^i)^2$, where $T$ is the sign of transposition. Let $t_f > s_2 > s_1$ and $\theta_2 > \theta_f > 0$ be given numbers with $s_2 + \theta_2 < t_f$. Let the function $f(t, x, y) = (f^1(t, x, y), \dots, f^n(t, x, y))^T$, $(t, x, y) \in [s_1, t_f] \times \mathbb{R}^n \times \mathbb{R}^n$ be continuous and continuously differentiable with respect to $(x, y)$; there exists a number $M > 0$ such that for all $(t, x, y) \in [s_1, t_f] \times \mathbb{R}^n \times \mathbb{R}^n$ we have

$$|f(t, x, y)| + \left| \frac{\partial f}{\partial x}(t, x, y) \right| + \left| \frac{\partial f}{\partial y}(t, x, y) \right| \leq M.$$

Next, $\varphi(t) \in \mathbb{R}^n, t \in [\hat{\tau}, s_2]$, where $\hat{\tau} = s_1 - \theta_2$, is a given continuously differentiable function.

To each element $w = (t_0, \tau) \in W = [s_1, s_2] \times [\theta_f, \theta_2]$ we assign the delay differential equation

$$\dot{x}(t) = f(t, x(t), x(t - \tau)), t \in [t_0, t_f], \qquad (2.1)$$

with the initial condition

$$x(t) = \varphi(t), t \in [\hat{\tau}, t_0]. \qquad (2.2)$$
In particular these involving nonlinear delay systems [29]. It appears in many applications where parameter estimation is ubiquitous. Inverse problems arise when our objective is to recover information about a system from observations of the system. The system is usually closed which refers to the property that the information we wish to recover cannot be directly observed.

In this work, we consider an inverse problem governed by a nonlinear delay differential system. It consists in finding the initial time and the delay that appear in the considered system of nonlinear DDE. We propose an optimal control approach to solve this inverse problem.

The inverse problems are highly ill-posed in Hadamard's sense [13]. There is no standard way to solve these kinds of problems. Therefore, the ill-posedness of inverse problems means that it needs to be somewhat modified to be solved. In this context, the optimal control formulation presents an excellent alternative to handle many inverse problems for differential equations. Indeed, the optimal controls bring a powerful theoretical framework and allow deriving optimality conditions as well as convergence rates. This explains the extended studies of inverse problems using this formulation [7, 18, 19, 26].

In this work, we study the existence of the optimal solution and we derive the necessary optimal conditions. In the numerical case, we develop a method based on the evolutionary algorithm (EA) and parallel technique running on a multiprocessor machine with distributed memory [22]. Indeed, as it's well known that EA is readily lending itself to parallel processing because each individual in the population is independent of the others and therefore can be processed concurrently. Thus, through many concrete problems, we present some numerical experiments illustrated by the performance study in terms of the approximation quality and time-consuming.

We highlight that the proposed theoretical and numerical parts are complementary. Indeed, to face the ill-posedness of the considered inverse problem, in the theoretical part we have shown the existence of an optimal solution and derived the necessary optimal conditions. Obviously, the optimality conditions are given only the existence of a local solution. By complementary, in the numerical part, we propose a numerical approach based on (EA) which is a robust algorithm that converges to a global solution [6, 22]. Consequently, the proposed approach in its entirety is an efficient tool to solve the considered ill-posed problem.

The rest of the paper is organized as follows. The Section 2 is devoted to the statement of the inverse problem. In Section 3, we derive the necessary optimality conditions. The proof of the main theorem is presented in Section 4. The numerical approach, the developed parallel algorithm, the numerical experiments with discussions, and the analysis of the performance of the parallel algorithm are presented in Section 5. The conclusion and remarks are presented in Section 6.

## 2. STATEMENT OF THE INVERSE PROBLEM

Let $\mathbb{R}^n$ be the $n$-dimensional vector space of points $x = (x^1, \dots, x^n)^T$ with $|x|^2 = \sum_{i=1}^n (x^i)^2$, where $T$ is the sign of transposition. Let $t_f > s_2 > s_1$ and $\theta_2 > \theta_f > 0$ be given numbers with $s_2 + \theta_2 < t_f$. Let the function $f(t, x, y) = (f^1(t, x, y), \dots, f^n(t, x, y))^T$, $(t, x, y) \in [s_1, t_f] \times \mathbb{R}^n \times \mathbb{R}^n$ be continuous and continuously differentiable with respect to $(x, y)$; there exists a number $M > 0$ such that for all $(t, x, y) \in [s_1, t_f] \times \mathbb{R}^n \times \mathbb{R}^n$ we have

$$|f(t, x, y)| + \left| \frac{\partial f}{\partial x}(t, x, y) \right| + \left| \frac{\partial f}{\partial y}(t, x, y) \right| \leq M.$$

Next, $\varphi(t) \in \mathbb{R}^n, t \in [\hat{\tau}, s_2]$, where $\hat{\tau} = s_1 - \theta_2$, is a given continuously differentiable function.

To each element $w = (t_0, \tau) \in W = [s_1, s_2] \times [\theta_f, \theta_2]$ we assign the delay differential equation

$$\dot{x}(t) = f(t, x(t), x(t - \tau)), t \in [t_0, t_f], \quad (2.1)$$

with the initial condition

$$x(t) = \varphi(t), t \in [\hat{\tau}, t_0]. \quad (2.2)$$
**Definition 2.1.** Let $w = (t_0, \tau) \in W$, a function $x(t) = x(t; w) \in \mathbb{R}^n$, $t \in [\hat{\tau}, t_f]$ is called a solution of differential equation (2.1) with the initial condition (2.2) or a solution corresponding to the element $w$ if $x(t)$ satisfies the initial condition (2.2) is continuously differentiable on the interval $[t_0, t_f]$ and satisfies equation (2.1) everywhere.

By the step method and Gronwall inequality can be proved that for every element $w \in W$ there exists a unique solution $x(t; w)$ of (2.1) and (2.2) defined on the interval $[\hat{\tau}, t_f]$ and it is continuous with respect to $w$.

Let us introduce the set

$$Z = \{z \in \mathbb{R}^n : \exists w \in W, x(t_f; w) = z\}.$$

## 2.1. Inverse problem

Let $z \in Z$ be a given vector. Find element $w \in W$ such that the following condition holds $x(t_f; w) = z$.

The vector $z$, as rule, by distinct error is beforehand given. Thus, instead of the vector $z$ we have $\hat{z}$ (so called an observed vector) which is an approximation to the $z$ and, in general, $\hat{z} \notin Z$. Therefore it is natural to change posed inverse problem by the following approximate problem.

## 2.2. The approximate inverse problem

Find an element $w \in W$ such that the deviation

$$J(w) = \frac{1}{2} |x(t_f; w) - \hat{z}|^2,$$

takes the minimal value. Where $x(t)$, $t \in [t_0, t_f]$ is a solution of (2.1) and (2.2).

It is clear that the approximate inverse problem is equivalent to the following optimal control problem:

$$\left\{
\begin{array}{ll}
\min_{w \in W} J(w) := \frac{1}{2} |x(t_f; w) - \hat{z}|^2 & \\
\text{where } x(t) \text{ is solution of} & (2.3) \\
\dot{x}(t) = f(t, x(t), x(t-\tau)), t \in [t_0, t_f], & \\
x(t) = \varphi(t), t \in [\hat{\tau}, t_0]. &
\end{array}
\right.$$

The problem (2.3) is called optimal control problem corresponding to the approximate inverse problem. The element $w_0 = (t_{00}, \tau_0) \in W$ is called optimal if

$$J(w_0) = \min_{w \in W} J(w).$$

For the equivalence, if $w^* \in W$ is a solution of the approximate inverse problem, then there exists $x(t)$, $t \in [t_0, t_f]$ solution of (2.1) and (2.2) such that $J(w^*) \le J(w) \forall w \in W$. This means that $w$ is a solution of (2.3).

Conversely, if $\hat{w}$ is a solution of the optimal control problem (2.3), then we have

$$\begin{aligned}
J(\hat{w}) &\le J(w) \quad \forall w \in W, \\
&\le J(w^*)
\end{aligned}$$

which means that $\hat{w}$ is a solution of the approximate inverse problem.

From continuity of the function $J(w)$, $w \in W$ it follows the existence of an optimal element $w_0$.

## 3. NECESSARY OPTIMALITY CONDITIONS

In this section we give some results allow to give the necessary optimality conditions of the control problem (2.3).

**Theorem 3.1.** Let $w_0 = (t_{00}, \tau_0) \in (s_1, s_2) \times (\theta_f, \theta_2)$ be an optimal element. Then the following conditions hold:
**Definition 2.1.** Let $w = (t_0, \tau) \in W$, a function $x(t) = x(t; w) \in \mathbb{R}^n$, $t \in [\hat{\tau}, t_f]$ is called a solution of differential equation (2.1) with the initial condition (2.2) or a solution corresponding to the element $w$ if $x(t)$ satisfies the initial condition (2.2) is continuously differentiable on the interval $[t_0, t_f]$ and satisfies equation (2.1) everywhere.

By the step method and Gronwall inequality can be proved that for every element $w \in W$ there exists a unique solution $x(t; w)$ of (2.1) and (2.2) defined on the interval $[\hat{\tau}, t_f]$ and it is continuous with respect to $w$.

Let us introduce the set

$$Z = \{z \in \mathbb{R}^n : \exists w \in W, x(t_f; w) = z\}.$$

## 2.1. Inverse problem

Let $z \in Z$ be a given vector. Find element $w \in W$ such that the following condition holds $x(t_f; w) = z$.

The vector $z$, as rule, by distinct error is beforehand given. Thus, instead of the vector $z$ we have $\hat{z}$ (so called an observed vector) which is an approximation to the $z$ and, in general, $\hat{z} \notin Z$. Therefore it is natural to change posed inverse problem by the following approximate problem.

## 2.2. The approximate inverse problem

Find an element $w \in W$ such that the deviation

$$J(w) = \frac{1}{2} |x(t_f; w) - \hat{z}|^2,$$

takes the minimal value. Where $x(t)$, $t \in [t_0, t_f]$ is a solution of (2.1) and (2.2).

It is clear that the approximate inverse problem is equivalent to the following optimal control problem:

$$\left\{
\begin{array}{ll}
\min_{w \in W} J(w) := \frac{1}{2} |x(t_f; w) - \hat{z}|^2 & \\
\text{where } x(t) \text{ is solution of} & (2.3) \\
\dot{x}(t) = f(t, x(t), x(t-\tau)), t \in [t_0, t_f], & \\
x(t) = \varphi(t), t \in [\hat{\tau}, t_0]. &
\end{array}
\right.$$

The problem (2.3) is called optimal control problem corresponding to the approximate inverse problem. The element $w_0 = (t_{00}, \tau_0) \in W$ is called optimal if

$$J(w_0) = \min_{w \in W} J(w).$$

For the equivalence, if $w^* \in W$ is a solution of the approximate inverse problem, then there exists $x(t)$, $t \in [t_0, t_f]$ solution of (2.1) and (2.2) such that $J(w^*) \le J(w) \forall w \in W$. This means that $w$ is a solution of (2.3).

Conversely, if $\hat{w}$ is a solution of the optimal control problem (2.3), then we have

$$\begin{array}{ll}
J(\hat{w}) \le J(w) & \forall w \in W, \\
\le J(w^*) &
\end{array}$$

which means that $\hat{w}$ is a solution of the approximate inverse problem.

From continuity of the function $J(w)$, $w \in W$ it follows the existence of an optimal element $w_0$.

## 3. NECESSARY OPTIMALITY CONDITIONS

In this section we give some results allow to give the necessary optimality conditions of the control problem (2.3).

**Theorem 3.1.** Let $w_0 = (t_{00}, \tau_0) \in (s_1, s_2) \times (\theta_f, \theta_2)$ be an optimal element. Then the following conditions hold:
(1) the condition for the optimal initial moment $t_{00}$

$$ \psi(t_{00})[\dot{\varphi}(t_{00}) - f(t_{00}, \varphi(t_{00}), \varphi(t_{00} - \tau_0))] = 0; $$

(2) the condition for the optimal delay parameter $\tau_0$

$$ \int_{t_{00}}^{t_f} \psi(t) f_y(t, x_0(t), x_0(t - \tau_0)) \dot{x}_0(t - \tau_0) dt = 0, $$

where $x_0(t) = x(t; w_0)$; Here $\psi(t)$ is a solution of the system

$$ \begin{aligned} \dot{\psi}(t) = & -\psi(t) \frac{\partial f}{\partial x}(t, x_0(t), x_0(t - \tau_0)) \\ & - \psi(t + \tau_0) f_y(t + \tau_0, x_0(t + \tau_0), x_0(t)), t \in [t_{00}, t_f] \end{aligned} \quad (3.1) $$

with the initial condition

$$ \psi(t_f) = -(x_0(t_f) - \hat{z})^T, \psi(t) = 0, t > t_f. \quad (3.2) $$

**Remark 3.2.** It is clear that the condition (2) is equivalent to the following relation

$$ \int_{t_{00}}^{t_{00}+\tau_0} \psi(t) f_y(t, x_0(t), \varphi(t-\tau_0)) \dot{\varphi}(t-\tau_0) dt \\ + \int_{t_{00}+\tau_0}^{t_1} \psi(t) \frac{\partial f}{\partial y}(t, x_0(t), x_0(t-\tau_0)) f(t-\tau_0, x_0(t-\tau_0), x_0(t-2\tau_0)) dt = 0. $$

**Theorem 3.3.** Let $w_0 = (s_1, \theta_f) \in W$ be an optimal element. Then the following conditions hold:

(3) the condition for the optimal initial moment $s_1$

$$ \psi(s_1)[\dot{\varphi}(s_1) - f(s_1, \varphi(s_1), \varphi(s_1 - \theta_f))] \le 0, $$

(4) the condition for the delay parameter $\tau_0$

$$ \int_{s_1}^{t_f} \psi(t) \frac{\partial f}{\partial y}(t, x_0(t), x_0(t - \theta_f)) \dot{x}_0(t - \tau_0) dt \ge 0. $$

Here $\psi(t)$ satisfies the equation (3.1) and the condition (3.2), where $t_{00} = s_1$ and $\tau_0 = \theta_f$.

**Theorem 3.4.** Let $w_0 = (s_2, \theta_2) \in W$ be an optimal element. Then the following conditions hold:

(5) the condition for the optimal initial moment $s_2$

$$ \psi(s_2)[\dot{\varphi}(s_2) - f(s_2, \varphi(s_2), \varphi(s_2 - \theta_2))] \ge 0, $$

(6) the condition for the optimal delay parameter $\tau_0$

$$ \int_{s_2}^{t_f} \psi(t) \frac{\partial f}{\partial y}(t, x_0(t), x_0(t - \theta_2)) \dot{x}_0(t - \tau_0) dt \le 0. $$

Here $\psi(t)$ satisfies the equation (3.1) and the condition (3.2), where $t_{00} = s_2$ and $\tau_0 = \theta_2$.
(1) the condition for the optimal initial moment $t_{00}$

$$ \psi(t_{00})[\dot{\varphi}(t_{00}) - f(t_{00}, \varphi(t_{00}), \varphi(t_{00} - \tau_0))] = 0; $$

(2) the condition for the optimal delay parameter $\tau_0$

$$ \int_{t_{00}}^{t_f} \psi(t) f_y(t, x_0(t), x_0(t - \tau_0)) \dot{x}_0(t - \tau_0) dt = 0, $$

where $x_0(t) = x(t; w_0)$; Here $\psi(t)$ is a solution of the system

$$ \begin{aligned} \dot{\psi}(t) = & -\psi(t) \frac{\partial f}{\partial x}(t, x_0(t), x_0(t - \tau_0)) \\ & - \psi(t + \tau_0) f_y(t + \tau_0, x_0(t + \tau_0), x_0(t)), t \in [t_{00}, t_f] \end{aligned} \quad (3.1) $$

with the initial condition

$$ \psi(t_f) = -(x_0(t_f) - \hat{z})^T, \psi(t) = 0, t > t_f. \quad (3.2) $$

**Remark 3.2.** It is clear that the condition (2) is equivalent to the following relation

$$ \int_{t_{00}}^{t_{00}+\tau_0} \psi(t) f_y(t, x_0(t), \varphi(t-\tau_0)) \dot{\varphi}(t-\tau_0) dt \\ + \int_{t_{00}+\tau_0}^{t_1} \psi(t) \frac{\partial f}{\partial y}(t, x_0(t), x_0(t-\tau_0)) f(t-\tau_0, x_0(t-\tau_0), x_0(t-2\tau_0)) dt = 0. $$

**Theorem 3.3.** Let $w_0 = (s_1, \theta_f) \in W$ be an optimal element. Then the following conditions hold:

(3) the condition for the optimal initial moment $s_1$

$$ \psi(s_1)[\dot{\varphi}(s_1) - f(s_1, \varphi(s_1), \varphi(s_1 - \theta_f))] \le 0, $$

(4) the condition for the delay parameter $\tau_0$

$$ \int_{s_1}^{t_f} \psi(t) \frac{\partial f}{\partial y}(t, x_0(t), x_0(t - \theta_f)) \dot{x}_0(t - \tau_0) dt \ge 0. $$

Here $\psi(t)$ satisfies the equation (3.1) and the condition (3.2), where $t_{00} = s_1$ and $\tau_0 = \theta_f$.

**Theorem 3.4.** Let $w_0 = (s_2, \theta_2) \in W$ be an optimal element. Then the following conditions hold:

(5) the condition for the optimal initial moment $s_2$

$$ \psi(s_2)[\dot{\varphi}(s_2) - f(s_2, \varphi(s_2), \varphi(s_2 - \theta_2))] \ge 0, $$

(6) the condition for the optimal delay parameter $\tau_0$

$$ \int_{s_2}^{t_f} \psi(t) \frac{\partial f}{\partial y}(t, x_0(t), x_0(t - \theta_2)) \dot{x}_0(t - \tau_0) dt \le 0. $$

Here $\psi(t)$ satisfies the equation (3.1) and the condition (3.2), where $t_{00} = s_2$ and $\tau_0 = \theta_2$.
**Remark 3.5.** The conditions (4) and (5) are equivalent to the following relations

$$ \int_{s_1}^{s_1+\theta_f} \psi(t) \frac{\partial f}{\partial y}(t, x_0(t), \varphi(t-\theta_f)) \dot{\varphi}(t-\theta_f) \, dt \\ + \int_{s_1+\theta_f}^{t_1} \psi(t) \frac{\partial f}{\partial y}(t, x_0(t), x_0(t-\theta_f)) f(t-\theta_f, x_0(t-\theta_f), x_0(t-2\theta_f)) \, dt \geq 0. $$

and

$$ \int_{s_2}^{s_2+\theta_2} \psi(t) \frac{\partial f}{\partial y}(t, x_0(t), \varphi(t-\theta_2)) \dot{\varphi}(t-\theta_2) \, dt \\ + \int_{s_2+\theta_2}^{t_1} \psi(t) \frac{\partial f}{\partial y}(t, x_0(t), x_0(t-\theta_2)) f(t-\theta_2, x_0(t-\theta_2), x_0(t-2\theta_2)) \, dt \leq 0 $$

respectively.

**Theorem 3.6.** Let $w_0 = (s_1, \theta_2) \in W$ be an optimal element. Then the conditions (3) and (6) hold.

**Theorem 3.7.** Let $w_0 = (s_2, \theta_f) \in W$ be an optimal element. Then the conditions (5) and (4) hold.

## 4. PROOF OF MAIN THEOREMS

In this section we give the proof of the theorems presented in the preview section. We will give the proof of the Theorems 3.1 and 3.3, for other theorems they are proved by the analogous scheme.

On the convex set $D = \mathbb{R}_+ \times W$, where $\mathbb{R}_+ = [0, \infty)$, let us define the mapping

$$ P: D \to \mathbb{R} \qquad (4.1) $$

by the formula

$$ P(v) = \frac{1}{2}|x(t_f; w) - \hat{z}|^2 + \xi, v = (\xi, w) \in D. $$

It is clear that

$$ P(v_0) \leq P(v), \forall v \in D, $$

where $v_0 = (0, w_0)$.

Thus the point $v_0 = (0, w_0) \in D$ is a critical point (see [12, 17]), since $P(v_0) \in \partial P(D)$. Where $\partial P(D)$ is the boundary of the $P(D)$. Moreover, the mapping (4.1) is continuous.

There exists a small $\varepsilon_0 > 0$ such that for an arbitrary $\varepsilon \in (0, \varepsilon_0)$ and $\delta v = (\delta \xi, \delta w) \in V_{v_0} \subset D - v_0$, where $\delta w = (\delta t_0, \delta \tau)$, $V_{v_0} = [0, \alpha) \times (t_{00} - \alpha, t_{00} + \alpha) \times (\tau_0 - \alpha, \tau_0 + \alpha)$ we get $v_0 + \varepsilon \delta v \in D$.

Now we calculate a differential of the mapping (4.1) at the point $v_0$. We have,

$$ \begin{aligned} P(v_0 + \varepsilon \delta v) - P(v_0) &= \frac{1}{2}|x(t_f; w_0 + \varepsilon \delta w) - \hat{z}|^2 + \varepsilon \delta \xi - \frac{1}{2}|x(t_f; w_0) - \hat{z}|^2 \\ &= \frac{1}{2} \int_0^1 \frac{d}{ds}|x(t_f; w_0) + s(x(t_f; w_0 + \varepsilon \delta w) - x(t_f; w_0)) - \hat{z}|^2 ds + \varepsilon \delta \xi \\ &= \frac{1}{2} \sum_{i=1}^n \int_0^1 \frac{d}{ds}[x^i(t_f; w_0) + s(x^i(t_f; w_0 + \varepsilon \delta w) - x^i(t_f; w_0)) - \hat{z}^i]^2 ds + \varepsilon \delta \xi \\ &= \sum_{i=1}^n \int_0^1 [x^i(t_f; w_0) + s(x^i(t_f; w_0 + \varepsilon \delta w) - x^i(t_f; w_0)) - \hat{z}^i] [x^i(t_f; w_0 + \varepsilon \delta w) - \hat{z}^i] ds \end{aligned} $$
**Remark 3.5.** The conditions (4) and (5) are equivalent to the following relations

$$ \int_{s_1}^{s_1+\theta_f} \psi(t) \frac{\partial f}{\partial y}(t, x_0(t), \varphi(t-\theta_f)) \dot{\varphi}(t-\theta_f) \, dt \\ + \int_{s_1+\theta_f}^{t_1} \psi(t) \frac{\partial f}{\partial y}(t, x_0(t), x_0(t-\theta_f)) f(t-\theta_f, x_0(t-\theta_f), x_0(t-2\theta_f)) \, dt \geq 0. $$

and

$$ \int_{s_2}^{s_2+\theta_2} \psi(t) \frac{\partial f}{\partial y}(t, x_0(t), \varphi(t-\theta_2)) \dot{\varphi}(t-\theta_2) \, dt \\ + \int_{s_2+\theta_2}^{t_1} \psi(t) \frac{\partial f}{\partial y}(t, x_0(t), x_0(t-\theta_2)) f(t-\theta_2, x_0(t-\theta_2), x_0(t-2\theta_2)) \, dt \leq 0 $$

respectively.

**Theorem 3.6.** Let $w_0 = (s_1, \theta_2) \in W$ be an optimal element. Then the conditions (3) and (6) hold.

**Theorem 3.7.** Let $w_0 = (s_2, \theta_f) \in W$ be an optimal element. Then the conditions (5) and (4) hold.

## 4. PROOF OF MAIN THEOREMS

In this section we give the proof of the theorems presented in the preview section. We will give the proof of the Theorems 3.1 and 3.3, for other theorems they are proved by the analogous scheme.

On the convex set $D = \mathbb{R}_+ \times W$, where $\mathbb{R}_+ = [0, \infty)$, let us define the mapping

$$ P: D \to \mathbb{R} \qquad (4.1) $$

by the formula

$$ P(v) = \frac{1}{2} |x(t_f; w) - \hat{z}|^2 + \xi, \quad v = (\xi, w) \in D. $$

It is clear that

$$ P(v_0) \leq P(v), \forall v \in D, $$

where $v_0 = (0, w_0)$.

Thus the point $v_0 = (0, w_0) \in D$ is a critical point (see [12, 17]), since $P(v_0) \in \partial P(D)$. Where $\partial P(D)$ is the boundary of the $P(D)$. Moreover, the mapping (4.1) is continuous.

There exists a small $\varepsilon_0 > 0$ such that for an arbitrary $\varepsilon \in (0, \varepsilon_0)$ and $\delta v = (\delta \xi, \delta w) \in V_{v_0} \subset D - v_0$, where $\delta w = (\delta t_0, \delta \tau)$, $V_{v_0} = [0, \alpha) \times (t_{00} - \alpha, t_{00} + \alpha) \times (\tau_0 - \alpha, \tau_0 + \alpha)$ we get $v_0 + \varepsilon \delta v \in D$.

Now we calculate a differential of the mapping (4.1) at the point $v_0$. We have,

$$ \begin{aligned} P(v_0 + \varepsilon \delta v) - P(v_0) &= \frac{1}{2} |x(t_f; w_0 + \varepsilon \delta w) - \hat{z}|^2 + \varepsilon \delta \xi - \frac{1}{2} |x(t_f; w_0) - \hat{z}|^2 \\ &= \frac{1}{2} \int_0^1 \frac{d}{ds} |x(t_f; w_0) + s(x(t_f; w_0 + \varepsilon \delta w) - x(t_f; w_0)) - \hat{z}|^2 ds + \varepsilon \delta \xi \\ &= \frac{1}{2} \sum_{i=1}^n \int_0^1 \frac{d}{ds} [x^i(t_f; w_0) + s(x^i(t_f; w_0 + \varepsilon \delta w) - x^i(t_f; w_0)) - \hat{z}^i]^2 ds + \varepsilon \delta \xi \\ &= \sum_{i=1}^n \int_0^1 [x^i(t_f; w_0) + s(x^i(t_f; w_0 + \varepsilon \delta w) - x^i(t_f; w_0)) - \hat{z}^i] [x^i(t_f; w_0 + \varepsilon \delta w) - \hat{z}^i] ds \end{aligned} $$
$$ 
\begin{aligned}
& -x^i(t_f; w_0)] ds = \sum_{i=1}^{n} [x^i(t_f; w_0) - \hat{z}^i] [x^i(t_f; w_0 + \varepsilon\delta w) - x^i(t_f; w_0)] \\
& + \frac{1}{2} [x^i(t_f; w_0 + \varepsilon\delta w) - x^i(t_f; w_0)]^2 = (x(t_f; w_0) - \hat{z})^T \Delta x(t_f; \varepsilon\delta w) \\
& + \frac{1}{2} |\Delta x(t_f; \varepsilon\delta w)|^2 + \varepsilon\delta\xi,
\end{aligned} 
\quad (4.2) $$

where

$$ \Delta x(t_f; \varepsilon\delta w) = x(t_f; w_0 + \varepsilon\delta w) - x(t_f; w_0). $$

Taking into account the relation

$$ \Delta x(t_f; \varepsilon\delta w) = \varepsilon\delta x(t_f; \delta w) + o(\varepsilon\delta w) $$

and the variation formula

$$ \delta x(t_f; \delta w) = Y(t_{00}; t_f)[\dot{\varphi}(t_{00}) - f(t_{00}, \varphi(t_{00}), \varphi(t_{00} - \tau_0))]\delta t_0 \\ - \left\{ \int_{t_{00}}^{t_f} Y(t; t_f) \frac{\partial f}{\partial y}(t, x_0(t), x_0(t - \tau_0)) \dot{x}_0(t - \tau_0) ds \right\} \delta\tau $$

(see [28]), from (4.2) by elementary transformation we obtain

$$ 
\begin{aligned}
P(v_0 + \varepsilon\delta v) - P(v_0) &= \varepsilon \left[ (x(t_f; w_0) - \hat{z})^T Y(t_{00}; t_f) [\dot{\varphi}(t_{00}) - f(t_{00}, \varphi(t_{00}), \varphi(t_{00} - \tau_0))] \delta t_0 \right. \\
&\quad \left. - \left\{ \int_{t_{00}}^{t_f} (x(t_1; w_0) - \hat{z})^T Y(t; t_f) \frac{\partial f}{\partial y}(t, x_0(t), x_0(t - \tau_0)) \dot{x}_0(t - \tau_0) dt \right\} \delta\tau + \delta\xi \right] + o(\varepsilon\delta w) \\
&= \varepsilon dP_{v_0}(\delta v) + o(\varepsilon\delta w).
\end{aligned} 
\quad (4.3) $$

Here, $dP_{v_0}(\delta v)$ denotes the differential of mapping (4.1); $Y(s; t_f)$ is the $n \times n$-matrix function satisfying the linear differential equation with advanced argument

$$ \frac{d}{ds} Y(t; t_f) = -Y(t; t_f) \frac{\partial f}{\partial x}(t, x_0(t), x_0(t - \tau_0)) - Y(t + \tau_0; t_f) \frac{\partial f}{\partial y}(t + \tau_0, x_0(t + \tau_0), x_0(t)), t \in [t_{00}, t_f] $$

and the conditions

$$ Y(t; t_f) = \begin{cases} I & \text{for } t = t_f \\ \Theta & \text{for } t > t_f, \end{cases} $$

$I$ is the identity matrix and $\Theta$ is the zero matrix;

$$ \lim_{\varepsilon \to 0} \frac{o(\varepsilon\delta w)}{\varepsilon} = 0 \text{ uniformly for } \delta v \in V_{v_0}. $$

From the necessary condition of criticality (see [12,17]) it follows that: there exists a number $\pi \neq 0$ such that

$$ \pi dP_{v_0}(\delta v) \leq 0, \forall \delta v \in \text{cone}(V_{v_0} - v_0). \quad (4.4) $$

Introduce the function

$$ \psi(t) = \pi(x(t_f; w_0) - \hat{z})^T Y(t; t_f) \quad (4.5) $$

as is easily seen, it satisfies the equation (3.1) and the condition

$$ \psi(t_f) = \pi(x(t_f; w_0) - \hat{z})^T. \quad (4.6) $$
$$ 
\begin{aligned}
& -x^i(t_f; w_0)] ds = \sum_{i=1}^{n} [x^i(t_f; w_0) - \hat{z}^i] [x^i(t_f; w_0 + \varepsilon\delta w) - x^i(t_f; w_0)] \\
& + \frac{1}{2} [x^i(t_f; w_0 + \varepsilon\delta w) - x^i(t_f; w_0)]^2 = (x(t_f; w_0) - \hat{z})^T \Delta x(t_f; \varepsilon\delta w) \\
& + \frac{1}{2} |\Delta x(t_f; \varepsilon\delta w)|^2 + \varepsilon\delta\xi,
\end{aligned} 
\quad (4.2) $$

where

$$ \Delta x(t_f; \varepsilon\delta w) = x(t_f; w_0 + \varepsilon\delta w) - x(t_f; w_0). $$

Taking into account the relation

$$ \Delta x(t_f; \varepsilon\delta w) = \varepsilon\delta x(t_f; \delta w) + o(\varepsilon\delta w) $$

and the variation formula

$$ \delta x(t_f; \delta w) = Y(t_{00}; t_f)[\dot{\varphi}(t_{00}) - f(t_{00}, \varphi(t_{00}), \varphi(t_{00} - \tau_0))]\delta t_0 \\ - \left\{ \int_{t_{00}}^{t_f} Y(t; t_f) \frac{\partial f}{\partial y}(t, x_0(t), x_0(t - \tau_0)) \dot{x}_0(t - \tau_0) ds \right\} \delta\tau $$

(see [28]), from (4.2) by elementary transformation we obtain

$$ 
\begin{aligned}
P(v_0 + \varepsilon\delta v) - P(v_0) &= \varepsilon \left[ (x(t_f; w_0) - \hat{z})^T Y(t_{00}; t_f) [\dot{\varphi}(t_{00}) - f(t_{00}, \varphi(t_{00}), \varphi(t_{00} - \tau_0))] \delta t_0 \right. \\
&\quad \left. - \left\{ \int_{t_{00}}^{t_f} (x(t_1; w_0) - \hat{z})^T Y(t; t_f) \frac{\partial f}{\partial y}(t, x_0(t), x_0(t - \tau_0)) \dot{x}_0(t - \tau_0) dt \right\} \delta\tau + \delta\xi \right] + o(\varepsilon\delta w) \\
&= \varepsilon dP_{v_0}(\delta v) + o(\varepsilon\delta w).
\end{aligned} 
\quad (4.3) $$

Here, $dP_{v_0}(\delta v)$ denotes the differential of mapping (4.1); $Y(s; t_f)$ is the $n \times n$-matrix function satisfying the linear differential equation with advanced argument

$$ \frac{d}{ds} Y(t; t_f) = -Y(t; t_f) \frac{\partial f}{\partial x}(t, x_0(t), x_0(t - \tau_0)) - Y(t + \tau_0; t_f) \frac{\partial f}{\partial y}(t + \tau_0, x_0(t + \tau_0), x_0(t)), t \in [t_{00}, t_f] $$

and the conditions

$$ Y(t; t_f) = \begin{cases} I & \text{for } t = t_f \\ \Theta & \text{for } t > t_f, \end{cases} $$

$I$ is the identity matrix and $\Theta$ is the zero matrix;

$$ \lim_{\varepsilon \to 0} \frac{o(\varepsilon\delta w)}{\varepsilon} = 0 \text{ uniformly for } \delta v \in V_{v_0}. $$

From the necessary condition of criticality (see [12,17]) it follows that: there exists a number $\pi \neq 0$ such that

$$ \pi dP_{v_0}(\delta v) \leq 0, \forall \delta v \in \text{cone}(V_{v_0} - v_0). \quad (4.4) $$

Introduce the function

$$ \psi(t) = \pi(x(t_f; w_0) - \hat{z})^T Y(t; t_f) \quad (4.5) $$

as is easily seen, it satisfies the equation (3.1) and the condition

$$ \psi(t_f) = \pi(x(t_f; w_0) - \hat{z})^T. \quad (4.6) $$
Taking into account (4.3), (4.5) and (4.6) from inequality (4.4) we obtain

$$ \psi(t_{00})[\dot{\varphi}(t_{00}) - f(t_{00}, \varphi(t_{00}), \varphi(t_{00} - \tau_0))]\delta t_0 - \left\{ \int_{t_{00}}^{t_f} \psi(t) \frac{\partial f}{\partial y}(t, x_0(t), x_0(t - \tau_0)) \dot{x}_0(t - \tau_0) dt \right\} \delta \tau + \pi \delta \xi \le 0, \\ \delta v \in \text{cone}(V_{v_0} - v_0). \tag{4.7} $$

The condition $\delta v \in \text{cone}(V_{v_0} - v_0)$ is equivalent to conditions:

$$ \delta\xi \in \mathbb{R}_{+}, \delta t_0 \in \mathbb{R}, \delta\tau \in \mathbb{R}. $$

Let $\delta t_0 = \delta \tau = 0$ in (4.7), then we obtain

$$ \pi \delta \xi \le 0, \forall \delta \xi \in \mathbb{R}_{+}. $$

This implies $\pi < 0$, therefore we can take $\pi = -1$, i.e., $\psi$ satisfies the condition (3.2).

Let $\delta\tau = \delta\xi = 0$ then, taking into account that $\delta t_0 \in \mathbb{R}$ from (4.7) we obtain the condition (1).

Let $\delta t_0 = \delta \xi = 0$ then, taking into account that $\delta \tau \in \mathbb{R}$ from (4.7) we obtain the condition (2).

Finally, we note that in order to prove the Theorem 3.3 it suffices to replace the set $V_{v_0}$ by the set

$$ [0, \alpha) \times [s_1, s_1 + \alpha) \times [\theta_1, \theta_0 + \alpha) $$

see the proof of the Theorem 3.1.

## 5. NUMERICAL APPROACH

In this section we shall present the numerical approach to solve the problem (2.3). First of all, we replace the constrained problem by its appropriate discrete optimal control problem. We subdivide the interval $[t_0, t_f]$ into $N$ sub-intervals with knots $t_0 < t_1 < \dots < t_f$ and $t_k = k\Delta t_k$, where $\Delta t_k$ is the mesh size of $k$th sub-interval. Where $k \in \{0, 1, \dots, N\}$.

The discrete problem associates to the continuous problem (2.3) is given by

$$ \left\{ \begin{array}{l} \min J(W, x_N) \\ x_{k+1} = \Psi(\Delta t_k, f(t_k, x_k, x_k(t_k - \tau))), \end{array} \right. \tag{5.1} $$

where $\Psi$ is a method of approximation of the delay differential equations (2.1) and (2.2). In this paper, the state problem will be solved by the solver dde23 using the popular approach based on the Runge-Kutta triple BS(2,3) [27].

### 5.1. Numerical algorithm

Evolutionary Algorithms are based on a model of natural, biological evolution, which was primarily developed by Holland [14]. It is essentially a searching method based on the Darwinian principles of biological evolution that explains the adaptive change of species by the principle of natural selection, which favors those species for survival and further evolution that are best adapted to their environmental conditions. It has been successfully applied to various optimization problems [6–8, 19].

In the evolutionary framework, a new generation of individuals is produced using the simulated genetic operations selection, crossover, and mutation. The evolution is based only on fitness. The fitness of an individual is measured only indirectly by its growth rate in comparison to others. Selection is the ability of those individuals that have outlasted the struggle for existence to bring their genetic information to the next generation. The crossover is the mechanism of reproduction based on the combination of the selected individuals. While the mutation is characterized by modifying one or more genes of an individual with some probability. The object of the mutation operator is to restore the lost or unexplored genetic material into the population to stop the premature convergence of the evolutionary algorithm to sub-optimal solutions.
Taking into account (4.3), (4.5) and (4.6) from inequality (4.4) we obtain

$$ \begin{aligned} & \psi(t_{00})[\dot{\varphi}(t_{00}) - f(t_{00}, \varphi(t_{00}), \varphi(t_{00} - \tau_0))]\delta t_0 - \left\{ \int_{t_{00}}^{t_f} \psi(t) \frac{\partial f}{\partial y}(t, x_0(t), x_0(t - \tau_0)) \dot{x}_0(t - \tau_0) dt \right\} \delta \tau + \pi \delta \xi \le 0, \\ & \delta v \in \text{cone}(V_{v_0} - v_0). \end{aligned} \quad (4.7) $$

The condition $\delta v \in \text{cone}(V_{v_0} - v_0)$ is equivalent to conditions:

$$ \delta\xi \in \mathbb{R}_{+}, \delta t_0 \in \mathbb{R}, \delta\tau \in \mathbb{R}. $$

Let $\delta t_0 = \delta \tau = 0$ in (4.7), then we obtain

$$ \pi \delta \xi \le 0, \forall \delta \xi \in \mathbb{R}_{+}. $$

This implies $\pi < 0$, therefore we can take $\pi = -1$, i.e., $\psi$ satisfies the condition (3.2).

Let $\delta\tau = \delta\xi = 0$ then, taking into account that $\delta t_0 \in \mathbb{R}$ from (4.7) we obtain the condition (1).

Let $\delta t_0 = \delta \xi = 0$ then, taking into account that $\delta \tau \in \mathbb{R}$ from (4.7) we obtain the condition (2).

Finally, we note that in order to prove the Theorem 3.3 it suffices to replace the set $V_{v_0}$ by the set

$$ [0, \alpha) \times [s_1, s_1 + \alpha) \times [\theta_1, \theta_0 + \alpha) $$

see the proof of the Theorem 3.1.

## 5. NUMERICAL APPROACH

In this section we shall present the numerical approach to solve the problem (2.3). First of all, we replace the constrained problem by its appropriate discrete optimal control problem. We subdivide the interval $[t_0, t_f]$ into $N$ sub-intervals with knots $t_0 < t_1 < \dots < t_f$ and $t_k = k\Delta t_k$, where $\Delta t_k$ is the mesh size of $k$th sub-interval. Where $k \in \{0, 1, \dots, N\}$.

The discrete problem associates to the continuous problem (2.3) is given by

$$ \left\{ \begin{array}{l} \min J(W, x_N) \\ x_{k+1} = \Psi(\Delta t_k, f(t_k, x_k, x_k(t_k - \tau))), \end{array} \right. \quad (5.1) $$

where $\Psi$ is a method of approximation of the delay differential equations (2.1) and (2.2). In this paper, the state problem will be solved by the solver dde23 using the popular approach based on the Runge-Kutta triple BS(2,3) [27].

### 5.1. Numerical algorithm

Evolutionary Algorithms are based on a model of natural, biological evolution, which was primarily developed by Holland [14]. It is essentially a searching method based on the Darwinian principles of biological evolution that explains the adaptive change of species by the principle of natural selection, which favors those species for survival and further evolution that are best adapted to their environmental conditions. It has been successfully applied to various optimization problems [6-8, 19].

In the evolutionary framework, a new generation of individuals is produced using the simulated genetic operations selection, crossover, and mutation. The evolution is based only on fitness. The fitness of an individual is measured only indirectly by its growth rate in comparison to others. Selection is the ability of those individuals that have outlasted the struggle for existence to bring their genetic information to the next generation. The crossover is the mechanism of reproduction based on the combination of the selected individuals. While the mutation is characterized by modifying one or more genes of an individual with some probability. The object of the mutation operator is to restore the lost or unexplored genetic material into the population to stop the premature convergence of the evolutionary algorithm to sub-optimal solutions.
The evolutionary algorithms differ from other methods of search and optimization in several practices. (1) Evolutionary algorithms, search from a population of possible solutions instead of a single one. (2) The fitness or cost function used to solve the redundancy has no necessity for continuity of the derivatives, so practically any fitness function can be chosen for optimizing. (3) Evolutionary algorithms use random operators across the process, including reproduction, crossover, and mutation. (4) Evolutionary algorithms are blind since no specified information about the intended problem is required to get the final solution.

The disadvantage of EA is the high cost of function evaluations, notably for large structures and larger population sizes, where many design iterations are necessary to reach the global optimum.

To circumvent this drawback, an improved evolutionary algorithm is considered to increase the performance of the optimal control problem based on parallel computing. Indeed, as it's well known, the EA readily lends itself to parallel processing because each individual in the population is independent of the others and therefore can be processed concurrently.

Furthermore, many parallel evolutionary algorithm implementations for solving optimization problems have demonstrated their success [21–24] as seen in the literature. For the implementation, we use a parallel simulation executed on a multiprocessor machine with distributed memory. There are various manners to exploit the parallelism of the EA algorithm. We opt for the implementation based on the master-slave approach. The interprocessor communication was realized via standard MPI (Message Passing Interface) [16]. In this procedure, each processor can directly access its memory and must explicitly communicate with other processors to access the data in their memories [16]. MPI's objective is to give a standard for writing message-passing programs. The ed algorithm is summarized in the Algorithm 1. Denotes by $m$ the size of the unknown parameter vector $W$.

**Algorithm 1. Parallel EA.**

1. Initialize MPI environment and configuration parameters.

2. Calculate $np$, the number of individuals to be evaluated by each processor.

3. Choose the bounds $x_i \in \mathbb{R}^m$ and $x_f \in \mathbb{R}^m$.

4. $t \leftarrow 0$

The master processor performs the following task: generate randomly an initial population $P(t)$, where all individual $W_i(t)$ are between $x_i$ and $x_f$.

5. $t \leftarrow t + 1$

Each processor performs the following task:

for each element $W_i(t)$ of the population $P(t)$ with $i = 1, \dots, np$ do:

(a) Given $W_i(t)$, Solve the state equation.

(b) Compute $J(W_i(t), x_i(t))$

The calculated data are exchanged by the MPI protocol.

6. If $|J(W_{best}(t), x_{best}(t))| < \epsilon$, then go to 8.

7. Run the genetic operation by the master processor:

(a) Selection

(b) Crossover

(c) Mutation

8. Collect all computed results and stop the MPI environment.

## 5.2. Numerical experiments

In order to establish the efficiency of the proposed approach. We consider some concrete problems governed by nonlinear DDE in the form (2.1). A good test problem must have an analytical solution. Finding attractive and nontrivial DDEs whose solution is known is sometimes a hard task. Each of the following examples emphasizes one or more of the peculiar aspects of DDEs that make their numerical solution challenging. In all sequel examples, we run the Algorithm 1 using the parameters listed in Table 1. All tests are executed under the same
The evolutionary algorithms differ from other methods of search and optimization in several practices. (1) Evolutionary algorithms, search from a population of possible solutions instead of a single one. (2) The fitness or cost function used to solve the redundancy has no necessity for continuity of the derivatives, so practically any fitness function can be chosen for optimizing. (3) Evolutionary algorithms use random operators across the process, including reproduction, crossover, and mutation. (4) Evolutionary algorithms are blind since no specified information about the intended problem is required to get the final solution.

The disadvantage of EA is the high cost of function evaluations, notably for large structures and larger population sizes, where many design iterations are necessary to reach the global optimum.

To circumvent this drawback, an improved evolutionary algorithm is considered to increase the performance of the optimal control problem based on parallel computing. Indeed, as it's well known, the EA readily lends itself to parallel processing because each individual in the population is independent of the others and therefore can be processed concurrently.

Furthermore, many parallel evolutionary algorithm implementations for solving optimization problems have demonstrated their success [21–24] as seen in the literature. For the implementation, we use a parallel simulation executed on a multiprocessor machine with distributed memory. There are various manners to exploit the parallelism of the EA algorithm. We opt for the implementation based on the master-slave approach. The interprocessor communication was realized via standard MPI (Message Passing Interface) [16]. In this procedure, each processor can directly access its memory and must explicitly communicate with other processors to access the data in their memories [16]. MPI's objective is to give a standard for writing message-passing programs. The ed algorithm is summarized in the Algorithm 1. Denotes by $m$ the size of the unknown parameter vector $W$.

**Algorithm 1. Parallel EA.**

1. Initialize MPI environment and configuration parameters.

2. Calculate $np$, the number of individuals to be evaluated by each processor.

3. Choose the bounds $x_i \in \mathbb{R}^m$ and $x_f \in \mathbb{R}^m$.

4. $t \leftarrow 0$

The master processor performs the following task: generate randomly an initial population $P(t)$, where all individual $W_i(t)$ are between $x_i$ and $x_f$.

5. $t \leftarrow t + 1$

Each processor performs the following task:

for each element $W_i(t)$ of the population $P(t)$ with $i = 1, \dots, np$ do:

(a) Given $W_i(t)$, Solve the state equation.

(b) Compute $J(W_i(t), x_i(t))$

The calculated data are exchanged by the MPI protocol.

6. If $|J(W_{best}(t), x_{best}(t))| < \epsilon$, then go to 8.

7. Run the genetic operation by the master processor:

(a) Selection

(b) Crossover

(c) Mutation

8. Collect all computed results and stop the MPI environment.

## 5.2. Numerical experiments

In order to establish the efficiency of the proposed approach. We consider some concrete problems governed by nonlinear DDE in the form (2.1). A good test problem must have an analytical solution. Finding attractive and nontrivial DDEs whose solution is known is sometimes a hard task. Each of the following examples emphasizes one or more of the peculiar aspects of DDEs that make their numerical solution challenging. In all sequel examples, we run the Algorithm 1 using the parameters listed in Table 1. All tests are executed under the same
<table><caption>TABLE 1. The (EA) parameters.</caption><tr><td>Population size N</td><td>16</td></tr><tr><td>Selection type</td><td>Roulette-wheel</td></tr><tr><td>Crossover type</td><td>Barycentric with random coefficients</td></tr><tr><td>Crossover probability P<sub>c</sub></td><td>0.6</td></tr><tr><td>Type of mutation</td><td>Non uniform</td></tr><tr><td>Mutation probability P<sub>m</sub></td><td>0.5%</td></tr></table>

machine footnote Calculator SGI, cluster-type calculation, consists of 888 Xeon cores clocked at 2.66 GHz with 2 GB of memory per core (http://www.ccipl.univ-nantes.fr), using the same solver.

### Example 1

We consider the following with constant and multiple delays. For $t_0 \le t \le 1$

$$ \left\{ \begin{aligned} \dot{x}_1(t) &= x_5(t - \tau_1) + x_3(t - \tau_1) \\ \dot{x}_2(t) &= x_1(t - \tau_1) + x_2(t - \tau_2) \\ \dot{x}_3(t) &= x_3(t - \tau_1) + x_1(t - \tau_2) \\ \dot{x}_4(t) &= x_5(t - \tau_1)x_4(t - \tau_1) \\ \dot{x}_5(t) &= x_1(t - \tau_1), \end{aligned} \right. \qquad (5.2) $$

with history

$$ \begin{array}{lll} x_1(t) = \exp(t+1), & x_2(t) = \exp(t+0.5), & x_3(t) = \sin(t+1), \\ x_4(t) = x_1(t), & x_5(t) = x_1(t) & \text{for } t \le t_0. \end{array} $$

In this case, we shall find $x_1(t)$, $x_2(t)$, $x_3(t)$, $x_4(t)$, $x_5(t)$ for $0 \le t \le 1$ and the control parameter $W = (t_0, \tau_1, \tau_2)$.

The desired parameter is

$$ W = (t_0, \tau_1, \tau_2) = (0, 1, 0.5) $$

and the exact solution is

$$ \left\{ \begin{aligned} \dot{x}_1(t) &= \exp(t) - \cos(t) + e \quad 0 \le t \le 1 \\ \dot{x}_2(t) &= \begin{cases} 2\exp(t) + \exp(0.5) - 2 & 0 \le t \le 0.5 \\ \exp(t) + 2(t-5) + t\exp(0.5) - 2t \\ +1.5\exp(0.5) - 3 & 0.5 \le t \le 1 \end{cases} \\ \dot{x}_3(t) &= \begin{cases} \exp(t+0.5) - \cos(t) + 1 - \exp(0.5) + \sin(1), & 0 \le t \le 0.5 \\ -\cos(t) + \exp(t-0.5) - \sin(t-0.5) + (t+0.5)e \\ -\exp(0.5) + \sin(1), & 0.5 \le t \le 1 \end{cases} \\ \dot{x}_4(t) &= 0.5 \exp(2t) - 0.5 + e, \quad 0 \le t \le 1 \\ \dot{x}_5(t) &= \exp(t) + e - 1, \quad 0 \le t \le 1. \end{aligned} \right. \qquad (5.3) $$

In Figure 1, we present the comparison of the obtained solution using the Algorithm 1.

The evolution of the cost function is presented in the Figure 2.

In the Table 2 we present the comparison between the obtained values of the parameter $W = (t_0, \tau_1, \tau_2)$, and the exact one.

As we can see from these results the obtained approximation is in good agreement with the exact solution.
<table><caption>TABLE 1. The (EA) parameters.</caption><tr><td>Population size N</td><td>16</td></tr><tr><td>Selection type</td><td>Roulette-wheel</td></tr><tr><td>Crossover type</td><td>Barycentric with random coefficients</td></tr><tr><td>Crossover probability P<sub>c</sub></td><td>0.6</td></tr><tr><td>Type of mutation</td><td>Non uniform</td></tr><tr><td>Mutation probability P<sub>m</sub></td><td>0.5%</td></tr></table>

machine footnote Calculator SGI, cluster-type calculation, consists of 888 Xeon cores clocked at 2.66 GHz with 2 GB of memory per core (http://www.ccipl.univ-nantes.fr), using the same solver.

### Example 1

We consider the following with constant and multiple delays. For $t_0 \le t \le 1$

$$ \left\{ \begin{aligned} \dot{x}_1(t) &= x_5(t - \tau_1) + x_3(t - \tau_1) \\ \dot{x}_2(t) &= x_1(t - \tau_1) + x_2(t - \tau_2) \\ \dot{x}_3(t) &= x_3(t - \tau_1) + x_1(t - \tau_2) \\ \dot{x}_4(t) &= x_5(t - \tau_1)x_4(t - \tau_1) \\ \dot{x}_5(t) &= x_1(t - \tau_1), \end{aligned} \right. \qquad (5.2) $$

with history

$$ \begin{array}{lll} x_1(t) = \exp(t+1), & x_2(t) = \exp(t+0.5), & x_3(t) = \sin(t+1), \\ x_4(t) = x_1(t), & x_5(t) = x_1(t) & \text{for } t \le t_0. \end{array} $$

In this case, we shall find $x_1(t)$, $x_2(t)$, $x_3(t)$, $x_4(t)$, $x_5(t)$ for $0 \le t \le 1$ and the control parameter $W = (t_0, \tau_1, \tau_2)$.

The desired parameter is

$$ W = (t_0, \tau_1, \tau_2) = (0, 1, 0.5) $$

and the exact solution is

$$ \left\{ \begin{aligned} \dot{x}_1(t) &= \exp(t) - \cos(t) + e \quad 0 \le t \le 1 \\ \dot{x}_2(t) &= \begin{cases} 2\exp(t) + \exp(0.5) - 2 & 0 \le t \le 0.5 \\ \exp(t) + 2(t-5) + t\exp(0.5) - 2t \\ +1.5\exp(0.5) - 3 & 0.5 \le t \le 1 \end{cases} \\ \dot{x}_3(t) &= \begin{cases} \exp(t+0.5) - \cos(t) + 1 - \exp(0.5) + \sin(1), & 0 \le t \le 0.5 \\ -\cos(t) + \exp(t-0.5) - \sin(t-0.5) + (t+0.5)e \\ -\exp(0.5) + \sin(1), & 0.5 \le t \le 1 \end{cases} \\ \dot{x}_4(t) &= 0.5\exp(2t) - 0.5 + e, \quad 0 \le t \le 1 \\ \dot{x}_5(t) &= \exp(t) + e - 1, \quad 0 \le t \le 1. \end{aligned} \right. \qquad (5.3) $$

In Figure 1, we present the comparison of the obtained solution using the Algorithm 1.

The evolution of the cost function is presented in the Figure 2.

In the Table 2 we present the comparison between the obtained values of the parameter $W = (t_0, \tau_1, \tau_2)$, and the exact one.

As we can see from these results the obtained approximation is in good agreement with the exact solution.
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Comparison of the approximate and exact solution of the problem (5.2).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 2. Evolution of the cost function.</figcaption></figure>

<table><caption>TABLE 2. Comparison between the exact and the approximate parameter $W = (t_0, \tau_1, \tau_2)$.</caption><thead><tr><th></th><th>t<sub>0</sub></th><th>τ<sub>1</sub></th><th>τ<sub>2</sub></th></tr></thead><tbody><tr><td>Exact</td><td>0</td><td>1</td><td>0.5</td></tr><tr><td>Approximate</td><td>0.0008</td><td>0.99999</td><td>0.49999</td></tr></tbody></table>

### Example 2

In this example, we consider the model of the four-year life cycle of a population of lemmings [30]. The problem is defined as follows

$$ \begin{cases} \dot{x}(t) = rx(t)\left(1 - \frac{x(t-\tau)}{m}\right), & t_0 \le t \le 40 \\ x(t) = 19, & t \le t_0. \end{cases} \qquad (5.4) $$

The parameters $r = 3.5$ and $m = 19$. Note that with these values, the problem has a constant (steady-state) solution, $x(t) = 19$. The original model is proposed with the solution as history and perturbs the initial value to $x(0) = 19.00001$ so that $x(t)$ will move away from the steady state. Here we use $x(0) = 19.001$ so as to find out the cyclic behavior sooner. The exact control parameters in this case are $W = (\tau, t_0) = (0.74, 0)$.

In the Figure 3, we present the obtained solution compared to the exact one.
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Comparison of the approximate and exact solution of the problem (5.2).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 2. Evolution of the cost function.</figcaption></figure>

<table><caption>TABLE 2. Comparison between the exact and the approximate parameter $W = (t_0, \tau_1, \tau_2)$.</caption><thead><tr><th></th><th>t<sub>0</sub></th><th>τ<sub>1</sub></th><th>τ<sub>2</sub></th></tr></thead><tbody><tr><td>Exact</td><td>0</td><td>1</td><td>0.5</td></tr><tr><td>Approximate</td><td>0.0008</td><td>0.99999</td><td>0.49999</td></tr></tbody></table>

### Example 2

In this example, we consider the model of the four-year life cycle of a population of lemmings [30]. The problem is defined as follows

$$ \begin{cases} \dot{x}(t) = rx(t)\left(1 - \frac{x(t-\tau)}{m}\right), & t_0 \le t \le 40 \\ x(t) = 19, & t \le t_0. \end{cases} \qquad (5.4) $$

The parameters $r = 3.5$ and $m = 19$. Note that with these values, the problem has a constant (steady-state) solution, $x(t) = 19$. The original model is proposed with the solution as history and perturbs the initial value to $x(0) = 19.00001$ so that $x(t)$ will move away from the steady state. Here we use $x(0) = 19.001$ so as to find out the cyclic behavior sooner. The exact control parameters in this case are $W = (\tau, t_0) = (0.74, 0)$.

In the Figure 3, we present the obtained solution compared to the exact one.
<figure><img src="image_5.png" /><figcaption>FIGURE 3. The comparison between the approximate and the exact solution of the problem (5.4).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 4. The evolution of the cost function.</figcaption></figure>

<table><caption>TABLE 3. Comparison between the exact and the approximate parameter $W = (t_0, \tau)$.</caption><thead><tr><th></th><th>t<sub>0</sub></th><th>&tau;</th></tr></thead><tbody><tr><td>Exact</td><td>0</td><td>0.74</td></tr><tr><td>Approximate</td><td>0.00007</td><td>0.741345</td></tr></tbody></table>

In the Figure 4, we present the evolution of the cost function related to the iteration.

In the Table 3 we present the comparison between the obtained values and the exact parameters.

We can see from the Figures 3 and 4 and the Table 3, that the obtained results show the robustness of the proposed method despite the difficulty of the problem which presents a very high sensitivity to perturbation.

### Example 3

In this example, we consider the predator-prey models derived by introducing a resource limitation on the prey and assuming the birth rate of predators respond to adjustment in the magnitude of the population $x_1$ of prey and the population $x_2$ of predators just after a time delay $\tau$.

$$ \begin{cases} \dot{x}_1(t) = ax_1(t)\left(1 - \frac{x_1(t)}{m}\right) + bx_1(t)x_2(t), & t_0 \le t \le 100 \\ \dot{x}_2(t) = cx_2(t) + dx_1(t - \tau)x_2(t - \tau), & t_0 \le t \le 100 \\ x_1(t) = 80, \quad x_2(t) = 30, & t \le t_0. \end{cases} \quad (5.5) $$
<figure><img src="image_5.png" /><figcaption>FIGURE 3. The comparison between the approximate and the exact solution of the problem (5.4).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 4. The evolution of the cost function.</figcaption></figure>

<table><caption>TABLE 3. Comparison between the exact and the approximate parameter $W = (t_0, \tau)$.</caption><thead><tr><th></th><th>t<sub>0</sub></th><th>&tau;</th></tr></thead><tbody><tr><td>Exact</td><td>0</td><td>0.74</td></tr><tr><td>Approximate</td><td>0.00007</td><td>0.741345</td></tr></tbody></table>

In the Figure 4, we present the evolution of the cost function related to the iteration.

In the Table 3 we present the comparison between the obtained values and the exact parameters.

We can see from the Figures 3 and 4 and the Table 3, that the obtained results show the robustness of the proposed method despite the difficulty of the problem which presents a very high sensitivity to perturbation.

### Example 3

In this example, we consider the predator-prey models derived by introducing a resource limitation on the prey and assuming the birth rate of predators respond to adjustment in the magnitude of the population $x_1$ of prey and the population $x_2$ of predators just after a time delay $\tau$.

$$ \begin{cases} \dot{x}_1(t) = ax_1(t)\left(1 - \frac{x_1(t)}{m}\right) + bx_1(t)x_2(t), & t_0 \le t \le 100 \\ \dot{x}_2(t) = cx_2(t) + dx_1(t - \tau)x_2(t - \tau), & t_0 \le t \le 100 \\ x_1(t) = 80, \quad x_2(t) = 30, & t \le t_0. \end{cases} \qquad (5.5) $$
<figure><img src="image_7.png" /><figcaption>FIGURE 5. The population $x_2$ of predators related to population $x_1$ of prey: comparison between the exact and the approximate solution.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 6. Evolution of the cost function related to the iterations.</figcaption></figure>

<table><caption>TABLE 4. The exact and the approximate parameter $W = (t_0, \tau)$.</caption><thead><tr><th></th><th>t<sub>0</sub></th><th>&tau;</th></tr></thead><tbody><tr><td>Exact</td><td>0</td><td>1</td></tr><tr><td>Approximate</td><td>0.0135</td><td>0.9994</td></tr></tbody></table>

Suppose that the parameters $a = 0.25$, $b = -0.01$, $c = -1.00$, $d = 0.01$, and $m = 200$. In the Figure 5 we present the behavior of the population of predators $x_2(t)$ according to the population of prey $x_1(t)$ for $t_0 \le t \le 100$. In the Figure 6, we present the evolution of the cost function related to the iterations. The Table 4 presents the obtained values of parameters and the exact parameters.

This model is more complicated to the previous. This explains the little lack of precision in the obtained results. However, the results still in good agreement with the exact solution.

### Example 4

In this last example, we consider a cardiovascular model due to Ottesen [25] concerning the arterial pressure, $Pa(t) = x_1(t)$, the venous pressure, $Pv(t) = x_2(t)$, and the heart rate, $H(t) = x_3(t)$. Ottesen study conditions under which the delay causes qualitative differences in the solution and in particular, oscillations in $Pa(t)$. The
<figure><img src="image_7.png" /><figcaption>FIGURE 5. The population $x_2$ of predators related to population $x_1$ of prey: comparison between the exact and the approximate solution.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 6. Evolution of the cost function related to the iterations.</figcaption></figure>

<table><caption>TABLE 4. The exact and the approximate parameter $W = (t_0, \tau)$.</caption><thead><tr><th></th><th>t<sub>0</sub></th><th>&tau;</th></tr></thead><tbody><tr><td>Exact</td><td>0</td><td>1</td></tr><tr><td>Approximate</td><td>0.0135</td><td>0.9994</td></tr></tbody></table>

Suppose that the parameters $a = 0.25$, $b = -0.01$, $c = -1.00$, $d = 0.01$, and $m = 200$. In the Figure 5 we present the behavior of the population of predators $x_2(t)$ according to the population of prey $x_1(t)$ for $t_0 \le t \le 100$. In the Figure 6, we present the evolution of the cost function related to the iterations. The Table 4 presents the obtained values of parameters and the exact parameters.

This model is more complicated to the previous. This explains the little lack of precision in the obtained results. However, the results still in good agreement with the exact solution.

### Example 4

In this last example, we consider a cardiovascular model due to Ottesen [25] concerning the arterial pressure, $Pa(t) = x_1(t)$, the venous pressure, $Pv(t) = x_2(t)$, and the heart rate, $H(t) = x_3(t)$. Ottesen study conditions under which the delay causes qualitative differences in the solution and in particular, oscillations in $Pa(t)$. The
<table><caption>TABLE 5. Parameters of the cardiovascular model.</caption><thead><tr><th>Parameter</th><th>R</th><th>r</th><th>V<sub>str</sub></th><th>c<sub>a</sub></th><th>c<sub>v</sub></th><th>α<sub>s</sub></th><th>β<sub>s</sub></th><th>α<sub>p</sub></th><th>β<sub>p</sub></th></tr></thead><tbody><tr><td>Value</td><td>1.05</td><td>0.068</td><td>67.9</td><td>1.55</td><td>519</td><td>93</td><td>7</td><td>93</td><td>7</td></tr><tr><td>Parameter</td><td>α<sub>H</sub></td><td>γ<sub>H</sub></td><td>β<sub>H</sub></td><td>P<sub>0</sub></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Value</td><td>0.84</td><td>0</td><td>1.17</td><td>93</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

<table><caption>TABLE 6. The exact and the approximate parameter $W = (t_0, \tau)$.</caption><thead><tr><th></th><th>t<sub>0</sub></th><th>τ</th></tr></thead><tbody><tr><td>Exact</td><td>0</td><td>1</td></tr><tr><td>Approximate</td><td>-2.4384e-006</td><td>1.00000006</td></tr></tbody></table>

model is presented as follows

$$ 
\begin{aligned}
x_1(t) &= \frac{1}{c_a} \left( \frac{1}{R} (-x_1(t) + x_2(t)) + V_{str}x_3(t) \right), & t_0 \le t \le 350 \\
x_2(t) &= \frac{1}{c_v} \left( \frac{1}{R} x_1(t) - \left( \frac{1}{R} + \frac{1}{r} \right) x_2(t) \right), & t_0 \le t \le 350 \\
x_3(t) &= f(T_s, T_p), & t_0 \le t \le 350,
\end{aligned} 
\quad (5.6) $$

where

$$ 
\begin{aligned}
T_s &= \frac{1}{1 + \left( \frac{x_1(t-\tau)}{\alpha_s} \right)^{\beta_s}} \\
T_p &= \frac{1}{1 + \left( \frac{\alpha_p}{x_1(t)} \right)^{\beta_p}} \\
f(T_s, T_p) &= \frac{\alpha_H T_s}{1 + \gamma_H T_p} - \beta_H T_p.
\end{aligned} 
$$

For the history is given by

$$ 
\begin{aligned}
x_1(t) &= P_0, & t \le t_0 \\
x_2(t) &= \frac{1}{1 + \frac{R}{r}} P_0, & t \le t_0 \\
x_3(t) &= \frac{1}{RV_{str}\left(1 + \frac{r}{R}\right)}, & t \le t_0.
\end{aligned} 
$$

The used parameters are given in the Table 5.

For the sake of visibility, we will limit ourselves to the presentation of the behavior of the arterial pressure $P_a(t)$ in the Figure 7. The other obtained results has a similar quality.

The evolution of the cost function for this example is presented in the Figure 8.

The approximation of the free parameter $\tau$ and $t_0$ is presented in the Table 6.

We have shown the efficiency of the proposed approach through different problems with different kinds of difficulties. According to the presented figures and tables, the obtained results confirm the robustness of the
<table><caption>TABLE 5. Parameters of the cardiovascular model.</caption><thead><tr><th>Parameter</th><th>R</th><th>r</th><th>V<sub>str</sub></th><th>c<sub>a</sub></th><th>c<sub>v</sub></th><th>α<sub>s</sub></th><th>β<sub>s</sub></th><th>α<sub>p</sub></th><th>β<sub>p</sub></th></tr></thead><tbody><tr><td>Value</td><td>1.05</td><td>0.068</td><td>67.9</td><td>1.55</td><td>519</td><td>93</td><td>7</td><td>93</td><td>7</td></tr><tr><td>Parameter</td><td>α<sub>H</sub></td><td>γ<sub>H</sub></td><td>β<sub>H</sub></td><td>P<sub>0</sub></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Value</td><td>0.84</td><td>0</td><td>1.17</td><td>93</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

<table><caption>TABLE 6. The exact and the approximate parameter $W = (t_0, \tau)$.</caption><thead><tr><th></th><th>t<sub>0</sub></th><th>τ</th></tr></thead><tbody><tr><td>Exact</td><td>0</td><td>1</td></tr><tr><td>Approximate</td><td>-2.4384e-006</td><td>1.00000006</td></tr></tbody></table>

model is presented as follows

$$ 
\begin{aligned}
x_1(t) &= \frac{1}{c_a} \left( \frac{1}{R} (-x_1(t) + x_2(t)) + V_{str}x_3(t) \right), & t_0 \le t \le 350 \\
x_2(t) &= \frac{1}{c_v} \left( \frac{1}{R} x_1(t) - \left( \frac{1}{R} + \frac{1}{r} \right) x_2(t) \right), & t_0 \le t \le 350 \\
x_3(t) &= f(T_s, T_p), & t_0 \le t \le 350,
\end{aligned} 
\quad (5.6) $$

where

$$ 
\begin{aligned}
T_s &= \frac{1}{1 + \left( \frac{x_1(t-\tau)}{\alpha_s} \right)^{\beta_s}} \\
T_p &= \frac{1}{1 + \left( \frac{\alpha_p}{x_1(t)} \right)^{\beta_p}} \\
f(T_s, T_p) &= \frac{\alpha_H T_s}{1 + \gamma_H T_p} - \beta_H T_p.
\end{aligned} 
$$

For the history is given by

$$ 
\begin{aligned}
x_1(t) &= P_0, & t \le t_0 \\
x_2(t) &= \frac{1}{1 + \frac{R}{r}} P_0, & t \le t_0 \\
x_3(t) &= \frac{1}{RV_{str}\left(1 + \frac{r}{R}\right)}, & t \le t_0.
\end{aligned} 
$$

The used parameters are given in the Table 5.

For the sake of visibility, we will limit ourselves to the presentation of the behavior of the arterial pressure $P_a(t)$ in the Figure 7. The other obtained results has a similar quality.

The evolution of the cost function for this example is presented in the Figure 8.

The approximation of the free parameter $\tau$ and $t_0$ is presented in the Table 6.

We have shown the efficiency of the proposed approach through different problems with different kinds of difficulties. According to the presented figures and tables, the obtained results confirm the robustness of the
<figure><img src="image_9.png" /><figcaption>FIGURE 7. The behavior of the arterial pressure $P_a(t)$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 8. The evolution of the cost function.</figcaption></figure>

developed approach in terms of precision quality. In the following section, we will study the performance of the proposed algorithm in terms of time execution.

## 5.3. Performance of the parallel Algorithm 1 in term of execution time

In order to examine the performance of the parallel algorithm, we start by analyzing the time complexity of the sequential evolutionary algorithm. The sequential Complexity is given by

<table><thead><tr><th>Complexity of sequential EA</th><th>Complexity parallel EA</th></tr></thead><tbody><tr><td>O(N * ng * O(J<sub>fitness</sub>))</td><td>O(N * ng * <sup>O(J<sub>fitness</sub>)</sup>/<sub>np</sub>) + CC</td></tr></tbody></table>

where ng denotes the number of iterations, np denotes the number of used processors and CC is the communication cost. Comparison of consuming time-related to the number of processors np is not enough to analyze the performance of the parallel algorithm. Indeed, the communication cost may offset any gains in computation time. Therefore, to ensure that the parallel implementation has a better performance than a simple one we must
<figure><img src="image_9.png" /><figcaption>FIGURE 7. The behavior of the arterial pressure $P_a(t)$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 8. The evolution of the cost function.</figcaption></figure>

developed approach in terms of precision quality. In the following section, we will study the performance of the proposed algorithm in terms of time execution.

## 5.3. Performance of the parallel Algorithm 1 in term of execution time

In order to examine the performance of the parallel algorithm, we start by analyzing the time complexity of the sequential evolutionary algorithm. The sequential Complexity is given by

<table><thead><tr><th>Complexity of sequential EA</th><th>Complexity parallel EA</th></tr></thead><tbody><tr><td>O(N * ng * O(J<sub>fitness</sub>))</td><td>O(N * ng * <sup>O(J<sub>fitness</sub>)</sup>/<sub>np</sub>) + CC</td></tr></tbody></table>

where ng denotes the number of iterations, np denotes the number of used processors and CC is the communication cost. Comparison of consuming time-related to the number of processors np is not enough to analyze the performance of the parallel algorithm. Indeed, the communication cost may offset any gains in computation time. Therefore, to ensure that the parallel implementation has a better performance than a simple one we must
<figure><img src="image_11.png" /><figcaption>FIGURE 9. Time consuming with different used number of processors *np*.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 10. Parallel speedup and the efficiency related to the used number of processors.</figcaption></figure>

take into consideration the parallel speedup *S* and the efficiency defined as follows

$$S = \frac{T_s}{T_p}, \quad E = \frac{S}{np} \qquad (5.7)$$

where $T_s$ denotes the times of sequential algorithm and $T_p$ is the one of the parallel algorithm.

In the Figure 9, we present the behavior of the time consuming related to the used number of processor.

The parallel speedup *S* and the efficiency *E* related the used number of processors is presented in the Figure 10.

According to the presented results in Figures 9 and 10, we can see that the parallel version reduces considerably the execution time in comparison to the sequential version. Indeed, only passing from a single processor to two processors reduces the computing time by half. Obviously, the ideal would be that the time observed using *np* processors should be exactly the time observed using 1 processor divided by *np*. However, in practice, the difference increases with the number of processors. This is due to the communication time, which increases with the number of processors. This is confirmed by noting the passage from 8 to 16 processors which did not save time. On the contrary, the time has intensified. This is explained by the fact that communication time has become more important than the gain from the distribution of tasks. Thus, from Figure 10 we can see that the parallel speedup increases until 8 processors, then stagnates. While the efficiency we can observe that it
<figure><img src="image_11.png" /><figcaption>FIGURE 9. Time consuming with different used number of processors *np*.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 10. Parallel speedup and the efficiency related to the used number of processors.</figcaption></figure>

take into consideration the parallel speedup *S* and the efficiency defined as follows

$$S = \frac{T_s}{T_p}, \quad E = \frac{S}{np} \qquad (5.7)$$

where $T_s$ denotes the times of sequential algorithm and $T_p$ is the one of the parallel algorithm.

In the Figure 9, we present the behavior of the time consuming related to the used number of processor.

The parallel speedup *S* and the efficiency *E* related the used number of processors is presented in the Figure 10.

According to the presented results in Figures 9 and 10, we can see that the parallel version reduces considerably the execution time in comparison to the sequential version. Indeed, only passing from a single processor to two processors reduces the computing time by half. Obviously, the ideal would be that the time observed using *np* processors should be exactly the time observed using 1 processor divided by *np*. However, in practice, the difference increases with the number of processors. This is due to the communication time, which increases with the number of processors. This is confirmed by noting the passage from 8 to 16 processors which did not save time. On the contrary, the time has intensified. This is explained by the fact that communication time has become more important than the gain from the distribution of tasks. Thus, from Figure 10 we can see that the parallel speedup increases until 8 processors, then stagnates. While the efficiency we can observe that it
decreases as the number of processors increases. The critical number of slaves that maintain a predetermined efficiency and maximize the parallel speedup is 8.

## 6. CONCLUSION

In this work, we have proposed an optimal control approach for solving an inverse problem governed by a nonlinear delay differential equation. Based on some observed data we have determined the initial moment and delay parameter. Thus, the necessary conditions of optimality related to the associate control problem have been formulated and proved. To approximate the obtained optimal control problem we have discussed an approach based on an improved parallel evolutionary algorithm. Through some numerical experiments, we have studied the performance of the proposed algorithm. In particular, a useful analysis of the parallel algorithms is presented. We claim that the proposed approach can easily be adapted to solve other inverse problems of parameter identification.

*Acknowledgements.* The paper was developed when the first author was visiting professor at the Laboratory of Mathematics Jean Leray UMR6629 CNRS, Nantes University, in November-December 2021. All parallel programs were executed by using the calculator SGI of the "Centre de Calcul Intensif des Pays de la Loire (CCIPL): http://www.ccipl.univ-nantes.fr".

## REFERENCES

[1] L. Afraites, A. Hadri, A. Laghrib and M. Nachaoui, A weighted parameter identification pde-constrained optimization for inverse image denoising problem. The Visual Computer (2021) 1–16. DOI: 10.1007/s00371-021-02162-x.

[2] L. Afraites, C. Masnaoui and M. Nachaoui, Shape optimization method for an inverse geometric source problem and stability at critical shape. *Discrete Contin. Dyn. Syst. Ser. S* **15** (2022) 1–21.

[3] K.A. Berdawood, A. Nachaoui, M. Nachaoui and F. Aboud, An effective relaxed alternating procedure for cauchy problem connected with helmholtz equation. *Numer. Methods Part. Differ. Equ.* (2021). DOI: 10.1002/num.22793.

[4] G. Bocharov and A. Romanyukha, Mathematical modeling of the immune response during acute viral infections. In: Theoretical and Experimental Insights into Immunology. Springer (1992) 309–321.

[5] G. Bocharov and A. Romanyukha, Numerical treatment of the parameter identification problem for delay-differential systems arising in immune response modelling. *Appl. Numer. Math.* **15** (1994) 307–326.

[6] A. Chakib, A. Nachaoui and M. Nachaoui, Finite element approximation of an optimal design problem. *J. Appl. Comput. Math.* **11** (2012) 19–26.

[7] A. Chakib, A. Ellabib, A. Nachaoui and M. Nachaoui, A shape optimization formulation of weld pool determination. *Appl. Math. Lett.* **25** (2012) 374–379.

[8] A. Chakib, A. Nachaoui and M. Nachaoui, Approximation and numerical realization of an optimal design welding problem. *Numer. Methods Part. Differ. Equ.* **29** (2013) 1563–1586.

[9] P. Dvalishvili, A. Nachaoui and T. Tadumadze, Effects of the initial moment and several delays perturbations in the variation formulas for a solution of a functional differential equation with the continuous initial condition. *Georgian Math. J.* **27** (2020) 53–66.

[10] J. Eastman, J. Sass, J.M. Gomes, R.W. Dos Santos and E.M. Cherry, Using delay differential equations to induce alternans in a model of cardiac electrophysiology. *J. Theor. Biol.* **404** (2016) 262–272.

[11] J.E. Forde, *Delay differential equation models in mathematical biology*. Ph.D thesis, University of Michigan, ProQuest LLC, Ann Arbor, MI (2005).

[12] R.V. Gamkrelidze and G.L. Karatishvili, Extremal problems in linear topological spaces. *Math. Syst. Theory* **1** (1967) 229–256.

[13] J. Hadamard, Lectures on Cauchy's Problem in Linear Partial Differential Equations. Dover Publications, New York (1953).

[14] J.H. Holland, *Adaptation in Natural and Artificial Systems: An Introductory Analysis With Applications to Biology, Control, and Artificial Intelligence*. MIT Press (1992).

[15] M. Huang, J. Luo, L. Hu, B. Zheng and J. Yu, Assessing the efficiency of wolbachia driven aedes mosquito suppression by delay differential equations. *J. Theor. Biol.* **440** (2018) 1–11.

[16] J. Kepner, Parallel programming with MatlabMPI. Preprint astro-ph/0107406 (2001).

[17] G.L. Kharatishvili and T.A. Tadumadze, Formulas for the Variation of a Solution and Optimal Control Problems for Differential Equations with Retarded Arguments. *J. Math. Sci. (NY)* **140** (2007) 1–175.

[18] S. Lenhart, V. Protopopescu and J. Yong, Solving inverse problems of identification type by optimal control methods. I: AIP Conference Proceedings. Vol. 411. American Institute of Physics (1997) 87–94.

[19] M. Nachaoui, *Étude théorique et approximation numérique d'un problème inverse de transfert de la chaleur*. Ph.D. thesis, Université de Nantes (2011).
decreases as the number of processors increases. The critical number of slaves that maintain a predetermined efficiency and maximize the parallel speedup is 8.

## 6. CONCLUSION

In this work, we have proposed an optimal control approach for solving an inverse problem governed by a nonlinear delay differential equation. Based on some observed data we have determined the initial moment and delay parameter. Thus, the necessary conditions of optimality related to the associate control problem have been formulated and proved. To approximate the obtained optimal control problem we have discussed an approach based on an improved parallel evolutionary algorithm. Through some numerical experiments, we have studied the performance of the proposed algorithm. In particular, a useful analysis of the parallel algorithms is presented. We claim that the proposed approach can easily be adapted to solve other inverse problems of parameter identification.

*Acknowledgements.* The paper was developed when the first author was visiting professor at the Laboratory of Mathematics Jean Leray UMR6629 CNRS, Nantes University, in November-December 2021. All parallel programs were executed by using the calculator SGI of the "Centre de Calcul Intensif des Pays de la Loire (CCIPL): http://www.ccipl.univ-nantes.fr".

## REFERENCES

[1] L. Afraites, A. Hadri, A. Laghrib and M. Nachaoui, A weighted parameter identification pde-constrained optimization for inverse image denoising problem. The Visual Computer (2021) 1–16. DOI: 10.1007/s00371-021-02162-x.

[2] L. Afraites, C. Masnaoui and M. Nachaoui, Shape optimization method for an inverse geometric source problem and stability at critical shape. *Discrete Contin. Dyn. Syst. Ser. S* **15** (2022) 1–21.

[3] K.A. Berdawood, A. Nachaoui, M. Nachaoui and F. Aboud, An effective relaxed alternating procedure for cauchy problem connected with helmholtz equation. *Numer. Methods Part. Differ. Equ.* (2021). DOI: 10.1002/num.22793.

[4] G. Bocharov and A. Romanyukha, Mathematical modeling of the immune response during acute viral infections. In: Theoretical and Experimental Insights into Immunology. Springer (1992) 309–321.

[5] G. Bocharov and A. Romanyukha, Numerical treatment of the parameter identification problem for delay-differential systems arising in immune response modelling. *Appl. Numer. Math.* **15** (1994) 307–326.

[6] A. Chakib, A. Nachaoui and M. Nachaoui, Finite element approximation of an optimal design problem. *J. Appl. Comput. Math.* **11** (2012) 19–26.

[7] A. Chakib, A. Ellabib, A. Nachaoui and M. Nachaoui, A shape optimization formulation of weld pool determination. *Appl. Math. Lett.* **25** (2012) 374–379.

[8] A. Chakib, A. Nachaoui and M. Nachaoui, Approximation and numerical realization of an optimal design welding problem. *Numer. Methods Part. Differ. Equ.* **29** (2013) 1563–1586.

[9] P. Dvalishvili, A. Nachaoui and T. Tadumadze, Effects of the initial moment and several delays perturbations in the variation formulas for a solution of a functional differential equation with the continuous initial condition. *Georgian Math. J.* **27** (2020) 53–66.

[10] J. Eastman, J. Sass, J.M. Gomes, R.W. Dos Santos and E.M. Cherry, Using delay differential equations to induce alternans in a model of cardiac electrophysiology. *J. Theor. Biol.* **404** (2016) 262–272.

[11] J.E. Forde, *Delay differential equation models in mathematical biology*. Ph.D thesis, University of Michigan, ProQuest LLC, Ann Arbor, MI (2005).

[12] R.V. Gamkrelidze and G.L. Karatishvili, Extremal problems in linear topological spaces. *Math. Syst. Theory* **1** (1967) 229–256.

[13] J. Hadamard, Lectures on Cauchy's Problem in Linear Partial Differential Equations. Dover Publications, New York (1953).

[14] J.H. Holland, *Adaptation in Natural and Artificial Systems: An Introductory Analysis With Applications to Biology, Control, and Artificial Intelligence*. MIT Press (1992).

[15] M. Huang, J. Luo, L. Hu, B. Zheng and J. Yu, Assessing the efficiency of wolbachia driven aedes mosquito suppression by delay differential equations. *J. Theor. Biol.* **440** (2018) 1–11.

[16] J. Kepner, Parallel programming with MatlabMPI. Preprint astro-ph/0107406 (2001).

[17] G.L. Kharatishvili and T.A. Tadumadze, Formulas for the Variation of a Solution and Optimal Control Problems for Differential Equations with Retarded Arguments. *J. Math. Sci. (NY)* **140** (2007) 1–175.

[18] S. Lenhart, V. Protopopescu and J. Yong, Solving inverse problems of identification type by optimal control methods. I: AIP Conference Proceedings. Vol. 411. American Institute of Physics (1997) 87–94.

[19] M. Nachaoui, *Étude théorique et approximation numérique d'un problème inverse de transfert de la chaleur*. Ph.D. thesis, Université de Nantes (2011).
[20] M. Nachaoui, Parameter learning for combined first and second order total variation for image reconstruction. *Adv. Math. Models App.* **5** (2020) 53–69.

[21] A. Nachaoui, J. Abouchabaka and N. Rafalia, Parallel solvers for the depletion region identification in metal semiconductor field effect transistors. *Numer. Algorithms* **40** (2005) 187–199.

[22] M. Nachaoui, A. Chakib and A. Nachaoui, An efficient evolutionary algorithm for a shape optimization problem. *Appl. Comput. Math.* **19** (2020) 220–244.

[23] M. Nachaoui, L. Afraites and A. Laghrib, A regularization by denoising super-resolution method based on genetic algorithms. *Signal Process. Image Commun.* **99** (2021) 116505.

[24] A. Nachaoui, M. Nachaoui and B. Gasimov, Parallel numerical computation of an analytical method for solving an inverse problem. *Adv. Math. Models App.* **6** (2021) 162–173.

[25] J.T. Ottesen, Modelling of the baroreflex-feedback mechanism with time-delay. *J. Math. Biol.* **36** (1997) 41–63.

[26] A.A. Samarskii and P.N. Vabishchevich, Numerical Methods for Solving Inverse Problems of Mathematical Physics. Vol. 52 of *Inverse and Ill-posed Problems Series*. Walter de Gruyter GmbH & Co. KG, Berlin (2007).

[27] L.F. Shampine and S. Thompson, Solving DDEs in matlab. *Appl. Numer. Math.* **37** (2001) 441–458.

[28] T. Tadumadze and A. Nachaoui, Variation formulas of solution for a controlled functional-differential equation considering delay perturbation. *TWMS J. Appl. Eng. Math.* **1** (2011) 58–68.

[29] T. Tadumadze, A. Nachaoui and F. Aboud, On one inverse problem for the linear controlled neutral differential equation. In: *Qualitative Theory of Differential Equations QUALITDE-2018* (2018) 183–184.

[30] L. Tavernini, Continuous-Time Modeling and Simulation. Vol. 2. CRC Press (1996).

Subscribe to Open (S2O)
A fair and sustainable open access model

<figure><img src="image_13.png" /></figure>

This journal is currently published in open access under a **Subscribe-to-Open model (S2O)**. S2O is a transformative model that aims to move subscription journals to open access. Open access is the free, immediate, online availability of research articles combined with the rights to use these articles fully in the digital environment. We are thankful to our subscribers and sponsors for making it possible to publish this journal in open access, free of charge for authors.

**Please help to maintain this journal in open access!**

Check that your library subscribes to the journal, or make a personal donation to the S2O programme, by contacting
**subscribers@edpsciences.org**

More information, including a list of sponsors and a financial transparency report, available at: https://www.edpsciences.org/en/maths-s2o-programme
[20] M. Nachaoui, Parameter learning for combined first and second order total variation for image reconstruction. *Adv. Math. Models App.* **5** (2020) 53–69.

[21] A. Nachaoui, J. Abouchabaka and N. Rafalia, Parallel solvers for the depletion region identification in metal semiconductor field effect transistors. *Numer. Algorithms* **40** (2005) 187–199.

[22] M. Nachaoui, A. Chakib and A. Nachaoui, An efficient evolutionary algorithm for a shape optimization problem. *Appl. Comput. Math.* **19** (2020) 220–244.

[23] M. Nachaoui, L. Afraites and A. Laghrib, A regularization by denoising super-resolution method based on genetic algorithms. *Signal Process. Image Commun.* **99** (2021) 116505.

[24] A. Nachaoui, M. Nachaoui and B. Gasimov, Parallel numerical computation of an analytical method for solving an inverse problem. *Adv. Math. Models App.* **6** (2021) 162–173.

[25] J.T. Ottesen, Modelling of the baroreflex-feedback mechanism with time-delay. *J. Math. Biol.* **36** (1997) 41–63.

[26] A.A. Samarskii and P.N. Vabishchevich, Numerical Methods for Solving Inverse Problems of Mathematical Physics. Vol. 52 of *Inverse and Ill-posed Problems Series*. Walter de Gruyter GmbH & Co. KG, Berlin (2007).

[27] L.F. Shampine and S. Thompson, Solving DDEs in matlab. *Appl. Numer. Math.* **37** (2001) 441–458.

[28] T. Tadumadze and A. Nachaoui, Variation formulas of solution for a controlled functional-differential equation considering delay perturbation. *TWMS J. Appl. Eng. Math.* **1** (2011) 58–68.

[29] T. Tadumadze, A. Nachaoui and F. Aboud, On one inverse problem for the linear controlled neutral differential equation. In: *Qualitative Theory of Differential Equations QUALITDE-2018* (2018) 183–184.

[30] L. Tavernini, Continuous-Time Modeling and Simulation. Vol. 2. CRC Press (1996).

Subscribe to Open (S2O)
A fair and sustainable open access model

<figure><img src="image_13.png" /></figure>

This journal is currently published in open access under a **Subscribe-to-Open model (S2O)**. S2O is a transformative model that aims to move subscription journals to open access. Open access is the free, immediate, online availability of research articles combined with the rights to use these articles fully in the digital environment. We are thankful to our subscribers and sponsors for making it possible to publish this journal in open access, free of charge for authors.

**Please help to maintain this journal in open access!**

Check that your library subscribes to the journal, or make a personal donation to the S2O programme, by contacting
**subscribers@edpsciences.org**

More information, including a list of sponsors and a financial transparency report, available at: https://www.edpsciences.org/en/maths-s2o-programme