ABSTRACT

The purpose of the present paper is to apply the Pontryagin Minimum Principle to mathematical models of cancer growth. In [1] , I presented a discrete affine model T of cancer growth in the variables C for cancer, GF for growth factors and GI for growth inhibitors. One can sometimes find an affine vector field X on ${ℝ}^3$ whose time one map is T. It is to this vector field we apply the Pontryagin Minimum Principle. We also apply the Discrete Pontryagin Minimum Principle to the model T. So we prove that maximal chemo therapy can be optimal and also that it might not depending on the spectral properties of the matrix A, (see below). In section five we determine an optimal strategy for chemo or immune therapy.

Keywords:

Cancer, Models of Cancer Growth, Pontryagin Minimum Principle

1. Introduction

Our reference to optimal control theory of ODEs is [2] and to control theory of discrete systems [3]. There is a review of papers in optimal control theory applied to cancer [4]. The model we consider here is from [1] and is defined by

$A=\left(\begin{array}{ccc}1+\gamma & \alpha & \beta \\ \delta & 1+{\mu }_F& 0\\ \sigma & 0& 1+{\mu }_I\end{array}\right),T\left(y\right)=Ay+g$ (1)

where $y={\left(C\mathrm{,}GF\mathrm{,}GI\right)}^{\text{T}}$,$g={\left(g_C\mathrm{,}{g}_F\mathrm{,}{g}_I\right)}^{\text{T}}\in {ℝ}^3$, and T denotes transpose. Here $\alpha \in {ℝ}_{+}$,$\beta \in {ℝ}_{-}$,$\delta ,\sigma \in ℝ$,${\mu }_F,{\mu }_I<0$. We assume that $\alpha \delta +\beta \sigma \ne 0$.

The matrix A has characteristic polynomial $-p\left(\lambda \right)$ where

