1. Introduction

Differential dynamic systems with time delays play an important role in the modeling of real-life phenomena in various fields of applications. There exists an extensive literature on delay differential models in biology and biomedicine; see e.g., [6,7,16,18,24,26,27]. Though there is a vast literature on optimal control problems with time delays in control and state variables, so far only a few papers have applied the framework of optimal control with delays to biomedicine; cf. the recent papers [15,21,23]. The aim of this paper is to present two case studies in biomedicine which illustrate the application of delayed optimal control problems and demonstrate that there exist efficient numerical techniques to solve such problems.

In Section 2, we consider optimal control problems with multiple time delays in control and state variables The control process can be subject to mixed control-state constraints. We review the necessary optimality conditions that were derived in [10] in the form of a Pontryagin type Minimum Principle. It is assumed that the so-called commensurability assumption holds which requires that the time delays and the terminal time are integer multiples of a joint stepsize. This assumption also underlies the discretization and nonlinear programming techniques that are briefly reviewed in Section 3. In Section 4, we study the delay differential model for tumour-immune-response with chemo-immunotherapy in Rihan et al. [21]. Aside from the state delay in this model we introduce a time delay in the control variable representing the immune therapy. The delay accounts for the fact that the human immune system needs some time to respond to the immune therapy. In contrast to the $L^2$-type objectives in [21] we consider $L^1$-type objectives characterized by linearly appearing controls which seem to be more appropriate in the biological framework; cf. the remarks in [22]. The computations show that all controls are of bang-bang type. For the non-delayed problem, we can verify the second-order sufficient conditions in [19,17] and thus show that the computed solution provides a strict strong minimum.

Section 5 considers the delay differential model in [7] describing the spread of Hepatitis B virus (HBV). The dynamical model exhibits a delay in the state variables. We introduce a control variable into this model and formulate an optimal control problem using $L^1$-type objectives. Again, all controls are of bang-bang type. We show that the non-delayed bang-bang control provides a strict strong minimum, whereas the delayed controls are extremal solutions that satisfy the necessary conditions with high accuracy.

---

2. Optimal control problems with multiple time-delays in state and control variables

---

2.1. Problem statement

Let $x(t) \in \mathbb{R}^n$ denote the state variable and $u(t) \in \mathbb{R}^m$ the control variable at time $t \in [0,t_f]$ with fixed terminal time $t_f > 0$. The time-delays in the state and control variables are given by a constant vector $(\tau_1,\ldots,\tau_d)\in\mathbb{R}^d$ satisfying

$0 = :\tau_0 < \tau_1< \ldots < \tau_d.$

Thus $\tau_0$ represents the non-delayed variables. In [9,10] we have studied the following optimal control problem with multiple time-delays and mixed control-state constraints (MDOCP): determine a pair of functions $(x,u) \in W^{1,\infty}([0,t_f],\mathbb{R}^n) \times L^\infty([0,t_f],\mathbb{R}^m)$ that minimize the functional in Mayer form

$J(x,u) = g(x(t_f))$

(1)

subject to the delayed (retarded) differential equation, boundary conditions and mixed control-state inequality constraints

$\dot x(t) = f(t,x(t-\tau_0),\ldots,x(t-\tau_d),u(t-\tau_0),\ldots,u(t-\tau_d)), \;\; {\rm a.e.}\ \ t \in [0,t_f],$

(2)

$x(t) = x_0(t), \ \ \ t \in [-\tau_d,0],$

(3)

$u(t) = u_0(t), \ \ \ t \in [-\tau_d,0),$

(4)

$\psi(x(T)) = 0,$

(5)

$C(t,x(t-\tau_0),\ldots,x(t-\tau_d),u(t-\tau_0),\ldots,u(t-\tau_d)) \leq 0 , \; \; {\rm a.e.}\ \ t\in [0,t_f].$

(6)

The functions $g:\mathbb{R}^{n}\rightarrow \mathbb{R}, f:[0,t_f]\times\mathbb{R}^{(d+1)\cdot n}\times\mathbb{R}^{(d+1)\cdot m}\rightarrow \mathbb{R}^n$, $\psi :\mathbb{R}^{n}\rightarrow \mathbb{R}^q (0\leq q\leq n),$ and $C:[0,t_f]\times\mathbb{R}^{(d+1)\cdot n}\times\mathbb{R}^{(d+1)\cdot m}\rightarrow \mathbb{R}^p$ are assumed to be continuously differentiable, while the functions $x_0 : [-\tau_d,0] \rightarrow \mathbb{R}^n, u_0: [-\tau_d,0] \rightarrow \mathbb{R}^m$ only need to be continuous.

Without lack of generality we have assumed that the cost functional is given in Mayer form (1). It is well known that an objective in Bolza form,

$J(x,u) = g(x(t_f)) + \int\limits_0^{t_f} L(t,x(t-\tau_0),\ldots,x(t-\tau_d),u(t-\tau_0),\ldots,u(t-\tau_d))\,{\rm d}t,$

can be reduced to Mayer form by introducing an additional state variable $x_{n+1}$ defined by

$\dot x_{n+1}(t) = L(t,x(t-\tau_0),\ldots,x(t-\tau_d),u(t-\tau_0),\ldots,u(t-\tau_d)), \ \ \ x_{n+1}(0) = 0 .$

Then we have to minimize the functional $\tilde J(x,x_{n+1},u) = g(x(t_f)) + x_{n+1}(t_f)$.

In the following, we shall use the placeholder variables $y_0,y_1,\ldots,y_d$ for the delayed state variables and $v_0,v_1,\ldots,v_d$ for the delayed control variables. The delayed variables are defined by