$p\left(\lambda \right)={\lambda }^3-{\lambda }^2\left(3+\gamma +{\mu }_F+{\mu }_I\right)-\lambda (\alpha \delta +\beta \sigma -\left(1+{\mu }_F\right)\left(1+{\mu }_I\right)$ (2)

$-\left(1+\gamma \right)\left(2+{\mu }_F+{\mu }_I\right)+\alpha \delta \left(1+{\mu }_I\right)+\beta \sigma \left(1+{\mu }_F\right)-\left(1+\gamma \right)\left(1+{\mu }_F\right)\left(1+{\mu }_I\right)$ (3)

and this polynomial can have (i) three real roots or (ii) one real root and two imaginary roots. It turns out that asking $\mu ={\mu }_F={\mu }_I$ simplifies matters considerably. Then

$p\left(\lambda \right)=-\left(1+\mu -\lambda \right)\left({\lambda }^2-\left(2+\gamma +\mu \right)\lambda +\left(1+\gamma \right)\left(1+\mu \right)-\alpha \delta -\beta \sigma \right)$ (4)

So the eigenvalues are $1+\mu$ and

${\lambda }_{\pm }=1+\frac{\gamma +\mu }{2}\pm \frac{1}{2}\sqrt{{\left(\gamma -\mu \right)}^2+4\left(\alpha \delta +\beta \sigma \right)}$ (5)

We assume that $1+\mu >0,\mu <0$. In case (i) if the eigenvalues of A, ${\lambda }_{+}\mathrm{,}{\lambda }_{-}\mathrm{,}\tilde{\lambda }$ are positive and distinct and we assume this, then you can find an affine vector field X on ${ℝ}^3$ such that the time one map is d

${\Phi }_1^X=T$ (6)

see [1] and below. In case (ii) if the eigenvalues of A are $1+\mu ,a\pm ib,a>0,1+\mu >0,b\ne 0$ and we assume this, you can also find an affine vector field X on ${ℝ}^3$ such that

${\Phi }_1^X=T$ (7)

In case (i) define a matrix of eigenvectors

$D=\left(\begin{array}{ccc}1+\mu -{\lambda }_{+}& 1+\mu -{\lambda }_{-}& 0\\ -\delta & -\delta & -\beta \\ -\sigma & -\sigma & \alpha \end{array}\right)$ (8)

and in case (ii) define

$U=\left(\begin{array}{ccc}1+\mu -a& -b& 0\\ -\delta & 0& -\beta \\ -\sigma & 0& \alpha \end{array}\right)$ (9)

Then in case (i)

$\Lambda =D^{-1}AD=\left(\begin{array}{ccc}{\lambda }_{+}& 0& 0\\ 0& {\lambda }_{-}& 0\\ 0& 0& \tilde{\lambda }\end{array}\right)$ (10)

if the eigenvalues are distinct and positive and in case (ii)

$\Lambda =U^{-1}AU=\left(\begin{array}{ccc}a& b& 0\\ -b& a& 0\\ 0& 0& 1+\mu \end{array}\right)$ (11)

see [1]. Here $D^{-1}$ denotes the inverse to D. To find an X in case (i) define when the eigenvalues are real, positive and distinct, an affine vector field

$Y\left(x\right)=\left(\begin{array}{ccc}\ln {\lambda }_{+}& 0& 0\\ 0& \ln {\lambda }_{-}& 0\\ 0& 0& \ln \tilde{\lambda }\end{array}\right)x+d$ (12)

Then with the right choice of d

${\Phi }_1^Y\left(x\right)=D^{-1}ADx+D^{-1}g$ (13)

hence if we let

$X\left(y\right)=DY\left(D^{-1}\left(y\right)\right)$ (14)

we get

${\Phi }_1^X=T$ (15)

And to find an X in case (ii) when $a>0,b\ne 0,1+\mu >0$, define

$Y\left(x\right)=\left(\begin{array}{ccc}{a}_1& b_1& 0\\ -b_1& a_1& 0\\ 0& 0& \ln \left(1+\mu \right)\end{array}\right)+\tilde{d}$ (16)

Then with the right choice of $a_1\mathrm{,}{b}_1\in ℝ\mathrm{,}\tilde{d}\in {ℝ}^3$ we get

${\Phi }_1^Y\left(x\right)=U^{-1}AU\left(x\right)+U^{-1}g$ (17)

hence with

$X\left(y\right)=UY\left(U^{-1}\left(y\right)\right)$ (18)

we have that the time one map is

${\Phi }_1^X\left(y\right)=T\left(y\right)$ (19)

See [1] for details of the above and also below.

In section two we solve the problem:

minmize $C\left(T\right)$,$T>0$ (fixed) subject to

${\left(\begin{array}{c}C\\ GF\\ GI\end{array}\right)}^{\prime }=X\left(\begin{array}{c}C\\ GF\\ GI\end{array}\right)+e_{3}u\left(t\right)$ (20)

$u\left(t\right)\in \left[0,g_I^0\right],g_I^0>0$, by first solving it for Y and then infer the solution for X. Here

$e_3=B={\left(0,0,1\right)}^{\text{T}}$ (21)

and $u\left(t\right)$ is piecewise continuous. In section three we apply the discrete Pontryagin minimum principle to the difference equation

$x_{k+1}=A{x}_k+B{u}_k+g$

where

$x_k={\left(C\left(k\right)\mathrm{,}GF\left(k\right)\mathrm{,}GI\left(k\right)\right)}^{\text{T}}$ (22)

$u_k\in \left[0,g_I^0\right],g_I^0>0$. $k=0,1,\cdots ,N-1$, where $N\in \mathbb{N},N\ge 2,x_0=x$, with the objective to minimize $C\left(N\right)$. There are again two cases to consider (i) and (ii) above. If $\mu ={\mu }_F={\mu }_I$ and the eigenvalues are positive and distinct, maximal chemo therapy is optimal. But in case (ii) it is not always optimal.

When ${\mu }_F\ne {\mu }_I$ and the eigenvalues are real and distinct, we produce a counter example to maximal chemo therapy being optimal, see section four. Some solid tumors grow like Gompertz functions, see [5]. There are several important monographs in mathematics and medicine, see [6] - [11]. [12] - [19] are my latest papers on cancer and mathematics. [20] is our reference to roots of cubic polynomials. [21] proves continuous dependence of roots of a polynomial on the coefficients of the polynomial.

In section five we consider optimality of the discrete model T when ${\mu }_F\ne {\mu }_I$. Here we also determine optimal control of the map T.

2. Optimal Control of X

The purpose of this section is to minimize $C\left(T\right)$, subject to

${\left(\begin{array}{c}C\\ GF\\ GI\end{array}\right)}^{\prime }=X\left(\begin{array}{c}C\\ GF\\ GI\end{array}\right)+e_{3}u\left(t\right)$ (23)

$T>0$ fixed and with $\mu ={\mu }_F={\mu }_I$. Let us consider (ii) first. We assume that there are a real eigenvalue $1+\mu$ and two imaginary eigenvalues $a\pm ib,a>0,1+\mu >0,b\ne 0$.

Now define the two by two matrix

$L_1=\left(\begin{array}{cc}{a}_1& b_1\\ -b_1& a_1\end{array}\right)$ (24)

where

$b_1={\tan }^{-1}\left(\frac{b}{a}\right)$ (25)

$a_1=\ln \sqrt{a^2+b^2}$ (26)

This will imply, that

$\exp \left(\begin{array}{cc}{a}_1& b_1\\ -b_1& a_1\end{array}\right)=\left(\begin{array}{cc}a& b\\ -b& a\end{array}\right)=L$ (27)

Also let

$\tilde{B}=\left(\begin{array}{ccc}{a}_1& b_1& 0\\ -b_1& a_1& 0\\ 0& 0& \ln \left(1+\mu \right)\end{array}\right)$ (28)

and define the vector field

$Y\left(x,v\right)=\tilde{B}x+d+U^{-1}{e}_{3}v=\tilde{Y}\left(x\right)+U^{-1}{e}_{3}v$ (29)

which is affine when $v=0$, where $e_1\mathrm{,}{e}_2\mathrm{,}{e}_3$ is the canonical basis in ${ℝ}^3$. Also $x,d\in {ℝ}^3,v\in \left[0,g_I^0\right],g_I^0>0$. Put

$\tilde{X}\left(y\right)=U\tilde{Y}{U}^{-1}\left(y\right)$ (30)

Let

$X\left(y,v\right)=\tilde{X}\left(y\right)+e_{3}v=UY\left(U^{-1}\left(y\right),v\right)$ (31)

Define $d_1\mathrm{,}{d}_2\mathrm{,}{d}_3$ by

$U^{-1}{\left(g\right)}_{\mathrm{1,2}}=L_1^{-1}\left(L-\text{id}\right)d_{\mathrm{1,2}}$ (32)

and

$\frac{\mu }{\ln \left(1+\mu \right)}{d}_3={\left(U^{-1}g\right)}_3$ (33)

Then

${\Phi }_1^X=T$ (34)

when $v=0$. To this vector field with $v\in \left[\mathrm{0,}{g}_I^0\right]$ associate the Hamiltonian

$H\left(x\mathrm{,}p\mathrm{,}v\mathrm{,}t\right)=p^{\text{T}}Y\left(x\mathrm{,}v\right)$ (35)

Then we have the adjoint equations

${p^{\prime }}_1=-\frac{\partial H}{\partial x_1}=-p_1{a}_1+p_2{b}_1$ (36)

${p^{\prime }}_2=-\frac{\partial H}{\partial x_2}=-p_1{b}_1-p_2{a}_1$ (37)

${p^{\prime }}_3=-\frac{\partial H}{\partial x_3}=-p_3\ln \left(1+\mu \right)$ (38)

So

${p^{\prime }}_{1,2}=\left(\begin{array}{cc}-a_1& b_1\\ -b_1& -a_1\end{array}\right)\left(\begin{array}{c}{p}_1\\ p_2\end{array}\right)$ (39)

which has flow

$p_{\mathrm{1,2}}\left(t\right)=\exp \left(-a_{1}t\right)\left(\begin{array}{cc}\cos \left(b_{1}t\right)& \sin \left(b_{1}t\right)\\ -\sin \left(b_{1}t\right)& \cos \left(b_{1}t\right)\end{array}\right)p_{\mathrm{1,2}}\left(0\right)$ (40)

Define

$d\left(t\right)=L_1^{-1}\left(\exp \left(L_{1}t\right)-\text{id}\right)d_{\mathrm{1,2}}$ (41)

Then the flow of Y is for $v=0$

${\Phi }^Y{\left(t\mathrm{,}x\right)}_{\mathrm{1,2}}=\exp \left(L_{1}t\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)+d\left(t\right)$ (42)

${\Phi }^Y{\left(t\mathrm{,}x\right)}_3=\exp \left(\ln \left(1+\mu \right)t\right)x_3+\frac{\exp \left(\ln \left(1+\mu \right)t\right)-1}{\ln \left(1+\mu \right)}{d}_3$ (43)

Define

$S_1\left(x\right)=\left(1+\mu -a\right)x_1-b{x}_2$ (44)

Then we have the transversality conditions

$p_1\left(T\right)=\frac{\partial S_1}{\partial x_1}=1+\mu -a$ (45)

$p_2\left(T\right)=\frac{\partial S_1}{\partial x_2}=-b$ (46)

$p_3\left(T\right)=\frac{\partial S_1}{\partial x_3}=0$ (47)

This means that $p_3\left(t\right)=0$ and

$p_{1,2}\left(T\right)=\exp \left(-a_{1}T\right)\left(\begin{array}{cc}\cos \left(b_{1}T\right)& \sin \left(b_{1}T\right)\\ -\sin \left(b_{1}T\right)& \cos \left(b_{1}T\right)\end{array}\right)p_{1,2}\left(0\right)=\left(\begin{array}{c}1+\mu -a\\ -b\end{array}\right)$ (48)

The two by two matrix in this equation has inverse

$\left(\begin{array}{cc}\cos \left(b_{1}T\right)& -\sin \left(b_{1}T\right)\\ \sin \left(b_{1}T\right)& \cos \left(b_{1}T\right)\end{array}\right)$ (49)

So

$p_{1,2}\left(0\right)=\exp \left(a_{1}T\right)\left(\begin{array}{cc}\cos \left(b_{1}T\right)& -\sin \left(b_{1}T\right)\\ \sin \left(b_{1}T\right)& \cos \left(b_{1}T\right)\end{array}\right)\left(\begin{array}{c}1+\mu -a\\ -b\end{array}\right)$ (50)

But then we have

$p_{1,2}\left(t\right)=\exp \left(a_1\left(T-t\right)\right)\left(\begin{array}{cc}\cos \left(b_1\left(t-T\right)\right)& \sin \left(b_1\left(t-T\right)\right)\\ -\sin \left(b_1\left(t-T\right)\right)& \cos \left(b_1\left(t-T\right)\right)\end{array}\right)\left(\begin{array}{c}1+\mu -a\\ -b\end{array}\right)$ (51)

Notice that we have

$U_{31}=b\beta$ (52)

$U_{32}=\left(1+\mu -a\right)\beta$ (53)

If

$H\left(x^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}p\left(t\right)\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}t\right)\le H\left(x^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}p\left(t\right)\mathrm{,}u\mathrm{,}t\right)$ (54)