$\label{delayed-variables} y_{\delta}(t) = x(t-\tau_{\delta}), \ \ \ v_{\delta}(t) = u(t-\tau_{\delta})\ \ (\delta = 0,1,\ldots,d).$

(7)

Note that we do not necessarily assume an equal number of state and control delays. The case of an unequal number of delays in state and control variables is included in this formulation as we admit that

$\frac{\partial h}{\partial y_\delta} = 0 \ \ \mbox{or} \ \ \frac{\partial h}{\partial v_\delta} = 0, \ \ h \in \{f,C,L\}, \ \ \mbox{for some}\ \ \delta\in\{0,\ldots, d\}.$

---

2.2. Minimum principle: First-order necessary conditions

A Pontryagin-type minimum principle for problem (MDOCP) has been derived in [9,10]. The main result requires that all positive time delays $\tau_1,\ldots,r_d$ can be expressed as integer multiples of a sufficiently small positive constant (stepsize).

Assumption 2.1 (Commensurability Condition). Assume that there exist a constant $h>0$ and integers $k_1,\ldots,k_d,N$ with

$\tau_\delta = k_\delta h (\delta = 1,\ldots,d) \ \ \mathit{and} \ \ t_f = N\,h .$

(8)

In view of $0 = \tau_0 < \tau_1< \ldots < \tau_d$ we have $0 < k_1 < \ldots < k_d$. Then in analogy to the non-delayed case we define the Hamiltonian function by

$\label{MDOCP-freeHam} H(t,y_0,\ldots,y_d,v_0,\ldots,v_d,\lambda) = \lambda\, f(t,y_0,\ldots,y_d,v_0,\ldots,v_d), \ \ \ \lambda \in\mathbb{R}^n ,$

(9)

where the adjoint variable $\lambda \in \mathbb{R}^n$ is a row vector. The augmented Hamiltonian function is defined by adjoining the mixed control-state constraint (6) to the Hamiltonian using a multiplier $\mu \in \mathbb{R}^p$ (row vector):

$\label{MDOCP-augHam} \begin{array}{rcl} {\mathcal H}(t,y_0,\ldots,y_d,v_0,\ldots,v_d,\lambda,\mu)& = &H(t,y_0,\ldots,y_d,v_0,\ldots,v_d,\lambda) \\ && + \mu \, C(t,y_0,\ldots,y_d,v_0,\ldots,v_d). \end{array}$

(10)

For ease of notation we refrain from denoting an optimal pair

$(x,u) \in W^{1,\infty}([0,t_f],\mathbb{R}^n) \times L^\infty([0,t_f],\mathbb{R}^m)$

by a hat or a similar symbol. We require the following regularity condition for the active control-state constraints.

Assumption 2.2 (Regularity Condition). Let $(x,u)$ be a locally optimal pair and let

$J_0(t): = \{j\in\{1,\ldots,p\}\;|\; C_j(t, x(t-\tau_0),\ldots, x(t-\tau_d), u(t-\tau_0),\ldots, u(t-\tau_d)) = 0\}$

denote the set of active indices for the inequality constraints (6). Assume that the gradients

$\frac{\partial C_j(t, x(t-\tau_0),\ldots, x(t-\tau_d), u(t-\tau_0),\ldots, u(t-\tau_d))} {\partial (v_0,\ldots, v_d)},\ \ \ j\in J_0(t),$

(11)

are linearly independent.

The following theorem summarizes the first-order necessary conditions for optimality for the control problem (MDOCP) [10].

Theorem 2.3. (Minimum Principle for Optimal Control Problems with Multiple Time-Delays [10]): Let $(x,u)$ be a locally optimal pair for (MDOCP) with delays $0 = \tau_0 < \tau_1 < \ldots \tau_d$ that satisfies the commensurability condition (8) and the regularity condition 2.2. Then there exist an adjoint (costate) function $\lambda\in{W}^{1,\infty}([0,t_f],\mathbb{R}^n)$, a number $\lambda_0 \geq 0$, a multiplier function $\mu \in {L}^\infty ([0,t_f],\mathbb{R}^p)$ and a multiplier $\nu\in\mathbb{R}^q$, such that the following conditions hold for a.e. $t\in[0,t_f]$:

1. Advanced Adjoint Differential Equation:

$\label{MP-adjoint} \dot{\lambda}(t) = -\sum\limits_{\delta = 0}^d\chi_{[0,t_f-\tau_\delta]}(t){\mathcal H}_{y_\delta}(t+\tau_\delta),$

(12)

where ${\mathcal H}_{y_\delta}[t] = {\mathcal H}_{y_\delta}(t,x(t-\tau_0),\ldots, x(t-\tau_d), u(t-\tau_0),\ldots,u(t-\tau_d),\lambda(t),\mu(t))$ and $\chi_{[0,t_f-\tau_\delta]}$ is the characteristic function of the interval $[0,t_f-\tau_\delta]$.

2. Transversality Condition:

$\label{MP-trans} \lambda(t_f) = \lambda_0 g_x( x(t_f))+\nu\, \psi_x( x(t_f)).$

(13)

3. Minimum Condition for the Hamiltonian:

$\label{MP-min} \begin{split} & \sum\limits_{\delta = 0}^d \chi_{[0,t_f-\tau_\delta]}(t) H [t+\tau_\delta]\\ & \leq H(t,\ldots,u, u (t-\tau_1),\ldots, u (t-\tau_d),\lambda(t))\\ & \ \ \ + \sum\limits_{\delta = 1}^{d-1} \chi_{[0,t_f-\tau_\delta]}(t) \,H(t+\tau_\delta,\ldots, u (t+\tau_\delta-\tau_{\delta-1}),u, u(t+\tau_\delta-\tau_{\delta+1}),\ldots)\\ & \ \ \ +\chi_{[0,t_f-\tau_d]}(t) \,H(t+\tau_d,\ldots, u(t+\tau_d-\tau_1),\ldots, u (t+\tau_d-\tau_{d-1}),u,\lambda(t)) \end{split}$

(14)

for all $u\in\mathbb{R}^m$ satisfying

$\begin{split} C(t, x(t-\tau_0),\ldots, x(t-\tau_d), &u(t-\tau_0),\ldots, u (t-\tau_{\delta-1}),u, \\ & u(t-\tau_{\delta+1}),\ldots, u(t-\tau_d)) \; \leq 0 \ \ \mathit{for} \ \ \delta = 0,\ldots,d, \end{split}$

where $H[t] = H(t, x(t-\tau_0),\ldots,x(t-\tau_d), u(t-\tau_0),\ldots, u(t-\tau_d),\lambda(t))$.

4. Local Minimum Condition for the Augmented Hamiltonian Function:

$\label{MP-locmin} \sum\limits_{\delta = 0}^d \chi_{[0,t_f-\tau_\delta]}(t) {\mathcal H}_{v_\delta}[t+\tau_\delta] = 0.$

(15)

5. Non-negativity of Multiplier and Complementarity Condition: for $t\in [0,t_f]$,

$\label{MP-nmcc} \mu(t)\geq 0 , \;\; \mu(t)C(t, x(t-\tau_0),\ldots,x(t-\tau_d,u(t-\tau_0), \ldots , u(t-\tau_d)) = 0.$

(16)

---

3. Numerical discretization methods

Similar to the case of non-delayed differential equations, we can employ integration methods of Runge-Kutta type or multistep methods, e.g., the Euler method and trapezoidal rule, to discretize the delay differential equation

$\dot x(t) = f(t,x(t-\tau_0),\ldots,x(t-\tau_d),u(t-\tau_0),\ldots,u(t-\tau_d)).$

Any integration method based on an equidistant discretization scheme utilizes a uniform step size $h > 0$. Due to the presence of time-delays it is crucial to match the delays $\tau_1,\ldots,\tau_d$ to the grid. This is ensured by the commensurability condition (8) in Assumption 2.1. For this purpose, let $h>0$ be a step size satisfying (8), i.e.

$\tau_\delta = k_\delta h \ \ (\delta = 0,\ldots,d), \ \ t_f = Nh,$

with integers $0 = k_0 < k_1<\ldots <k_d$ and $N$. Note that this grid can be refined by use of any integer fraction of $h$, This defines an equidistant discretization mesh with grid points $t_i = ih$ for $i = 0,1,...,N$.

Let $x_i \in \mathbb{R}^n$ and $u_i \in \mathbb{R}^m$ denote approximations of $x(t_i)$ and $u(t_i)$ at the grid points $t_i$ for $i = 0,1,\ldots,N$. For convenience, we shall use the abbreviations

$f_i = f(t_i,x_i,x_{i-k_1},\ldots,x_{i-k_{\delta}},u_i,u_{i-k_1},\ldots,u_{i-k_{\delta}}).$

The initial value profiles $x_0(\cdot)$ and $u_0(\cdot)$ provide the values

$\label{NLP-initial-x-u} x_{-i} = x_0(-ih) \ \ (i = 0,..,k_d),\ \ u_{-i} = u_0(-ih) \ \ (i = 1,..,k_d).$

(17)

Since the focus in this paper is not on discussing various numerical methods, we present only two integration methods that can be easily implemented. The simplest method is the first order method of Euler which is defined by the recursion

$\label{Euler} x_{i+1} = x_i + h f_i \,, \ \ \ i = 0,1,\ldots,N-1 .$

(18)

The trapezoidal rule is an implicit method of second order:

$\label{Trapezoidal} x_{i+1} = x_i + \frac{1}{2} h \,( f_i + f_{i+1}) \,, \ \ \ i = 0,1,\ldots,N-1 .$

(19)

Then for the Euler method and the optimization variable

$z: = (u_0,x_1,u_1,x_2,...,u_{N-1},x_N) \in \mathbb{R}^{N(m+n)}$

we obtain the following nonlinear programming problem (NLP) with equality and inequality constraints:

$\label{NLP_cost} \mbox{Minimize} \ \ \ J(z) = g(x_N)$

(20)

subject to

$x_{i+1} = x_{i}+hf(t_i,x_{i-k_0},\ldots,x_{i-k_d},u_{i-k_0},\ldots,u_{i-k_d}),\ \ \ i = 0,\ldots,N-1,$

(21)

$C(t_i,x_{i-k_0},\ldots,x_{i-k_d},u_{i-k_0},\ldots,u_{i-k_d})\leq 0,\ \ \ \ \ i = 0,\ldots,N-1,$

(22)

$\psi(x_N) = 0 ,$

(23)

and initial values (17). Using the trapezoidal method (19) we simply replace the equations (21) by the equations defined in (19).

Let $\lambda = (\lambda_0,\lambda_1,\ldots,\lambda_{N-1}) \in \mathbb{R}^{n\cdot N}, \, \lambda_i \in \mathbb{R}^n \;(i = 0,\ldots,N-1)$, be the Lagrange multipliers for equations (21) and let $\mu = (\mu_0,\mu_1,\ldots,\mu_{N-1}) \in \mathbb{R}^{p\cdot N}, \, \mu_i \in \mathbb{R}^p \ (i = 0,\ldots,N-1)$, be the multipliers for the inequality constraints (22) and $\nu_N \in \mathbb{R}^q$ be the multiplier for the boundary condition (23). In [9,10] we have discussed the Karush-Kuhn-Tucker (KKT) necessary optimality conditions for the (NLP) using the Euler scheme (18) and showed that the property of consistency holds. This means that the Lagrange multipliers provide approximations for the adjoint variable $\lambda(t)$, the multiplier $\mu(t)$ and $\nu$ according to

$\label{approximation} \lambda(t_i) \approx \lambda_i \in \mathbb{R}^n ,\ \ \ \mu(t_i) \approx \mu_i/h \in \mathbb{R}^p \ \ (i = 0,...,N-1), \ \ \ \nu_N \approx \nu .$

(24)

This follows from the fact that the Lagrange multipliers $\lambda_i$ satisfy the advanced adjoint equations using the same discretization scheme in a backward mode.

To solve the optimization problem (NLP) in (20)-(22) numerically, we employ the Applied Modeling Programming Language (AMPL) developed by Fourer, Gay and Kernighan [8] which can be linked to the interior-point optimization solver IPOPT developed by Wächter et al. [28] or to the SQP solver WORHP by Büskens and Gerdts [4]. Every solver provides the Lagrange multipliers and therefore gives access to approximations of adjoint variables and multiplier functions for the control problem (MDOCP) according to (24). Thus we can test whether the numerical solution is an extremal solution which satisfies the necessary optimality conditions in Theorem 2.3.

---

4. Optimal control of chemo-immuno-therapy

---

4.1. Optimal control problem

We consider the delay differential model in Rihan et al. [21] that proposes a chemo-immuno-therapy of cancer. The authors introduce a time delay only in the state variable and present a stability analysis of drug free steady states. We shall extend the model by including also a control delay in the control $u_2$ of immune therapy. The delay accounts for the fact that the human immune system takes some time to respond to the immune therapy. The state and variables have the following meaning:

$E$: concentration of effector cells (plasma B cells, producing antibodies).

$T$: concentration of tumour cells.

$N$: concentration of healthy cells.

$U$: concentration of cytostatic agent for chemotherapy.

$u_1$: dose control for chemotherapy,

$u_2$: dose control for immune therapy of the effector cells.

Denoting the state delay by $\tau_1$ and the control delay by $\tau_2$, the dynamical system is given by

$\label{chemo-immune-system} \begin{array}{lcl} \dot E(t)& = &\sigma+\left(\tfrac{\rho}{\eta+ T(t-\tau_1)}-\mu_e\right)E(t-\tau_1)T(t-\tau_1) \\ &&-(\delta+a_1(1-e^{-U(t)}))E(t) + u_2(t-\tau_2)s_1,\\ \dot T(t)& = &\left(r_2(1-\beta T(t))-n_TE(t)-c_1N(t)-a_2(1-e^{-U(t)})\right)T(t),\\ \dot N(t)& = &\left(r_3(1-\beta_2 N(t))-c_2T(t)-a_3(1-e^{-U(t)})\right)N(t),\\ \dot U(t)& = &u_1(t)-d_1U(t). \end{array}$

(25)

The initial values and initial functions for the delayed state and control variables are as follows:

$\label{initial-values} \begin{array}{rclrcll} E(0)& = &E_0 = 0.3,&E(t)& = &E_0 &\forall \; -\tau_1\leq t\leq 0, \\ T(0)& = &T_0 = 300,&T(t)& = &T_0 &\forall \; -\tau_1\leq t\leq 0, \\ N(0)& = &N_0 = 0.9,&u_2(t)& = &0 &\forall \; -\tau_2\leq t < 0.\\ U(0)& = &U_0 = 0.0. && && \end{array}$

(26)

We shall consider the control constraints

$\label{control-constraint} 0 \leq u_k(t)\leq u_{k,\max} \ \ \ \forall \; t\in [0,t_f] \ \ (k = 1,2).$

(27)

Let us denote the state and control variables by

$x = (E,T,N,U) \in \mathbb{R}^4, \ \ \ u = (u_1,u_2) \in \mathbb{R}^2.$

For notational convenience, we simplify the notations (7) for the delayed state and control variables. In the context of the dynamical system (25) it is more convenient to consider the delayed state variables $y_1, y_2$ and control variable $v_2$ defined by

$\label{delayed-variables-Rihan} \begin{array}{l} y_1(t) = x_1(t-\tau_1) = E(t-\tau_1), \;\, y_2(t) = x_2(t-\tau_1) = T(t-\tau_1) , \\ v_2(t) = u_2(t-\tau_2). \end{array}$

(28)

With these notations the dynamical system (25) can be written as

$\label{dynamics} \dot x(t) = f(x(t),y_1(t),y_2(t),u(t),v_2(t)).$

(29)

Then the optimal control problem is as follows: determine a control function $u = (u_1, u_2) \in L^\infty([0,t_f],\mathbb{R}^2)$ that minimizes the objective functional

$\label{objective} J_p(x,u) = \int_0^{t_f}\,(T(t)-E(t)+ B_1(u_1(t))^p+B_2(u_2(t))^p)\,{\rm d}t \ \ \ \ (p = 1,2)$

(30)

subject to the dynamic constraints (25), initial conditions (26) and control constraints (27). The objective functional (30) represents a trade-off between minimizing the tumour cells and the total doses of the cytotoxic and immunologic agents on one hand and maximizing the plasma cells on the other hand. The constants $B_1 > 0, B_2 >0$ are appropriate weights which are listed in Table 1 together with the system parameters.

Rihan et al. [21] consider only the $L^2$-type functional $J_2(x,u)$ in (30) which is quadratic in the control variable $u$. $L^2$-type functionals are often used in economics to describe, e.g., production costs, but are mostly not appropriate in a biological framework; cf. the remarks in [22]. The $L^1$ functional $J_1(x,u)$ incorporates the total amount of drugs used as a penalty and thus appears to be more realistic. For that reason, we shall mainly focus on the functional $J_1(x,u)$ in the sequel.

Now we apply the necessary optimality conditions in the form of a Minimum Principle as stated in Theorem 2.3. Denoting the adjoint variable by the row vector $\lambda = (\lambda_E,\lambda_T,\lambda_N,\lambda_U) \in \mathbb{R}^4$, the Hamiltonian for the objective $J_1(x,u)$ and the control system (29) is given by

$\label{Hamiltonian-Rihan} H(x,y_1,y_2,u,v_2,\lambda) = T-E+ B_1 u_1 + B_2 u_2 + \lambda f(x,y_1,y_2,u,v_2) .$

(31)

According to Theorem 2.3 (1), the advanced adjoint equations are given by

$\label{adjoint-equation-Rihan} \begin{array}{l} \dot\lambda_E(t) = - H_E[t] - \chi_{\,[0,t_f-\tau_1]\,}(t)\, H_{y_1}[t+\tau_1] , \\ \dot\lambda_T(t) = - H_T[t] - \chi_{\,[0,t_f-\tau_1]\,}(t)\, H_{y_2}[t+\tau_1], \\ \dot\lambda_N(t) = - H_N[t] , \dot\lambda_U(t) = - H_U[t] . \end{array}$

(32)

We do not write out the adjoint variables explicitly, since the adjoint variables can be computed as Lagrange multipliers of the discretized control problem as explained in the preceding section. Due to the free terminal state, the transversality condition (13) is

$\label{transversality} \lambda(t_f) = (0,0,0,0).$

(33)

The optimal control $u(t)$ minimizes the sum of Hamiltonians in (14). Since both controls appear linearly in the Hamiltonian, the minimizing controls are determined by the switching functions

$\label{switching-function} \begin{array}{l} \phi_1(t) = H_{u_1}[t] = B_1+\lambda_U(t),\\ \phi_2(t) = H_{u_2}(t)+\chi_{[0,t_f-\tau_2]}(t)H_{v_2}[t+\tau_2] = B_2 +\chi_{[0,t_f-\tau_2]}(t)\lambda_E(t+\tau_2)s_1, \end{array}$

(34)

according to the control law

$\label{control-law} u_k(t) = \left\{ \begin{array}{lcl} 0,&\mbox{if}&\phi_k(t) > 0 \\ u_{k,\rm max},&\mbox{if}&\phi_k(t) < 0 \\ {\rm singular},&\mbox{if}&\phi_k(t) = 0 \;\; \forall \; t \in I_s \subset [0,t_f] \\ \end{array} \right\}, \ \ \ k = 1,2.$

(35)

Singular controls will not be discussed further, since our computations only yield bang-bang controls. Due to the transversality condition $\lambda(t_f) = 0$ the switching functions satisfy $\phi_k(t_f) = B_k > 0$ for $k = 1,2$. Hence, the control law (35) shows that $u_k(t) = 0$ holds on a terminal interval $[t_k, t_f]$ for $k = 1,2$. Parameters for the subsequent computations are given in the Table 1.

---

4.2. Optimal solution of the non-delayed control problem

First, we present the solution for the non-delayed control problem with $\tau_1 = \tau_2 = 0$ and the functional $J_1(x,u)$. Recall the upper control bounds $u_{1,\rm max} = u_{2,\rm max} = 1$, the terminal time $t_f = 30$ (days) and the weights $B_1 = 5$ and $B_2 = 10$ from Table 1. Applying AMPL/IPOPT with $N = 3000$ grid points and the trapezoidal rule (19) we find the following bang-bang controls $u_k(t)$ with only one switch at $t_k$,

$\label{control-structure-rihan} u_k(t) = \left\{ \begin{array}{lcl} 1&\mbox{for}&0 \leq t < t_k \\ 0&\mbox{for}&t_k \leq t \leq t_f \end{array} \right\} \ \ \ (k = 1,2), \ \ \ 0 < t_1 < t_2 < t_f \,.$

(36)

To obtain a refinement of the solution, we solve the Induced Optimization Problem (IOP) with the switching times $t_1$ and $t_2$ as optimization variables; cf. [17,19]). The arc-parametrization method [17] and the optimal control package NUDOCCCS due to Büskens [2] yield the following numerical results