for all $u\in \mathbb{U}=\left[\mathrm{0,}{g}_I^0\right]$, then $\left(x^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\right)$ is optimal, see Equations (55) to (60). But this amounts to the inequality

$\left(p_1\left(t\right)b\beta u_{*}+p_2\left(t\right)\left(1+\mu -a\right)\beta u_{*}\right)\frac{1}{\det \left(U\right)}$ (55)

$=\exp \left(a_1\left(T-t\right)\right)(\left(\cos \left(b_1\left(t-T\right)\right)\left(1+\mu -a\right)-b\sin \left(b_1\left(t-T\right)\right)\right)b\beta u_{*}$ (56)

$+\left(-\sin \left(b_1\left(t-T\right)\right)\left(1+\mu -a\right)-b\cos \left(b_1\left(t-T\right)\right)\right)\left(\left(1+\mu -a\right)\beta u_{*}\right))\frac{1}{\det \left(U\right)}$ (57)

$=\exp \left(a_1\left(T-t\right)\right)\sin \left(b_1\left(t-T\right)\right)\left(-b^2-{\left(1+\mu -a\right)}^2\right)\beta u_{*}\frac{1}{\det \left(U\right)}$ (58)

$=\exp \left(a_1\left(T-t\right)\right)\left(\alpha \delta +\beta \sigma \right)\beta u_{*}\sin \left(b_1\left(t-T\right)\right)\frac{1}{\det \left(U\right)}$ (59)

$\le \exp \left(a_1\left(T-t\right)\right)\sin \left(b_1\left(t-T\right)\right)\frac{\beta u}{-b}$ (60)

where $\text{det}\left(U\right)$ denotes the determinant of U. We have used that

$1+\mu -a=1+\mu -\left(1+\frac{\gamma +\mu }{2}\right)=\frac{\mu -\gamma }{2}$ (61)

and

${\lambda }_{+}=1+\frac{\gamma +\mu }{2}+\frac{1}{2}\sqrt{{\left(\gamma -\mu \right)}^2+4\left(\alpha \delta +\beta \sigma \right)}=a+ib$ (62)

so that

${\left(1+\mu -a\right)}^2+b^2=-\left(\alpha \delta +\beta \sigma \right)$ (63)

and also that $\det U=-b\left(\alpha \delta +\beta \sigma \right)$. So if

$\frac{\sin \left(b_1\left(t-T\right)\right)}{-b}>0$ (64)

then

$u^{*}\left(t\right)=g_I^0$ (65)

is optimal and if the reverse inequality holds then

$u^{*}\left(t\right)=0$ (66)

is optimal. We shall now consider the case where all eigenvalues ${\lambda }_{+}\mathrm{,}{\lambda }_{-}\mathrm{,}\tilde{\lambda }$ are real, positive and distinct and that $\mu ={\mu }_F={\mu }_I$. Here we let

$Y\left(x,u\right)=\tilde{\Lambda }x+d+D^{-1}{e}_{3}u$ (67)

where

$\tilde{\Lambda }=\text{diag}\left(\ln {\lambda }_{+},\ln {\lambda }_{-},\ln \left(1+\mu \right)\right)$ (68)

Also define d in

$\frac{{\lambda }_{+}-1}{\ln {\lambda }_{+}}{d}_1={\left(D^{-1}g\right)}_1$ (69)

$\frac{{\lambda }_{-}-1}{\ln {\lambda }_{+}}{d}_2={\left(D^{-1}g\right)}_2$ (70)

$\frac{\mu }{\ln \left(1+\mu \right)}{d}_3={\left(D^{-1}g\right)}_3$ (71)

$d={\left(d_1\mathrm{,}{d}_2\mathrm{,}{d}_3\right)}^{\text{T}}\in {ℝ}^3$. Now define the Hamiltonian

$H\left(x,p,u,t\right)=p^{\text{T}}Y\left(x,u\right)$ (72)

Then we get the adjoint equations

${p^{\prime }}_1=-\frac{\partial H}{\partial x_1}=-\ln {\lambda }_{+}{p}_1$ (73)

${p^{\prime }}_2=-\frac{\partial H}{\partial x_2}=-\ln {\lambda }_{-}{p}_2$ (74)

${p^{\prime }}_3=-\frac{\partial H}{\partial x_3}=-\ln \left(1+\mu \right)p_3$ (75)

see [2]. Now define

$S_1\left(x\right)=\left(1+\mu -{\lambda }_{+}\right)x_1+\left(1+\mu -{\lambda }_{-}\right)x_2$ (76)

Then we have the transversality conditions

$p\left(T\right)=\frac{\partial S_1}{\partial x}=\left(\begin{array}{c}1+\mu -{\lambda }_{+}\\ 1+\mu -{\lambda }_{-}\\ 0\end{array}\right)$ (77)

by [2] and see below. Now observe, that