$\begin{array}{rclrclrcl} J_1(x,u) & = & 1399.02 ,& t_1 & = & 3.93 031 ,&t_2 & = & 9.76 562, \\ E(t_f) & = & 0.64 03 03 ,&T(t_f) & = & 0.18 07 26, &N(t_f) & = & 0.90 49 68 , \\ U(t_f)& = & 2.9 69 62. &&&& \end{array}$

The initial values of the adjoint variables are

$\lambda_E(0) = -770.13 , \; \lambda_T(0) = 2.9980 , \; \lambda_N(0) = -0.027548, \; \lambda_U(0) = -281.11.$

The non-delayed solution is shown in Figure 1. A common strategy in medical practise is the administration of a pulse therapy or a blockwise application of drugs. Such a strategy is promoted by the controls in Figure 1.

Now we show that the second-order sufficient conditions in [19], Chapter 7, are satisfied for the bang-bang control (36). For that purpose, we have to check two further conditions. First, notice that the objective $J_1(x,u)$ becomes a function $J_1(t_1,t_2)$ of the two switching times $t_1,t_2$, if we assume the control structure (36). The Hessian of $J_1(t_1,t_2)$ is computed as the positive definite $2\times 2$ matrix

$D^2 J_1 (t_1,t_2) = \left( \begin{array}{cc} 19.167 & 11.120\\ 11.120 & 10.887 \end{array} \right).$

Furthermore, as can be seen in Figure 2, the following strict bang-bang property with respect to the Minimum Principle holds for $k = 1,2$:

$\label{strict-bang-bang} \phi_k(t) < 0 \ \ \ \forall \; 0 \leq t < t_k, \ \ \ \dot\phi_k(t_k) > 0,\ \ \ \phi_k(t) > 0 \ \ \ \forall \; t_k < t \leq t_f \, .$

(37)

Hence, the solution shown in Figure 1 provides a strict strong minimum.

We briefly compare the solutions for the functionals $J_1(x,u)$ and $J_2(x,u)$. The controls $u_1$ and $u_2$ for the functional $J_2(x,u)$ are continuous, since the strict Legendre-Clebsch condition holds and the Hamiltonian has a unique minimum with respect to $u_1$ and $u_2$. Figure 3 displays a comparison of the controls $u_1$ and $u_2$ for both functionals. The state variables for the functional $J_2(x,u)$ are very similar to those shown in Figure 1 and thus are not displayed here. The functional value is $J_2(x,u) = 1392.88$ versus $J_1(x,u) = 1399.02$ and the final state is computed as

$E(t_f) = 0.615 728, \; T(t_f) = 0.108 124 , \; N(t_f) = 0.903 899, \; U(t_f) = 3.20 922 .$