$\det \left(D\right)=\left({\lambda }_{+}-{\lambda }_{-}\right)\left(\alpha \delta +\beta \sigma \right)$ (78)

Because

$1+\mu -{\lambda }_{\pm }=\frac{\mu -\gamma }{2}\mp \frac{1}{2}\sqrt{{\left(\mu -\gamma \right)}^2+4\left(\alpha \delta +\beta \sigma \right)}$ (79)

we find

$\left(1+\mu -{\lambda }_{+}\right)\left(1+\mu -{\lambda }_{-}\right)=-\left(\alpha \delta +\beta \sigma \right)$ (80)

We shall need

$D^{-1}=\left(\begin{array}{ccc}-\alpha \delta -\beta \sigma & -\alpha \left(1+\mu -{\lambda }_{-}\right)& -\beta \left(1+\mu -{\lambda }_{-}\right)\\ \alpha \delta +\beta \sigma & \alpha \left(1+\mu -{\lambda }_{+}\right)& \beta \left(1+\mu -{\lambda }_{+}\right)\\ 0& -\sigma \left({\lambda }_{+}-{\lambda }_{-}\right)& \delta \left({\lambda }_{+}-{\lambda }_{-}\right)\end{array}\right)\frac{1}{\det D}$ (81)

We have

$\left(\begin{array}{c}C\left(T\right)\\ 0\\ 0\end{array}\right)=D\left(\begin{array}{c}{x}_1\left(T\right)\\ x_2\left(T\right)\\ x_3\left(T\right)\end{array}\right)$ (82)

hence the definition of $S_1$. Now

$p\left(T\right)=\frac{\partial S_1}{\partial x}$ (83)

thus

$p_1\left(T\right)={p_1\left(0\right)\exp \left(-\ln {\lambda }_{+}t\right)|}_{t=T}=p_1\left(0\right)\exp \left(-\ln {\lambda }_{+}T\right)=1+\mu -{\lambda }_{+}$ (84)

$p_2\left(T\right)={p_2\left(0\right)\exp \left(-\ln {\lambda }_{-}t\right)|}_{t=T}=p_2\left(0\right)\exp \left(-\ln {\lambda }_{-}T\right)=1+\mu -{\lambda }_{-}$ (85)

$p_3\left(T\right)={p_3\left(0\right)\exp \left(-\ln \left(1+\mu \right)t\right)|}_{t=T}=p_3\left(0\right)\exp \left(-\ln \left(1+\mu \right)T\right)=0$ (86)

So

$p_1\left(t\right)=\left(1+\mu -{\lambda }_{+}\right)\exp \left(\ln {\lambda }_{+}\left(T-t\right)\right)$ (87)

$p_2\left(t\right)=\left(1+\mu -{\lambda }_{-}\right)\exp \left(\ln {\lambda }_{-}\left(T-t\right)\right)$ (88)

$p_3\left(t\right)=0$ (89)

If

$H\left(x^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}p\left(t\right)\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}t\right)\le H\left(x^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}p\left(t\right)\mathrm{,}u\mathrm{,}t\right)$ (90)

for all $u\in \mathbb{U}$,$\left(x^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\right)$ is optimal. And this is equivalent to

$p{\left(t\right)}^{\text{T}}{D}^{-1}{e}_3{u}^{*}\left(t\right)\le p{\left(t\right)}^{\text{T}}{D}^{-1}{e}_{3}u$ (91)

which again is equivalent to

$\left(1+\mu -{\lambda }_{+}\right)\exp \left(\ln {\lambda }_{+}\left(T-t\right)\right)\left(-\beta \right)\left(1+\mu -{\lambda }_{-}\right)\frac{u^{*}\left(t\right)}{\det \left(D\right)}$ (92)

$+\left(1+\mu -{\lambda }_{-}\right)\exp \left(\ln {\lambda }_{-}\left(T-t\right)\right)\beta \left(1+\mu -{\lambda }_{+}\right)\frac{u^{*}\left(t\right)}{\det \left(D\right)}$ (93)

$=u^{*}\left(t\right)\beta \frac{-\alpha \delta -\beta \sigma }{\left({\lambda }_{+}-{\lambda }_{-}\right)\left(\alpha \delta +\beta \sigma \right)}\left(-\exp \left(\ln {\lambda }_{+}\left(T-t\right)\right)+\exp \left(\ln {\lambda }_{-}\left(T-t\right)\right)\right)$ (94)

$\le u\beta \left(\exp \left(\ln {\lambda }_{+}\left(T-t\right)\right)-\exp \left(\ln {\lambda }_{-}\left(T-t\right)\right)\right)\frac{1}{{\lambda }_{+}-{\lambda }_{-}}$ (95)

for all $u\in \mathbb{U}$. It follows that

$u^{*}\left(t\right)=g_I^0$ (96)

is optimal. We have the two Hamiltonians

$H^Y\left(x,p,u,t\right)=p^{\text{T}}Y\left(x,u\right)=\langle p,Y\left(x,u\right)\rangle$ (97)

$H^X\left(y,q,v,t\right)=q^{\text{T}}X\left(y,v\right)=\langle q,X\left(y,v\right)\rangle$ (98)