---

4.3. Numerical solution of the delayed control problem

We choose the state delay $\tau_1 = 1.5$ and the control delay $\tau_2 = 3$. To obtain a rather precise reference solution, we apply AMPL/IPOPT with $N = 6000$ grid points and tolerance $tol = 10^{-8}$. As in the non-delayed case we obtain a bang-bang control $u(t) = (u_1(t),u_2(t))$, where each $u_k(t)$ has only one switch at $t_k$:

$\label{control-structure-rihan-delayed} u_k(t) = \left\{ \begin{array}{lcl} 1&\mbox{for}&0 \leq t < t_k \\ 0&\mbox{for}&t_k \leq t \leq t_f \end{array} \right\} \ \ (k = 1,2), \ \ 0 < t_1 < t_2 < t_f \,.$

(38)

We obtain the numerical results

$\begin{array}{rclrclrcl} J_1(x,u)& = & 2126.69 ,&t_1 & = & 4.692 ,&t_2 & = & 10.42 , \\ E(t_f) & = & 0.6612 58 ,&T(t_f) & = & 0.1362 62 , &N(t_f) & = & 0.9027 47 , \\ U(t_f) & = & 3.555 46 . &&& \end{array}$

The initial values of the adjoint variables are

$\lambda_E(0) = -485.41 , \; \lambda_T(0) = 2.2403, \; \lambda_N(0) = -0.022090,\; \lambda_U(0) = -248.50 .$

Using the Euler method (18) with the same number $N = 6000$ grid points, the numerical results are less accurate by two decimals. The control and state trajectories are shown in Figure 4. Figure 5 displays the controls and the switching functions in a neighborhood of the switching times. The zoom into the controls confirms that the control law (35) is precisely satisfied and that the strict bang-bang property (37) holds as well for the delayed solution. Unfortunately, we can not check any kind of sufficient conditions for the delayed solution, since numerically verifiable sufficient conditions are not available in the literature.

Finally, as in the non-delayed case we briefly compare the solutions for the functionals $J_1(x,u)$ and $J_2(x,u)$. The controls $u_1$ and $u_2$ for the functional $J_2(x,u)$ are continuous, since the strict Legendre-Clebsch condition holds and the Hamiltonian has a unique minimum with respect to $u_1$ and $u_2$. Figure 6 displays a comparison of the controls $u_1$ and $u_2$ for both functionals.

---

4.4. Numerical solution of the delayed control problem with mixed control-state constraint $U(t)+u_2(t) \leq 3$

We add the following mixed control-state constraint to the delayed optimal control problem:

$\label{mixed-constraint} U(t) + u_2(t) \leq 3 \ \ \ \forall \; t \in [0,t_f] .$

(39)

This constraint means that sum of the cytotoxic agent and the immune dose is bounded from above. Here we consider the augmented Hamiltonian

$\label{augmented-Hamiltonian} {\mathcal H}(x,y_1,y_2,u,v_2,\lambda,\mu) = H(x,y_1,y_2,u,v_2,\lambda) + \mu (U+u_2),$

(40)

where the mixed constraint is adjoined to the Hamiltonian (31) by a multiplier $\mu \geq 0$. The local minimum condition (15) yields

$\label{local-minimum-condition} 0 = \mathcal{H}_{u_2}[t]+\chi_{[0,t_f-\tau_2]}(t)\, \mathcal{H} _{v_2}[t+\tau_2] = \phi_2(t) + \mu(t),$

(41)

where $\phi_2(t) = B_2 + \chi_{ [0,tf-\tau_2]}(t)\,\lambda_E(t+\tau_2) s_1$ is the switching function defined in (34). The multiplier satisfies the complementarity condition $\mu(t)(U(t)+u_2(t)-3) = 0$ for $t \in [0,t_f]$. Hence, on a boundary arc with $U(t)+u_2(t) = 3$ for $t \in [t_1,t_2]$ we obtain an explicit formula of the multiplier in view of (41):

$\label{multiplier} \mu(t) = - \phi_2(t) = - B_2 - \chi_{[0,t_f-\tau_2]}(t)\lambda_E(t+\tau_2)s_1 \ \ \ \forall \; t \in [t_1,t_2].$

(42)

Computations show that the control $u_2(t)$ is constant on a boundary arc and thus we obtain by differentiation