where $\langle \cdot \mathrm{,}\cdot \rangle$ denotes the canonical inner product. It follows that when the eigenvalues are $a\pm ib,1+\mu >0,a>0,b\ne 0$

$y\left(t\right)=Ux\left(t\right)$ (99)

$p\left(t\right)=U^{\text{T}}q\left(t\right)$ (100)

if

$y\left(0\right)=Ux\left(0\right)$ (101)

$p\left(T\right)=U^{\text{T}}q\left(T\right)$ (102)

where $y\left(t\right)$ is an integral curve of X and $x\left(t\right)$ an integral curve of Y. Because

$q^{\prime }=-\frac{\partial H^X}{\partial y}=-{\left(U^{-1}\right)}^{\text{T}}{\tilde{B}}^{\text{T}}{U}^{\text{T}}q$ (103)

hence

$U^{\text{T}}{q}^{\prime }=-{\tilde{B}}^{\text{T}}{U}^{\text{T}}q$ (104)

and

$p^{\prime }=-{\tilde{B}}^{\text{T}}p$ (105)

We now get, that

$H^X\left(y\left(t\right),q\left(t\right),v\left(t\right),t\right)=\langle q\left(t\right),X\left(y\left(t\right),v\left(t\right)\right)\rangle$ (106)

$=\langle q\left(t\right),UY\left(U^{-1}y\left(t\right),v\left(t\right)\right)\rangle$ (107)

$=\langle U^{\text{T}}\left(q\left(t\right)\right),Y\left(U^{-1}\left(y\left(t\right)\right),v\left(t\right)\right)\rangle$ (108)

$=\langle p\left(t\right)\mathrm{,}Y\left(x\left(t\right)\mathrm{,}v\left(t\right)\right)\rangle$ (109)

$=H^Y\left(x\left(t\right),p\left(t\right),v\left(t\right),t\right)$ (110)

So

$H^X\left(y\left(t\right),q\left(t\right),v\left(t\right),t\right)=H^Y\left(x\left(t\right),p\left(t\right),v\left(t\right),t\right)$ (111)

When the eigenvalues of A are real, distinct and positive

$y\left(t\right)=Dx\left(t\right)$ (112)

$p\left(t\right)=D^{\text{T}}q\left(t\right)$ (113)

if

$y\left(0\right)=Dx\left(0\right)$ (114)

$p\left(T\right)=D^{\text{T}}q\left(T\right)$ (115)

We also have

$q^{\prime }=-{\left(D^{-1}\right)}^{\text{T}}{\Lambda }^{\text{T}}{D}^{\text{T}}q$ (116)

and

${\left(D^{\text{T}}q\right)}^{\prime }=-{\Lambda }^{\text{T}}{D}^{\text{T}}q$ (117)

thus

$p^{\prime }=-{\Lambda }^{\text{T}}p$ (118)

Hence

$H^X\left(y\left(t\right),q\left(t\right),v\left(t\right),t\right)=H^Y\left(x\left(t\right),p\left(t\right),v\left(t\right),t\right)$ (119)

We need the following theorem, which is well known.

Theorem 1 Now The following statements about a C² function from

$f\mathrm{<mo>\; :\; </mo>}{ℝ}^n\to ℝ$ (120)

to the reals, where $n$ is a positive integer, $n\in \mathbb{N}$ are equivalent:

(i) $f\left(\lambda x+\left(1-\lambda \right)y\right)\le \lambda f\left(x\right)+\left(1-\lambda \right)f\left(y\right)$ where $\lambda \in \left[\mathrm{0,1}\right]$ ;

(ii) $f\left(y\right)\ge f\left(x\right)+\frac{\partial f}{\partial x}\left(x\right)\left(y-x\right)$ ;

(iii) ${\displaystyle {\sum }_{i,j=1,\cdots ,n}}{\gamma }_i\frac{{\partial }^{2}f}{\partial x_i\partial x_j}{\gamma }_j\ge 0$

where $x\mathrm{,}y\mathrm{,}\gamma \in {ℝ}^n$.

$\left(x\left(t\right)\mathrm{,}u\left(t\right)\right)$ is admissible by definition if $0\le u\left(t\right)\le g_I^0$ and

$x^{\prime }\left(t\right)=Y\left(x\left(t\right)\mathrm{,}u\left(t\right)\right),x\left(0\right)\text{=}{x}_0\in {ℝ}^3$ (121)

To see that

$S_1\left(x^{\mathrm{<mo>\; *\; </mo>}}\left(T\right)\right)\le S_1\left(x\left(T\right)\right)$ (122)

for $\left(x^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\right)$ optimal candidate and $\left(x\left(t\right)\mathrm{,}u\left(t\right)\right)$ admissible argue as in [2]

$\Delta =S_1\left(x^{*}\left(T\right)\right)-S_1\left(x\left(T\right)\right)$ (123)

$=S_1\left(x^{*}\left(T\right)\right)-S_1\left(x\left(T\right)\right)$ (124)