$0 = \dot U(t) = u_1(t) - d_1 U(t) = u_1(t) - d_1\,(3-u_2(t)).$

Since we have $u_2(t) = 1$ on a boundary arc, the control $u_1(t)$ on the boundary arc is given by

$u_1(t) = d_1\,(3-u_2(t)) = 0.02 \ \ \ (d_1 = 0.01).$

Using the trapezoidal method (19) with $N = 3000$ grid points we find the control structure

$\label{control-structure-rihan-delayed-mixed} u_1(t) = \left\{ \begin{array}{lcl} 1&\mbox{for}&0 \leq t < t_1 \\ 0.02&\mbox{for}&t_1 \leq t < t_2 \\ 0&\mbox{for}&t_2 \leq t \leq t_f \end{array} \right\},\ \ u_2(t) = \left\{ \begin{array}{lcl} 1&\mbox{for}&0 \leq t < t_3 \\ 0&\mbox{for}&t_3 \leq t \leq t_f \end{array} \right\}$

(43)

with $0 < t_1 < t_2 < t_3 < t_f$ and the boundary arc $[t_1,t_2]$. We obtain the numerical results:

$\begin{array}{rclrclrcl} J_1(x,u)& = & 2236.06 ,&t_1 & = & 2.045 ,&t_2 & = & 9.95 , \\ t_3 & = & 10.98,&E(t_f) & = & 0.725 265 ,&T(t_f) & = & 0.100 546 , \\ N(t_f) & = & 0.919 108 ,&U(t_f) & = & 1.63 720 . \end{array}$

---

5. Optimal control of a delay model of Hepatitis B virus infection

---

5.1. Optimal control model

Eikenberry et al. [7] report that currently about two billion people - roughly $30 \%$ of the human population - have been infected by Hepatitis B virus (HBV). The disease has attracted considerable attention from mathematical biologists who have developed various models to study the HBV dynamics. Eikenberry et al. [7] present a dynamical model with state variables

$x: \ \ {\rm number~of~healthy~cells,}\\ p: \ \ {\rm number~of ~exposed~cells,}\\ y: \ \ {\rm number~of~infected~cells,} \\ v: \ \ {\rm free~virion~load. }$

The model (4.1)-(4.4) in [7] does not yet involve a control variable. We choose the control variable $u$ as the effect of treatment which corresponds to the coefficient $\gamma$ in the dynamic equation (4.4) in [7]. Denoting the time by $t \in [0,t_f]$ with fixed final time $t_f > 0$ and the delay in the state variable by $\tau \geq 0$, the dynamic system (4.1)-(4.4) in [7] reads as follows:

$\label{dynamics-Eikenberry} \begin{array}{rl} \dot x(t) & = r \,x(t) \left(1-\frac{T(t)}{K}\right) -d\,x(t) - \beta\, v(t) \frac{x(t)}{T(t)}\; , \\ \dot p(t) & = -d\,p(t) + \beta\, v(t) \frac{x(t)}{T(t)} - \beta e^{-d\,\tau} v(t-\tau)\,\frac{x(t-\tau)}{T(t-\tau)} , \\ \dot y(t)& = \beta e^{-d\,\tau} v(t-\tau)\,\frac{x(t-\tau)}{T(t-\tau)} - a\,y(t) , \\ \dot v(t) & = k(1-u(t))\,y(t) - \mu \,y(t) . \end{array}$

(44)

The variable $T$ denotes the total number of cells defined by

$T = x + p + y \, .$

The delay $\tau$ appears in all three variables $x, p, y$. Hence, the initial conditions are given by initial functions for $x, p, y$ and an initial value for $v$:

$\label{initial-condition-Eikenberry} x(t) = x_0 , \; p(t) = p_0 , \; y(t) = y_0 \ \ \mbox{for}\ \ -\tau \leq t \leq 0, \ \ v(0) = v_0 \,.$

(45)

We impose the control constraint

$\label{control-constraint-Eikenberry} 0 \leq u(t) \leq 1 \ \ \ \forall \;\; t \in [0,t_f].$

(46)

Denoting the state vector by $X: = (x,p,y,v) \in \mathbb{R}^4$ and the delayed variable by $Y$, where $Y(t) = X(t-\tau)$, the dynamical system can be written as

$\label{ODE-X} \dot X = f(X,Y, u)$

(47)

with initial functions and conditions given in (45).

The optimal control problem then consists in determining a control function $u \in L^1([0,t_f],\mathbb{R})$ that minimizes the cost functional

$\label{cost-functional} J(X,u) = \int_0^{t_f} (-x(t) + B \, u(t) ) \,{\rm d}t \ \ \ ( B > 0 ),$

(48)

subject to the dynamics (44) with initial conditions (45) and the control constraint (46). The objective functional represents a trade-off between maximizing the number of healthy cells and minimizing the treatment cost.

---

5.2. Necessary optimality conditions: Minimum principle

We briefly discuss the necessary optimality conditions in Theorem 2.1. The Hamiltonian is given by

$\label{Hamiltonian-Eikenberry} H(X,Y,u,\lambda) = - x + Bu + \lambda \, f(X,Y,u),\ \ \ \ \lambda = (\lambda_x, \lambda_p, \lambda_y, \lambda_v) \in \mathbb{R}^4.$

(49)

We do not explicitly write out the advanced adjoint equation (12):

$\label{adjoint-eq-Eikenberry} \dot\lambda(t) = - H_X[t] - \chi_{\,[0,t_f-\tau]\,}(t)\, H_Y[t+\tau] .$

(50)

The control variable $u$ appears linearly in the Hamiltonian and does not involve a delay. Hence, defining the switching functions by