$\begin{array}{l}+{\displaystyle {\int }_0^T}\left(H^Y\left(x^{*}\left(t\right),p\left(t\right),u^{*}\left(t\right),t\right)-p{\left(t\right)}^{\text{T}}{x}^{*}{}^{\prime }\left(t\right)\right)\\ -\left(H^Y\left(x\left(t\right),p\left(t\right),u\left(t\right),t\right)-p{\left(t\right)}^{\text{T}}{x}^{\prime }\left(t\right)\right)\text{d}t\end{array}$ (125)

We have the following inequality, which follows from theorem 1.

$\begin{array}{l}{H}^Y\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}p\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}t\right)-H^Y\left(x\mathrm{,}p\mathrm{,}u\mathrm{,}t\right)\\ \le \frac{\partial H^Y}{\partial x}\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}p\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}t\right)\left(x^{\mathrm{<mo>\; *\; </mo>}}-x\right)+\frac{\partial H^Y}{\partial u}\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}p\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}t\right)\left(u^{\mathrm{<mo>\; *\; </mo>}}-u\right)\end{array}$ (126)

We also have

$p^{\prime }=-\frac{\partial H^Y}{\partial x}\left(x^{*},p,u^{*},t\right)$ (127)

We can now estimate

$\Delta \le S_1\left(x^{\mathrm{<mo>\; *\; </mo>}}\left(T\right)\right)-S_1\left(x\left(T\right)\right)+{\displaystyle {\int }_0^T}{p^{\prime }}^{\text{T}}\left(x-x^{\mathrm{<mo>\; *\; </mo>}}\right)+p^{\text{T}}\left(x^{\prime }-x^{\mathrm{<mo>\; *\; </mo>}}{}^{\prime }\right)\text{d}t$ (128)

$\begin{array}{l}+{\displaystyle {\int }_0^T}\frac{\partial H^Y}{\partial u}\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}p\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}t\right)\left(u^{\mathrm{<mo>\; *\; </mo>}}-u\right)\text{d}t\\ \le {\displaystyle {\int }_0^T}\frac{\text{d}}{\text{d}t}\left(p^{\text{T}}\left(x-x^{\mathrm{<mo>\; *\; </mo>}}\right)\right)\left(t\right)\text{d}t+S_1\left(x^{\mathrm{<mo>\; *\; </mo>}}\left(T\right)\right)-S_1\left(x\left(T\right)\right)\end{array}$ (129)

$=p{\left(T\right)}^{\text{T}}\left(x\left(T\right)-x^{\mathrm{<mo>\; *\; </mo>}}\left(T\right)\right)+S_1\left(x^{\mathrm{<mo>\; *\; </mo>}}\left(T\right)\right)-S_1\left(x\left(T\right)\right)$ (130)

$=\frac{\partial S_1}{\partial x}\left(x^{\mathrm{<mo>\; *\; </mo>}}\left(T\right)\right)\left(x\left(T\right)-x^{\mathrm{<mo>\; *\; </mo>}}\left(T\right)\right)+S_1\left(x^{\mathrm{<mo>\; *\; </mo>}}\left(T\right)\right)-S_1\left(x\left(T\right)\right)$ (131)

$\le 0$ (132)

because $S_1$ is convex, by theorem 1. We have used, that we have arranged, that

$\frac{\partial H}{\partial u}\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}p\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}t\right)\left(u^{\mathrm{<mo>\; *\; </mo>}}-u\right)\le 0$ (133)

for all $u\in \mathbb{U}$, by the mean value theorem.

We have optimality.

3. Optimal Control of T

In this section we consider the problem: minimize $C\left(N\right)$ subject to

$y_{k+1}=A{y}_k+B{u}_k+g=f\left(y_k,u_k\right)$ (134)

$k=1,\cdots ,N-1,N\in \mathbb{N},N\ge 2,u\left(k\right)\in \mathbb{U}=\left[0,g_I^0\right],g_i^0>0,\mu ={\mu }_F={\mu }_I$ where A is as in the introduction.

Also

$y_k=\left(C\left(k\right)\mathrm{,}GF\left(k\right)\mathrm{,}GI\left(k\right)\right)\mathrm{,}g\in {ℝ}^3$ (135)

Here

$B={\left(0,0,1\right)}^{\text{T}},g={\left(g_C,g_F,g_I\right)}^{\text{T}}$ (136)

Assume (i) the eigenvalues of A are real and distinct.

In the Discrete Pontryagin Minimal Principle applied to T you define the Hamiltonian by (138) and then you find

$\frac{\partial H}{\partial u_k}\left(x_k^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{u}_k^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\lambda }_k\right)\left(u_k^{\mathrm{<mo>\; *\; </mo>}}-u_k\right)$ (137)

and minimize it to find the optimal control $u_k^{\mathrm{<mo>\; *\; </mo>}}$. It is optimal due to computations (157) to (163) below.

Define then the Hamiltonian

$H\left(x_k,u_k,{\lambda }_k\right)={\lambda }_k^{\text{T}}\left(\tilde{A}{x}_k+\hat{B}{u}_k+D^{-1}g\right)$ (138)