$\label{switching-function-Eikenberry} \phi(t) = H_{u}[t] = B -\lambda_v(t)\,k\,y(t) ,$

(51)

the minimizing control is characterized by the control law

$\label{control-law-Eikenberry} u(t) = \left\{ \begin{array}{lcl} 0,&\mbox{if}&\phi(t) > 0 \\ 1,&\mbox{if}&\phi(t) < 0 \\ {\rm singular},&\mbox{if}&\phi(t) = 0 \;\; \forall \; t \in I_s \subset [0,t_f] \\ \end{array} \right\}.$

(52)

Singular controls will not be discussed further, because we only found bang-bang controls. The following parameters from [7], page 294 below, will be used in our computations:

$\begin{array}{llll} a = 0.011,&d = 0.0039, &\beta = 4.8\cdot 10^{-5},&k = 200, \\ K = 2,&r = 1 ,&\mu = 0.693. \end{array}$

(53)

The state variable $X = (x,p,y,v)$ is scaled by $10^{-11}$ so that we can choose, e.g., the following initial conditions:

$\label{initial-condition-Eikenberry-2} x(t) = 1.4 , \ \ p(t) = 0.3, \ \ y(t) = 0.2 \ \ \ \ \forall \; -\tau \leq t \leq 0, \ \ \ \ v(0) = 500.$

(54)

The time horizon is $t_f = 500$ (days) and the weight parameter in the objective (48) is taken as $B = 0.05 \,.$

---

5.3. Comparison of solutions for several delays

We compare the solutions for the delays $\tau = 0$ (non-delayed solution), $\tau = 10$ and $\tau = 15$. Applying AMPL/IPOPT with $N = 5000$ grid points and using the trapezoidal rule (19), we find a bang-bang control $u(t)$ with only one switch at $t_1$,

$\label{control-structure-Eikenberry} u(t) = \left\{ \begin{array}{lcl} 1&\mbox{for}&0 \leq t < t_1 \\ 0&\mbox{for}&t_1 \leq t \leq t_f \end{array} \right\} .$

(55)

In the non-delayed case, a refinement of the solution is obtained by solving the Induced Optimization Problem (IOP) with respect to the switching time [17,19]. We get the numerical results:

$\begin{array}{lcrclrcl} \tau = 0 &:& J(X,u) = 893.072,&t_1 = 261.70 , \\ \tau = 10 &:& J(X,u) = 913.388,&t_1 = 293.50 , \\ \tau = 15 &:& J(X,u) = 923.032,&t_1 = 304.10 . \end{array}$

A comparison of the controls and switching functions for the delays $\tau = 0, 10, 15$ is shown in Figure 8. The bang-bang control for $\tau = 0$ provides a strict strong minimum, since second-order sufficient conditions (SSC) in [17,19] are satisfied. The numerical test of SSC proceeds as follows. Since the bang-bang control (55) has only one switch at $t_1$, the objective functional becomes a function $J = J(t_1)$ of the scalar optimization variable $t_1$. One verifies numerically that the second derivative is positive: ${\rm d}^2J/{\rm d}t_1^2 = 0.005028 > 0$. Moreover, the following strict bang-bang property [17,19] for the switching function $\phi(t)$ holds; cf. Figure 8, left:

$\phi(t) < 0 \ \ \mbox{for} \;\; 0 \leq t < t_1 \,, \ \ \dot\phi(t_1) > 0 \,, \ \ \phi(t) > 0 \ \ \mbox{for} \;\; t_1 < t \leq t_f = 500 \,.$

Note that the strict bang-bang property is also satisfied for the delayed control with delays $\tau = 10$ and $\tau = 15$. However, as in the preceding section we can not conclude that the delayed controls in Figure 8 provide a strict strong minimum. Figure 9 displays a comparison of the state variables for delays $\tau = 0, \tau = 10, \tau = 15$.

---

6. Conclusion

We presented two applications of delayed optimal control problems in biomedicine. In the first case study, we extended the delay differential model of tumour-immune-response in Rihan et al. [21] by including a time delay in the control variable $u_2$ which represents the immune therapy. The delay is due to the delayed response of the human immune system to the immune therapy. Rihan et al. [21] considered a $L^2$-type objective which is quadratic in the control variables. From a numerical point of view, the control solution in [21] remained a bit obscure.

Therefore, we improved the results in this paper in two regards. First, we considered a more realistic $L^1$-type objective which is linear in the two control variables. Secondly, we applied the discretization and nonlinear programming methods [10] (see Section 3) to obtain extremal solutions that satisfy the necessary optimality conditions in Theorem 2.1 with high accuracy. The computations showed that both controls $u_1$ and $u_2$ are of bang-bang type with only one switch from the upper bound $u_k(t) = u_{k,max}$ to the zero control $u_k(t) = 0$ for $k = 1,2$. Apparently, it is much easier to administer the therapy protocol induced by a bang-bang control then applying a treatment plan resulting from a $L^2$-type objective; cf. Figure 3. In the non-delayed case we could show that the bang-bang controls are indeed optimal, since they satisfy the second-order sufficient conditions in [19,17]. To our knowledge, sufficient conditions for delayed bang-bang controls are not available in the literature. We have also studied the solution under the mixed control-state constraint (39) which combines the cytostatic agent $U(t)$ and the immune control $u_2(t)$. The computations gave very accurate extremal solutions.

The second delay differential model, which describes the spread of Hepatitis B virus, was taken from Eikenberry at al. [7]. We introduced a control variable into the originally uncontrolled model and considered $L^1$-type objectives. For different delays we obtained only bang-bang controls as in the first case study. Sufficient optimality conditions [19,17] could only be verified for the non-delayed bang-bang control.

---