where ${\lambda }_k\in {ℝ}^3$ and (139)

$\hat{B}=D^{-1}B$ (140)

$\tilde{A}=D^{-1}AD$ (141)

Then we have the adjoint equation

${\lambda }_{k-1}=\frac{\partial H}{\partial x_k}\left(x_k,u_k,{\lambda }_k\right)={\tilde{A}}^{\text{T}}{\lambda }_k$ (142)

Inductively

${\lambda }_{N-k-1}={\left({\tilde{A}}^{\text{T}}\right)}^k{\lambda }_{N-1}$ (143)

In particular

${\lambda }_0={\left({\tilde{A}}^{\text{T}}\right)}^{N-1}{\lambda }_{N-1}$ (144)

For $k=0$ we have

$\begin{array}{l}H\left(x_0^{*},u_0^{*},{\lambda }_0\right)={\left({\lambda }_0\right)}^{\text{T}}\left(\tilde{A}{x}_0^{*}+\hat{B}{u}_0^{*}+D^{-1}g\right)\\ \le H\left(x_0^{*},u,{\lambda }_0\right)={\left({\lambda }_0\right)}^{\text{T}}\left(\tilde{A}{x}_0^{*}+\hat{B}u+D^{-1}g\right)\end{array}$ (145)

which is equivalent to

${\left({\lambda }_0\right)}^{\text{T}}\hat{B}{u}_0^{\mathrm{<mo>\; *\; </mo>}}\le {\left({\lambda }_0\right)}^{\text{T}}\hat{B}u$ (146)

Define

$S_1\left(x\right)=F\left(x\right)=\left(1+\mu -{\lambda }_{+}\right)x_1+\left(1+\mu -{\lambda }_{-}\right)x_2$ (147)

Now

${\lambda }_0={\left({\tilde{A}}^{\text{T}}\right)}^{N-1}{\lambda }_{N-1}$ (148)

where

${\lambda }_{N-1}=\frac{\partial F}{\partial x}\left(x_N\right)$ (149)

Thus

${\left({\lambda }_0\right)}^{\text{T}}\hat{B}{u}_0^{*}={\left({\left({\tilde{A}}^{\text{T}}\right)}^{N-1}{\lambda }_{N-1}\right)}^{\text{T}}\hat{B}{u}_0^{*}={\left({\lambda }_{N-1}\right)}^{\text{T}}{\tilde{A}}^{N-1}\hat{B}{u}_0^{*}={\frac{\partial F}{\partial x}}^{\text{T}}\left(x_N\right){\tilde{A}}^{N-1}\hat{B}{u}_0^{*}$ (150)

Assume that (i) holds and ${\lambda }_{+},{\lambda }_{-},1+\mu \in ℝ$ are distinct, when $\mu ={\mu }_F={\mu }_I$. Then

$\hat{B}=\left(\begin{array}{c}{D}_{31}\\ D_{32}\\ D_{33}\end{array}\right)\frac{1}{\det \left(D\right)}=\left(\begin{array}{c}-\beta \left(1+\mu -{\lambda }_{-}\right)\\ \beta \left(1+\mu -{\lambda }_{+}\right)\\ \delta \left({\lambda }_{+}-{\lambda }_{-}\right)\end{array}\right)\frac{1}{\det \left(D\right)}=D^{-1}{e}_3$ (151)

So now we get, that (146) amounts to

$\left(1+\mu -{\lambda }_{+}\mathrm{,1}-\mu -{\lambda }_{-}\mathrm{,0}\right)\left(\begin{array}{ccc}{\lambda }_{+}^{N-1}& 0& 0\\ 0& {\lambda }_{-}^{N-1}& 0\\ 0& 0& {\left(1+\mu \right)}^{N-1}\end{array}\right)\hat{B}$ (152)

$=\frac{-\beta \left(-\alpha \delta -\beta \sigma \right){\lambda }_{+}^{N-1}}{\det \left(D\right)}+\frac{\beta \left(-\alpha \delta -\beta \sigma \right){\lambda }_{-}^{N-1}}{\det \left(D\right)}$ (153)

$=\frac{\alpha \delta +\beta \sigma }{\left({\lambda }_{+}-{\lambda }_{-}\right)\left(\alpha \delta +\beta \sigma \right)}\beta \left({\lambda }_{+}^{N-1}-{\lambda }_{-}^{N-1}\right)$ (154)

$=\frac{\beta }{{\lambda }_{+}-{\lambda }_{-}}\left({\lambda }_{+}^{N-1}-{\lambda }_{-}^{N-1}\right)\le 0$ (155)

Similarly for $k=1,\cdots ,N-2$

${\lambda }_k^{\text{T}}\hat{B}{u}_k^{*}=\left(1+\mu -{\lambda }_{+},1+\mu -{\lambda }_{-},0\right)\left(\begin{array}{ccc}{\lambda }_{+}^{N-k-1}& 0& 0\\ 0& {\lambda }_{-}^{N-k-1}& 0\\ 0& 0& {\left(1+\mu \right)}^{N-k-1}\end{array}\right)\hat{B}{u}_k^{*}<0$ (156)