###### More Information

**Submitted:** 22 October 2019 | **Approved:** 09 December 2019 | **Published:** 10 December 2019

**How to cite this article:**: Sweilam NH, Tharwat AA, Abd El Moniem NK. Different optimization strategies for the optimal control of tumor growth. Arch Cancer Sci Ther. 2019; 3: 052-062.

**DOI:** 10.29328/journal.acst.1001010

**Copyright:** © 2019 Sweilam NH, et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**Keywords:** Optimal control for differential equations; Tumor-immune model; Optimal control direct methods; Opti-mal control indirect methods; Nonlinear programming

# Different optimization strategies for the optimal control of tumor growth

####
Sweilam NH^{1}, Tharwat AA^{2} and Abd El Moniem NK^{3}*

^{1}Department of Mathematics, Cairo University, Faculty of Science, Giza, Egypt

^{2}Department of Management and Marketing, Colleague of BA, American University at UAE

^{3}Department of Information Technology, Cairo University, National Cancer Institute, Egypt

***Address for Correspondence:** Abd El Moniem NK, Department of Information Technology, Cairo University, National Cancer Institute, Egypt, Email: nermeen2000@hotmail.com; nermeen.elzaweli@nci.cu.edu.eg

In this article different numerical techniques for solving optimal control problems is introduced, the aim of this paper is to achieve the best accuracy for the Optimal Control Problem (OCP) which has the objective of minimizing the size of tumor cells by the end of the treatment. An important aspect is considered, which is, the optimal concentrations of drugs that not affect the patient’s health significantly. To study the behavior of tumor growth, a mathematical model is used to simulate the dynamic behavior of tumors since it is difficult to prototype dynamic behavior of the tumor. A tumor-immune model with four components, namely, tumor cells, active cytotoxic T-cells (CTLs), helper T-cells, and a chemotherapeutic drug is used. Two general categories of optimal control methods which are indirect methods and direct ones based on nonlinear programming solvers and interior point algorithms are compared. Within the direct optimal control techniques, we review three different solutions techniques namely (i) multiple shooting methods, (ii) trapezoidal direct collocation method, (iii) Hermit- Simpson’s collocation method and within the indirect methods we review the Pontryagin’s Maximum principle with both collocation method and the backward forward sweep method. Results show that the direct methods achieved better control than indirect methods.

Mutations is the major consequence of abnormal behavior of the genetic material DNA. Mutations affect normal growth and division i.e., mutations cause uncontrollable growth of cells, these mutations are caused by two reasons which are the signals telling cells to begin dividing are left on continuously or growth suppressing signals telling cells not divide are turned off [19]. The process starts from an evolutionary process which may give rise to abnormal DNA when a cell duplicates its genome due to defects in tumor suppressor or DNA mismatch repair genes. Tumors releases hormones that alter body function, they can grow and interfere with the digestive, nervous and circulatory systems, they can also invade nearly tissues and successfully spread to other parts of the body and grow, these tumors are of malignant type, whereas benign tumors are not cancerous and not spread to other parts of the body.

The immune system is composed of a wide variety of cells with different functions, most important types of them are innate, or nonspecific and specific, or adaptive, immune response [1], the similarity between innate and adaptive responses is that both must contact with target tumor cells in order to be able to kill them, whereas the difference between them is that the innate cells are always on patrol, and kill tumor cells not recognized them as self, so it is an early defense against pathogens. While adaptive, immune response must be primed to recognize antigen specific to the tumor cells. Cell-mediated immunity, also called cellular immunity, is mediated by T lymphocytes (also called T cells) are part of the adaptive immune response. as well as the Killer T cells are referred to as CTL “Cytotoxic T Lymphocyte” cells, or CD8+ T cells, which directly attack and eliminate infected cells, while, the CD4 helper T cells, assist other cells of the immune system during an infection.

The immune system alone usually fails to effectively fight the tumor for the following reasons:

The first is an insufficient immune response due to the tumor being poorly antigenic. This is described by Curilas [13] “too little of a good thing” Cytotoxic T cells are not sufficiently activated by the tumor, and therefore the response is minimal.

The second cause for failure is “too much of a bad thing” in which immunosuppressive factors damage otherwise capable immune system. Over the last few decades, the importance of the second paradigm has become clear; immune suppression is likely to be a significant factor when cancer is present in the host.

Both mathematical modeling and experimental “clinical” results are used to explore the significance of tumor–immune inter- actions. Previous work has included models of T cells with interleukin-2 (IL-2) [24], transforming growth factor beta (TGF-b) [25], regulatory T cells (Tregs) [28], and natural killer cells [14].

Several answers to the question how to choose best possible dose to reduce harm to healthy cells while beat cancer are introduced. These answers was introduced before based on empirical methods, empirical methods depend basically on drug holidays, which are rest periods give the health cells the opportunity to recover from the toxic attack of the drug these periods are specified by try and error methods. Apart from the empirical methods the mathematical modeling and optimal control can answer this question, for more details see [35].

In this article a base line model [37] is used based on tumor-immune model. A comparative study between the results of the Computation method used in our baseline article [37] which uses indirect method resulted in two point boundary value problem (TPBVP) and solved by collocation method is compared with another method also under the category of indirect method which also resulted in TPBVP but solved by other method named forward backward sweep method (FBWM), results of the base line article [37] also compared with other three different computational methods under the category of direct methods.

Previous comparative studies between direct and indirect control given in [38], where this paper two different approaches (direct and indirect) for the numerical solution of fractional optimal control problems (FOCPs) based on a spectral method using Chebyshev polynomial are presented. Moreover, in [18] a comparative study between singular arc method (without con- sidering inequality constraints) as indirect method and three different direct methods namely Hermite-Simpson’s collocation method, 5^{th} degree Gauss-Lobatto collocation method, Radau Psuedospectral collocation method using GPOP [34]. In [9] a comparative study between dynamic programming, indirect methods and direct methods is presented.

This article is organized as follows: in section 2 we consider a tumor growth model exhibiting the effect of tumor–immune interaction with chemotherapeutic. The model construction and assumptions is introduced in “Mathematical Model” section. In section 3, the optimization stratigetis is presented. In section 4, the solution of OC problem is discussed. In section 5 a comparative study for solving OC problem is presented. The conclusion is given is section 6.

##### Mathematical modeling and optimal control problem for a tumor-immune system

In this section, a mathematical model problem is introduced, for more details see [37] as a description of the phenomena of tumor-immune interaction and the prediction of the outcome of the application of chemotherapy regimen is examined by using optimal control [5,33]. We can consider optimal control problem as a type of optimization problem where the objective is to determine the inputs (equivalently, the trajectory, state or path), the control inputs (equivalently, the trajectory, state or path), the control input *u*∈(*t*) *R*^{m}, the initial time ${t}_{0}\in R$
and the terminal time, ${t}_{f}\in R$
(where $t\in \left[{t}_{0},{t}_{f}\right]$
the independent variable) of the dynamical system [34] optimize (i.e., minimize or maximize) a speciﬁed objective function (performance index) while satisfying any constraints on the motion of the system.

the Objective function (performance index) is represented by:

$J({u}^{\ast})\text{}=min\{J\left(u\left(t\right)\right):u\in U\},$

Where

$J\left(u\left(t\right)\right)={\displaystyle \underset{0}{\overset{{t}_{f}}{\int}}\left[{B}_{1}T+\frac{1}{2}{B}_{2}{u}^{2}\right]}dt\text{(2}\text{.1)}$

where *B _{1},Bsub>1* are positive constants representing the weights of the terms. The ﬁrst term represents the tumor cell populations and the second term represents the harmful effects of drug on body. The square of the control variable $\left({u}^{2}\left(t\right)\right)$
reﬂects the severity of the side effects of the drug imposed, for more details see [22,42] and the references cited therein. When chemotherapeutic drugs are administered in high dose, they are toxic to the human body, which justiﬁes the quadratic terms in the functional. So the functional given in Eq. (2.1) should be minimized.

The dynamical system is deﬁned by a set of ordinary differential equations (ODE’s):

$\frac{dT}{dt}={\text{r}}_{\text{1}}\text{T}\left(\text{1}-{\text{p}}_{\text{1}}\text{T}\right)-\alpha {\text{TI}}_{\text{H}}-{\text{q}}_{\text{1}}\text{DT,}$ $\frac{d{I}_{H}}{dt}=\beta {\text{I}}_{\text{H}}{\text{I}}_{\text{R}}-{\alpha}_{\text{2}}{\text{TI}}_{\text{H}}-{\text{dI}}_{\text{H}}-{\text{q}}_{\text{2}}{\text{DI}}_{\text{H}},$ $\frac{d{I}_{R}}{dt}={\text{r}}_{\text{2}}{\text{I}}_{\text{R}}\left(\text{1}-{\text{p}}_{\text{2}}{\text{I}}_{\text{R}}\right)-\beta {\text{I}}_{\text{H}}{\text{I}}_{\text{R}}-{\text{q}}_{\text{3}}{\text{DI}}_{\text{R}},$ $\frac{dD}{dt}=u(t)-\gamma D,\text{(2}\text{.2)}$

$\text{satisfying}\text{\hspace{1em}}\text{\hspace{0.22em}}T\left(0\right)\ge {T}_{0},\text{\hspace{0.17em}}{I}_{H}\left(0\right)\ge 0,\text{\hspace{0.17em}}{I}_{R}\left(0\right)\ge 0,\text{\hspace{0.17em}}D\left(0\right)\ge 0.$

Where

• *T* (*t*) is the numbers of tumor cell.

• I*I _{H}* (

*t*) is the active CTL cells (hunting CTL cells).

• I*I _{R}*(

*t*) is the helper T-cells (resting T-cells).

*D*(*t*) is the density of chemotherapeutic drug at time *t*.

The tumor-immune model of the base line paper is originally developed by de Pillis and Radunskaya [31], an optimal control for this model is introduced in [15]. This model considers interactions between tumor cell population and two types of immune cell populations which are helper (resting) T-cells and active (hunting) CTL cells. The damage done to each cell due to chemotherapy is subtracted from each cell population. According to [37] which is based on [14,15] the model assumptions are as follows:

• Hunting CTL cells are capable of killing tumor cells, the effect of CTL on tumor cells in represented by both the terms:

1. *αTI _{H}* which represents the loss of tumor cells it is proportional to the product of the of densities of tumor cells and active CTL cell (hunting) which is subtracted in the equation represents rate of change of tumor population $\frac{dT}{dt}$

2. *αT2 _{H}* is the loss in the active CTL cells due to encounters of tumor cells which is assumed to be proportional to the product of the densities of tumor cells and active CTL cells. Which is subtracted in the equation represents rate of change active CTL population $\frac{d{I}_{H}}{dt}$

• In the absence of active (hunting) CTL cells and chemotherapeutic drug both the tumor cell population and helper T-cell population are assumed to grow logistically.

• Mass-action kill rate assumes that all immune cells are similarly prone to communicate with any tumor cell: it assumes spatial homogeneity.

• *I _{H}* helper (resting) T-cells are not able to attack and destroy tumor cells directly but it convert CTLs into active (hunting) CTL helper either by releasing cytokines (interleukin-2) or by direct contact with them.

• *I _{R}* active (hunting) CTL cells have the role of attacking, destroying, or ingesting the tumor cells.

• Chemotherapeutic drug destroys tumor cells as well as helper T-cells and active CTL cells; that is, chemotherapeutic drug has a negative effect on both tumor cells and immune cells and this negative effect is expressed by subtracting this destroy from each cell population equation.

The model parameters are described as follows, for more details see [37] (Table 1).

##### Optimization strategies for OCP

The problem formulation which is given in (2.1),(2.2), based on open loop control (also referred to it in literatures as dynamic optimization ), i.e., feedback is not utilized (control input u(t) is independent of state). There are three main approaches to numerically solve continuous time OCP:

**1. Dynamic programming methods:** The optimal criterion in continuous time is based on the Hamilton-Jacobi-Bellman partial differential equation, for more details see [4], which is not in the scope of this paper.

**2. Indirect methods:** It take an approach optimize first then discretize, also relies on Pontryagin’s Maximum Principle (PMP), for more details see [32]. Typically, the optimal control problem is turned into TPBVP containing the same mathematical information as the original one by means of necessary conditions of optimality, for more details see [3], [38] and [39].

**3. Direct methods:** It take an approach discretize first then optimize, it can be applied without deriving the necessary condition of optimality. Direct methods are based on a finite dimensional parameterization of the infinite dimensional problem. The finite dimensional problem is typically solved using an optimization method, such as nonlinear programming (NLP) techniques. NLP problems can be solved to local optimality relying on the so called Karush-Kuhn-Tucker conditions (KKT), for more details see [3], if we are using KKT we can claim the first-order conditions of optimality. These conditions were first derived by Karush in 1939 [23], and later, in 1951, independently by Kuhn and Tucker [26].

##### Indirect methods

In indirect methods the necessary optimality conditions is derived by using Pontryagin’s maximum principle [32], by considering a simple optimal control problem.

$\mathrm{min}J\left(t,x\left(t\right),u\left(t\right)\right)=\Phi \left({t}_{f},x\left({t}_{f}\right)\right)+{\displaystyle \underset{{t}_{I}}{\overset{{t}_{F}}{\int}}[L\left(t,x\left(t\right),u\left(t\right)\right)}dt$ Subject to $\dot{x}(t)=f\left(\left(t,x\left(t\right),u\left(t\right)\right)\right)$ , dynamic constraints, ${\Psi}_{x}(tF)=0,$ final boundary condition,

where *J* is the objective function in Bolza form, where Bolza form is the sum of the Mayer term Φ(*t _{f}, x(t_{f}*)), and the Lagrange term, ${\int}_{tI}^{tF}[L(t,x(t),u(t))dt,$

the basic principle of Pontryagin’s maximum principle is defining the Hamiltonian, where Hamiltonian is a scalar function $\mathscr{H}:\left[{t}_{I},{t}_{F}\right]\times {\mathbb{R}}^{{n}_{x}}\times {\mathbb{R}}^{{n}_{u}}\times {\mathbb{R}}^{{n}_{x}}\to \mathbb{R},$ defined $H:[{t}_{I},{t}_{F}]\times {R}^{{n}_{x}}\times {R}^{{n}_{u}}\times {R}^{{n}_{x}}\to R$ defined as $H(t,x(t),u(t),\lambda (t))=L(t,x(t),u(t))+{\lambda}^{T}(t)f(t,x(t),u(t)),$ by defining the the auxiliary function $\Phi [{t}_{I},{t}_{F}]\times {R}^{{n}_{x}}\to {R}^{{n}_{x}}\to R,\text{by:}$

$\Phi (t,x)=\Phi (t,x(t))+{\mu}^{T}\Psi (x(t)),$

by setting first variation of the Lagrangian to zero where $\mathcal{L}:[{t}_{I},{t}_{F}]\times {\mathbb{R}}^{{n}_{x}}\times {\mathbb{R}}^{{n}_{u}}\times {\mathbb{R}}^{{n}_{q}}\times {\mathbb{R}}^{{n}_{x}}\to \mathbb{R}$ is defined as

$L\left(t,x(t),\mu ,\lambda (t)\right)=\left(\Phi \left({t}_{f},x({t}_{f}\right)\right)+{\mu}^{T}\Psi \left(x\left({t}_{F}\right)\right)+{\displaystyle \underset{{t}_{I}}{\overset{{t}_{F}}{\int}}\left[L\left(t,x\left(t\right),u\left(t\right)\right)+{\lambda}^{T}\left(t\right)\left(x-f\left(t,\dot{x}\left(t\right),u\left(t\right)\right)\right)\right]}dt\text{(3}\text{.2)}$

i.e., Integrating by parts the last term on the right side in Eq.(3.2), it yields:

$L(t,x(t),\mu ,\lambda (t))={[\phi (t,x(t))]}_{t={t}_{F}}-{\lambda}^{T}({t}_{F})x({t}_{F})+{\lambda}^{T}({t}_{I})x({t}_{I})+{\displaystyle {\int}_{{t}_{I}}^{{t}_{F}}[(H(t,x(t),u(t),\lambda (t)+{\dot{\lambda}}^{T}(t)x(t)]dt}$

We can conclude that the necessary optimality conditions for the unconstrained optimal is derived which stated as:

$\begin{array}{llll}\frac{d\lambda}{dt}\hfill & =\hfill & -\frac{\partial \mathscr{H}}{\partial x}\hfill & ,\text{adjointequations,}\hfill \\ \partial \mathscr{H}\hfill & =\hfill & 0\hfill & ,\text{controlequations,}\hfill \\ \lambda ({t}_{F})\hfill & =\hfill & {\left[\frac{\partial \phi}{\partial x}\right]}_{t={t}_{F}}\hfill & ,transversality\text{}conditions,\hfill \\ \lambda ({t}_{I})\hfill & =\hfill & 0\hfill & ,transversality\text{}conditions.\hfill \end{array}$

They are referred to as the Euler-Lagrange equations.

For more details see [11], which will be solved numerically. Many numerical methods which are based on the Euler- Lagrange differential equation (EL-DEQ) are available to solve the TPBVP. One may classify these numerical methods according to the particular approach used in [2].

In this article we introduce two approaches namely:

• Indirect Collocation Method.

• Forward-Backward Sweep Method (FBSM).

The purpose of examining two methods is to choose the method with the higher accuracy and set this method as the base method to compare against it. In the following we explained one of the most important methods for solving resultant system of the differential equations.

Indirect collocation method [36]: The Indirect Collocation Method is the merit of the techniques for solving TPBVP, finite difference method with continuous extension as well as collocation method are the mathematical tools for this method. As a difference methods it is based on Implicit Runge–Kutta (IRK) method which is equivalent to collocation method according to the following theorem.

Referring to the theorem of Guillou & Soule´1969, Wright 1970 which states that:

The collocation method defined as:

given s positive integer and *c _{1}, ..., c_{s}* distinct real numbers (typically between 0 and 1), the corresponding collocation polynomial

*u(x)*of degree s is defined by

$\begin{array}{lllll}\hfill & u({x}_{0})\hfill & =\hfill & {y}_{0}\hfill & (\text{initialvalue),}\hfill \\ \hfill & u\stackrel{\xb4}{}({x}_{0}+{c}_{i}h)\hfill & =\hfill & f({x}_{0}+c{}_{i}h,u\left({x}_{0}+{c}_{i}h\right),i=1,\mathrm{...},s.\hfill & \hfill \end{array}$

the numerical solution is given

$\begin{array}{ccc}{y}_{1}& =& u\left({x}_{0}+h\right).\end{array}$

is equivalent to the s-stage IRK-method which is defined as

$\begin{array}{l}{k}_{i}=f({x}_{0}+{c}_{i}h,{y}_{0}+h{\displaystyle \sum _{j=1}^{s}{a}_{ij}}{k}_{j})),\text{\hspace{1em}}i=1,\mathrm{...},s,\\ {y}_{1}={y}_{0}+h{\displaystyle \sum _{i=1}^{S}{b}_{i}}{k}_{i}.\end{array}$

with coefficients

${a}_{ij}={\displaystyle {\int}_{0}^{{c}_{i}}{\ell}_{j}}(t)dt,\text{\hspace{1em}\hspace{1em}}{b}_{j}={\displaystyle {\int}_{0}^{1}{\ell}_{j}}(t)dt,\text{\hspace{1em}\hspace{1em}}i,j=1,\mathrm{...},s,$

where the ${\ell}_{j}(t)$ are the Lagrange polynomials

${\ell}_{j}\left(t\right)=\prod _{k\ne j}\frac{\left(t-{c}_{k}\right)}{\left({c}_{j}-{c}_{k}\right)}.$

for theorem proof see [21].

**Forward-backward sweep method [27]:** In FBSM the initial value problem of the state equation is initialized by using an estimate for the control and costate variables and solved forward in time. Then the costate final value problem is solved backwards in time. An early reference to a technique that has the forward-backward flavor is [30]. FBSM method can be summarized by the following algorithm.

Information about convergence and stability of Runge-Kutta 4 ODE’s solver can be found in [20].

##### Direct methods

Sequential and simultaneous, are the two main classes of direct methods: for more details see [41]. Sequential methods only parameterize the control while simultaneous methods parameterize both the state and control.

The differences on how we discretize the problem and how the continuity between discritizated intervals are deﬁned result in different transcription method. In the next sections different methods depend on different discretization schemes and different continuity conditions are presented namely multiple shooting method, trapezoidal direct collocation and Hermite Simpson’s direct collocation.

**Multiple shooting method:** We start with the single shooting method, since single shooting is a special case of multiple shooting [12].

**Step1: Transcription:**

In single shooting system the time interval [*t _{0}, t_{f}*] is divided into equal sub intervals (segments) $[{t}_{i},{t}_{i+1}]$
such that $0={t}_{0}<{t}_{1}<\mathrm{...}<{t}_{N}={t}_{f},$
where

*N*is the total number of sub intervals.

Then the control vector is transformed into a parameterized finite dimensional control vector *u(t,q)* that depends on the finite dimensional parameter vector $q\in {R}^{N}$
there are several parametrization schemes for more details see [3,8], we assume piece wise constant control:

$u(t):={u}_{k}\text{\hspace{1em}}for\text{\hspace{1em}}t\in [{t}_{k},{t}_{k+1}),k=0,\mathrm{...},N-1.$

3. The initial value problem (IVP)

$\dot{x}\left(t\right)=f\left(x\left(t\right),u\left(t,q\right),t\right),\text{\hspace{1em}}\text{\hspace{0.22em}}x\left({t}_{0}\right)={x}_{0},\text{\hspace{1em}}\text{\hspace{0.22em}}t\in \left[{t}_{0},{t}_{f}\right],$

which is solved to yield the state vector $x\left(t,\text{}q\right)$
in the time interval [*t _{0}, t_{f}*].

4. The integral of the objective function is calculated together with the initial value problem solution (by using a quadrature formula). The Lagrangian part of the cost function is evaluated on each interval independently.

##### Optimization

1. Due to the numerical simulation the model equations are eliminated. The path constraints are also discretized. Thus the optimal control problem Eq.(3.1) is rewritten as:

$\begin{array}{lllll}min\text{}J\hfill & =\hfill & \phi (x({t}_{0}),q,{t}_{0},x({t}_{f}),q,{t}_{f})\hfill & +\hfill & {\displaystyle {\int}_{{t}_{0}}^{{t}_{f}}g}(x(t,q),u(t,q))dt,\hfill \\ subject\text{to}\hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & C(x({t}_{i},q),u({t}_{i},q))\hfill & \le \hfill & 0,\hfill \\ \hfill & \hfill & E(x({t}_{0}),q,{t}_{0},x({t}_{f}),q,{t}_{f})\hfill & =\hfill & 0.\hfill \end{array}$

Which is an NLP problem that is solved using the Interior Point Method (IP) in this paper.

**Continuity constraints:** Since the simulation is done over the whole time horizon this method does not have continuity constraints.

**Accuracy:** Determine the accuracy of the finite-dimensional approximation and if necessary repeat the transcription and optimization steps (i.e. go to step 1).

single shooting method has the following drawbacks:

• Convergence of the NLP solution is slow, because of the high nonlinear dependence of the objective and constraint functions on the variable *u*,

• NLP solver cannot initialize with an initial guess for *x ^{1}, ..., x^{N}* even it is available.

• Parallel evaluation of the states and objective functions is impossible because of the recursive elimination for both states and objective functions.

**Direct multiple shooting method [10]:** This method combines the advantages of simultaneous methods like collocation method with the main advantages of the single shooting method, so that it is sometimes called a hybrid method [17].

We follow the same steps in single shooting method except that we solve ODE’s in each interval [*t ^{i}, t^{i+1}*] independently, starting with an artificial initial value

*S*:

^{j}Solving these initial value problems numerically, we can obtain trajectory pieces ${x}_{i}\left(t;{s}_{i};{q}_{i}\right),$ where the extra arguments after the semicolon are introduced to denote the dependence on the interval’s initial values and controls. Simultaneously with the decoupled ODE’s solution, we also numerically compute the integrals is numerically computed by:

${l}_{i}({s}_{i},{q}_{i}):={\displaystyle \underset{{t}_{i}}{\overset{{t}_{i+1}}{\int}}L({x}_{i}({t}_{i};{s}_{i},{q}_{i}),{q}_{i})dt}.$

**Continuity constrains:** To enforce continuity between discretized intervals, the following NLP constraints is added at the interface of each sub-interval, this constraint is formed such that the propagated (as an illustration for propagation means consider a cannon be aimed such that the cannonball hit its target i.e., shoot) or integrated value of the state from the previous phase match the value of the state at the current state. The continuity depends on propagation which is approximate because propagation is using algebraic formulas based on numerical integration schemes (or discretization schemes).

**Figure 3.1:** Collocation Defect Constraints.

So Continuity Constrains is defined as:

${s}_{i+1}={x}_{i}\left({t}_{i+1};{s}_{i},{q}_{i}\right),$

We have the following NLP, but contains the extra variables *s ^{i}*, and has a block sparse structure.

$\begin{array}{llll}\mathrm{min}\hfill & J=\phi (x({t}_{0}),{q}_{0},{s}_{0},{t}_{0},x({t}_{f}),q,{t}_{f})\hfill & +\hfill & {\displaystyle \sum _{i=0}^{N-1}g}(x(t,{q}_{i},{s}_{i}),u(t,{q}_{i},{s}_{i})),\hfill \\ \hfill & \text{subjectto}\hfill & \hfill & \hfill \\ \hfill & {s}_{0}-{x}_{0}\hfill & =\hfill & 0,(\text{initialvalue),}\hfill \\ \hfill & {s}_{i+1}-{x}_{i}({t}_{i+1};{s}_{i},{q}_{i})\hfill & =\hfill & 0,(continuity),\hfill \\ \hfill & C(x({t}_{i},{q}_{i},{s}_{i}),u({t}_{i},{q}_{i},{s}_{i}))\hfill & \le \hfill & 0,(\text{discretized}\text{path}\text{constraints)},\hfill \\ \hfill & E(x({t}_{0}),{q}_{0},{s}_{0},{t}_{0},x({t}_{f}),{q}_{N},{s}_{N},{t}_{f})\hfill & =\hfill & 0,((\text{terminal}\text{constraints)}.\hfill \end{array}$

**Direct collocation:** Problem size, non-linearity and sparsity of the NLP resulting from direct transcription methods are the prices of moving from one method to another. In single shooting the NLP is highly nonlinear while its size is small. Multiple shooting is less nonlinear but larger and sparser. The direct collocation goes further in the same direction as it is less nonlinear, sparser but even larger, direct collocation method is introduced by Dickmanns [16] as a direct transcription for solving optimal control problems. It has the fundamental steps of direct transcription method, but differs slightly.

It similarly discretizes the state trajectory $\left[{t}_{0},{t}_{f}\right]$
into equal sub intervals (segments) $\left[{t}_{i},{t}_{i+1}\right]$
, such that $0={t}_{0}<{t}_{1}<\mathrm{...}<{t}_{N}={t}_{f}$
where *N* is the total number of sub intervals.

But it differs at

It further divide the interval $\left[{t}_{i},{t}_{i+1}\right]$ into K sub intervals $\left[{\tau}_{j},{\tau}_{j+1}\right],$ where

${\tau}_{j}={t}_{i}+{h}_{i}{\alpha}_{j},\text{\hspace{1em}}j=1,\mathrm{...},K,\text{\hspace{1em}}{h}_{i}={t}_{i+1}-{t}_{i},\text{(3}\text{.3)}$

It employs an interpolating function to approximate the state of the system. Usually the interpolating function is poly nomial.

Polynomials consider the following *K ^{th}*-degree piece wise polynomial:

$X(t)\approx {\displaystyle \sum _{k=0}^{K}{a}_{i}}{(t-{t}_{i})}^{k},\text{\hspace{0.17em}}t\in [{t}_{i},{t}_{i+1}],\text{(3}\text{.4)}$

suppose further that the coefficients (a^{0},..., a^{K}) of the piece wise polynomial are chosen to match the value of the function at the beginning of the step, i.e.,

$X\left({t}_{i}\right)=x\left({t}_{i}\right).$

finally, suppose we choose to match the derivative of the state at the points defined by Eq.(3.3), i.e.,

$\dot{X}({\tau}_{j})=f(x({\tau}_{j}),{\tau}_{j}),\text{\hspace{1em}\hspace{1em}}j=1,\mathrm{...},K.\text{(3}\text{.5)}$

Eq.(3.5) is called collocation condition because the approximation to the derivative is set equal to the right-hand side of the differential equation evaluated at each of the intermediate points $\left({\tau}_{1},\mathrm{...},{\tau}_{K}\right)$ .

– By setting *K = 2* in Eq.(3.4) we get the quadratic interpolation polynomial which results in trapezoidal method.

– By setting *K = 3* in Eq.(3.4) results in Hermite–Simpson’s method.

Trapezoidal method and Hermite–Simpson’s method are the two direct collocation methods we consider in this article.

**Defects constraints:** The difference between the first derivative of the interpolating polynomial at the midpoint of a segment and the first derivative calculated from the equations of motion at the segment midpoint is used as the defect. To have a good approximation of the actual states the defects must approach zero and thus interpolating polynomial is a good approximation, the defect is written in the following form:

${\zeta}_{j}=\dot{X}({\tau}_{j})-f(x({\tau}_{j}),{\tau}_{j}),\text{(3}\text{.6)}$

the defect constraint for the trapezoidal method (by setting $P\left(x\right)=x\left({\tau}_{j}\right)-\frac{{h}_{j}}{2}\left({f}_{j}+{f}_{j+1}\right)$ and ${x}_{i+1}=x\left({\tau}_{j+1}\right)$ at Eq.(3.10))

is defined as:

${\zeta}_{j}\left({\tau}_{j}\right)=x\left({\tau}_{j+1}\right)-x\left({\tau}_{j}\right)-\frac{{h}_{j}}{2}\left({f}_{j}+{f}_{j+1}\right)=0,$

the defect constraint for the Hermite–Simpson’s method (by setting $P(x)=x({\tau}_{j})-\frac{{h}_{k}}{6}({f}_{j}+4{\overline{f}}_{j}+{f}_{j+1})$ and $x{\text{}}_{i+1}\text{}=x\left(\text{}{\tau}_{i+1}\text{}\right)\text{})$ at Eq.(3.6) is defined as:

let ${\overline{\tau}}_{j}=\left({\tau}_{j}+{\tau}_{j+1}\right)/2,\text{\hspace{0.17em}}and\overline{\text{\hspace{0.17em}}u}({\overline{\tau}}_{j})=\left(u\left({\tau}_{j}\right)+u\left({\tau}_{+1}\right)\right)/2,$ $\begin{array}{llll}\text{x \u02dc}(\overline{{\tau}_{\text{j}}})\hfill & =\hfill & \frac{1}{2}[x({\tau}_{j})+x({\tau}_{j+1})]\hfill & +\frac{{h}_{k}}{8}({f}_{j}-{f}_{j+1}),\hfill \\ {\overline{f}}_{j}\hfill & =\hfill & f[\text{x \u02dc}({\overline{\tau}}_{\text{j}}),\text{u \xaf (}{\tau}_{\text{j}}),\text{p},{\overline{\tau}}_{\text{j}}]\hfill & ,\hfill \\ \zeta ({t}_{j})\hfill & =\hfill & x({\tau}_{j+1})-x({\tau}_{j})\hfill & +\frac{{h}_{k}}{6}({f}_{j}+4{\overline{f}}_{j}+{f}_{j+1})=0.\hfill \end{array}$

• The quadrature rule used for integration is consistent with the numerical method used for solving the differential equa- tion. If one is using a Runge-Kutta method for solving the differential equation, the cost would also be approximated using Runge-Kutta integration. In the case of an Hermite–Simpson’s collocation method, the integration rule is Her- mite–Simpson collocated quadrature rule [33].

• Thus the optimal control problem Eq.(3.1) is rewritten as: let *s _{i}* be values of state vector at at grid point,

*sˆ*which represent the states at the collocation points in each subinterval.

_{i}$\begin{array}{llll}\mathrm{min}\hfill & J=\phi (x({t}_{0}),{q}_{0},{s}_{0},{t}_{0},x({t}_{f}),q,{t}_{f})\hfill & +\hfill & {\displaystyle \sum _{i=0}^{N-1}g}(x(t,{q}_{i},{s}_{i},{\widehat{s}}_{i}),u(t,{q}_{i},{s}_{i},{\widehat{s}}_{i})),\hfill \\ \hfill & \text{subjectto}\hfill & \hfill & \hfill \\ \hfill & {s}_{0}-{x}_{0}\hfill & =0,\hfill & (\text{initialvalue),}\hfill \\ \hfill & {s}_{i+1}-{x}_{i}({t}_{i+1};{s}_{i},{\widehat{s}}_{i},{q}_{i})\hfill & =0,\hfill & (\text{continuity}),\hfill \\ \hfill & C(x({t}_{i},{q}_{i},{s}_{i},{\widehat{s}}_{i},),u({t}_{i},{q}_{i},{s}_{i},{\widehat{s}}_{i}))\hfill & \le 0,\hfill & (\text{discretizedpath}\text{constraints)},\hfill \\ \hfill & E(x({t}_{0}),{q}_{0},{s}_{0},{t}_{0},x({t}_{f}),{q}_{N},{s}_{N},{t}_{f})\hfill & =0.\hfill & \text{(terminal}\text{constraints)}.\hfill \end{array}$

Solution of the NLP problem

Numerical methods for solving NLPs fall into categories: gradient based (local) methods and heuristic (global) methods, gradient based (local) approach is now described in the following algorithm (the heuristic method is out of the interest of this article, for more details see [40]).

The convergence may depend on a given tolerance or depends on no further change in the objective function after several iterations. There are many ways to modify direction selection and step size, and this leads to many different algorithms. Some of common ones included Steepest Descent Methods, Conjugate Direction Method, Simplex Method, Interior Point Method (IP) and Sequential Quadratic Programming, more details on the IP methods will be given here.

**IP Method:** it is typically eliminate constraints of the form $\underset{\xaf}{h}\le h(x)\le \overline{h}$
by introducing slack variables additional equality constraints. After locating an interior point, i.e., a point where $\underset{\xaf}{x}<x<\overline{x}$
holds with strict inequality, which may require
reformulating some decision variables as parameters, they formulate a barrier problem of the form:

$\begin{array}{llllll}\mathrm{min}\hfill & f\left(x\right)\hfill & -\hfill & \mu {\displaystyle \sum _{i=1}^{{n}_{x}}l}og({\overline{x}}_{i}-{x}_{i})\hfill & +\hfill & log({x}_{i}-{\underset{\xaf}{x}}_{i})),\hfill \\ \hfill & \text{subjectto}\hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill & g\left(x\right)=0.\hfill & \hfill & \hfill \end{array}$

Sequences of barrier problems for increasing values of the barrier parameter µ are then solved with Newton-type methods, for more details see [7,40].

##### Accuracy

Since the solution obtained from the NLP problem is a discrete set of numbers we don’t know what happen between the discrete points, a continuous approximation to the discrete solution resulted from the NLP is needed, i.e., representing the solution (interpolation) [6].

**Spline representation:** To get the solution as continuous approximation $\left[\tilde{\text{y}}(t),\tilde{\text{u}}(t)\right]$
from the NLP solution we use spline representation (where spline is a sequence i.e., whole collection of polynomial) defined as follows:

$\left[\begin{array}{l}y(t)\\ u(t)\end{array}\right]=\left[\begin{array}{l}\tilde{y}(t)\\ \tilde{u}(t)\end{array}\right]={\displaystyle \sum _{i=1}^{{n}_{1}}{\alpha}_{i}}{\beta}_{i}(t),\text{(3}\text{.7)}$

where *n*^{1} = 2*M*, *M* is the number of mesh points, *y(t), u(t)* are the true state and control respectively, $\left[\tilde{\text{y}}(t),\tilde{\text{u}}(t)\right]$
are the approximate state and control respectively, the functions *β ^{i}(t)* form a basis for C

_{0}or C

_{1}cubic B-splines with

*n*

^{1}= 2

*M*, where

*M*is the number of mesh points. The coefficients αi in the state or control variable representation which is defined by different Interpolation of discrete solution depending on the discritization method, the spline approximation Eq. (3.7) must match the state at the grid points i.e.,

$\tilde{y}\left(t\right)={y}_{k},\text{\hspace{1em}\hspace{1em}}\text{\hspace{0.17em}}\text{\hspace{0.22em}}for\text{}\text{}k=1,\mathrm{...},M,$

the derivative of the spline approximation must match the right-hand side of the differential equations,

$\frac{d}{dt}\tilde{y}\left(t\right)={f}_{k},\text{\hspace{1em}\hspace{1em}}\text{\hspace{0.17em}}\text{\hspace{0.22em}}for\text{\hspace{0.17em}}k=1,\mathrm{...},M,$

we require the spline approximation in Eq.(3.7) to match the control at the grid points, i.e.,

$\tilde{u}\left({t}_{k}\right)={u}_{k},,\text{\hspace{1em}\hspace{1em}}\text{\hspace{0.17em}}\text{\hspace{0.22em}}for\text{}\text{}k=1,\mathrm{...},M,$

the state variable is C^{1} cubic function whereas the control variable is C^{0} linear or quadratic functions depending on the discritization method. Since both $\tilde{y}\left(t\right),\tilde{u}\left(t\right)$
are approximate to the true solution the degree to which this approximation approximate the true solution need to be determined which is known as discretization error and relative local error.

**Discretization error:** by assuming that the computed control is correct and assume the spline solutions $\tilde{y}\left(t\right),\tilde{u}\left(t\right)$
produced from the NLP, and the single interval ${t}_{k}\le t\le {t}_{k}+{h}_{k},$
from the state equation we can define

$y\left({t}_{k}+{h}_{k}\right)=y\left({t}_{k}\right)+{\int}_{{t}_{k}}^{{t}_{k}+{h}_{k}}\dot{y}\text{d}t=y\left({t}_{k}\right)+{\int}_{{t}_{k}}^{{t}_{k}+{h}_{k}}f\left(y,u,t\right)\text{\hspace{0.17em}}\text{d}t,$

where *y* and *u* are the true state and control values, we can consider the approximation:

$\widehat{y}({t}_{k}+{h}_{k})=y({t}_{k})+{\displaystyle {\int}_{{t}_{k}}^{{t}_{k}+{h}_{k}}f}[\tilde{y}(t),\tilde{u}(t),t]dt,\text{(3}\text{.8)}$

from Eq.(3.8) we can define the discretization error on the *K ^{th}* mesh iteration as:

${\eta}_{k}=\underset{i}{\mathrm{max}}\left\{{a}_{i}\left|{\tilde{y}}_{i}\left({t}_{k}+{h}_{k}\right)-{\stackrel{\u2322}{y}}_{i}\left({t}_{k}+{h}_{k}\right)\right|\right\},$

For $i=1,2,\mathrm{...},n,$
where the weights *a ^{i}* are chosen to appropriately normalize the error.

Relative local error: is the maximum relative error over all components *i* in the state equations *y˙_ f* evaluated in the interval k, and is defined as

${\u03f5}_{k}\approx \underset{i}{\mathrm{max}}\frac{{\eta}_{i,k}}{\left({w}_{i}+1\right)},$

where the scale weight ${w}_{i}=\underset{k=1}{\overset{M}{\mathrm{max}}}\left[\left|{\tilde{y}}_{i,k}|,|{\dot{\tilde{y}}}_{i,k}\right|\right],$
defines the maximum value for the *i*^{th} state variable or its derivative over the *M* grid points in the phase. An equivalent form can be defined as follows the absolute local error on a particular step by

${\eta}_{i,k}={\displaystyle {\int}_{{t}_{k}}^{{t}_{k+1}}\left|{\epsilon}_{i}(s)\right|}ds,\text{(3}\text{.9)}$

Where $\epsilon \left(t\right)=\dot{\tilde{y}}\left(t\right)-f\left[\tilde{y}\left(t\right),\tilde{u}\left(t\right),t\right]$ defines the error in the differential equation as a function of t. An accurate estimate for the integral in Eq.(3.9) is evaluated by using a standard quadrature method, (i.e., Romberg quadrature algorithm).

In this part we compare the results from the our base line article [37], with another indirect method solved by sweep method [27] and another three direct method namely direct multiple shooting, trapezoidal direct collocation and Hermite–Simpsons’s direct collocation.

**Figure 4.1:** Tumor Cell population using 3D plot (left figure) and 2D plot (right figure).

**Figure 4.2:** Helper T cells using 3D plot (left figure) and 2D plot (right figure).

**Figure 4.3:** Active CTL using 3D plot (left figure) and 2D plot (right figure).

**Figure 4.4:**Density of Chemo Therapy Drug (3D plot) left side, (2D plot) right side.

**Figure 4.5:** Control Plot using different methods.

**Figure 4.6:** Objective Function Using different methods.

**Figure 4.7:** State Error for Multiple Shooting Method.

**Figure 4.8:** State Error for Hermite Simpson Method.

**Figure 4.9:** State Error for Direct Multiple Shooting Method.

**Figure 4.10:** dx/dt - f(t,x,u) for Multiple Shooting.

**Figure 4.11:** dx/dt - f(t,x,u) for Hermite Simpson.

**Figure 4.12:** dx/dt - f(t,x,u) for trapezoid.

In all the previous figures the direct methods achieve better results than indirect methods. The values of controls and number of tumor cells made the big differences in the objective function of the theses methods. The number of tumor cells is lower in direct method than indirect methods also direct methods achieve better control than indirect method the following table shows the value of objective function for five different methods.

Since the direct methods achieve better results than indirect method this motivate us to give more insight into direct methods through new merits such as discretization error, absolute local error, number of NLP iterations, NLP time the following table gives the values of the previous merits for the methods of direct methods only.

Within the direct method apart from indirect methods the results are butter in both Hermite Simpson DT and direct multiple shooting is butter than the trapozidal since both two methods represented by polynomials, the Hermite Simpson’s DT is better than the Trapezoidal DT since Hermite Simpson’s DT uses a higher order polynomial than trapezoidal DT which is an agreement with the concept of (p-method), for more details see [29], the hermite Simpson is the most accurate method on the price of time.

The following plots shows the different discretization errors and state errors for the three method.

In this article five numerical techniques to study the tumor immune dynamic optimization model are discussed, these methods are indirect method in which the resulting TPBVP is studied by collocation method, and by the backward forward sweep method and direct methods which are multiple shooting method, trapezoid direct collocation method, Hermite-Simpson’s direct collocation. Many studies claim that indirect methods have drawbacks. The two major drawbacks are that the differential equations obtained are often difficult to solve due to strong nonliterary and instability and also define a guess initial solution for the Lagrangian multiplier. As a matter of fact, these variables do not have a physical meaning, this leads to problems in the definition of an initial guess solution to start the algorithm. Another problem that it is necessary to take into account in the resolution of a BVP is that the existence and uniqueness of solution is not guaranteed as in an initial value problem. Each problem may have a unique solution, several solutions or no solution at all. By comparing different methods both in direct and indirect methods we can claim that direct method can be used to get better results than indirect method and it is easy to implement. A question that is needed to be answered is how to choose the discretization method among different methods related to direct method to get the best results this question can be a future work.

- Gibbs WW. Untangling the roots of cancer. Scientific America. 2003.
- Abbas AK, Lichtman AH, Pillai S. Cellular and molecular immunology. Saunders Elsevier. 2014.
- Curiel T. Tregs and rethinking cancer immunotherapy. Journal of Clinical Investigation. 2007; 117: 1167-1174. PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17476346
- Kirschner D, P Panetta. Modeling immuno therapy of the tumor–immune interaction. Journal of Mathematical Biology. 1998; 37: 235-252. PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9785481
- Kirschner DE, TL Jackson, JC Arciero. A mathematical model of tumorimmune evasion and siRNA treatment. Discrete and continous dynamical systems series- B. 2003; 37: 39-58.
- K Leon, K Garcia, J Carneiro. A Lage. How regulatory CD25(+)CD4(+) T cells impinge on tumor immunobiology? On the existence of two alternative dynamical classes of tumors. Journal of Theoretical Biology. 2007; 247: 122-137.
- De Pillis LG, Radunskaya AE. Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations. Journal of Theoretical Biology. 2006; 238.
- Schattlera H, Urszula L. Optimal Control for Mathematical Models of Cancer Therapies. Springer Publishing Co., USA. 2015.
- Sharma S, Samanta GP. Dynamical Behaviour of a Tumor-Immune System with Chemotherapy and Optimal Control. Journal of Nonlinear Dynamics. 2013: 1-13, 2013.
- Sweilam NH, Al-Ajami TM. Legendre spectral-collocation method for solving some types of fractional optimal control problems. Journal of Advanced Research, 2015; 393-403. PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26257937
- García-Heras J, Soler M, Sáez FJ. A Comparison of Optimal Control Methods for Minimum Fuel Cruise at Constant Altitude and Course with Fixed Arrival Time. Procedia Engineering. 2014; 80:231-244.
- Rao AV, Benson DA, Darby C, Patterson MA, Francolin C, et al. Algorithm 902: GPOPS, A MATLAB software for solving multiple-phase optimal control problems using the gauss pseudospectral method. ACM Transactions on Mathematical Software (TOMS), 2010; 37: 1-39.
- Biral F, Bertolazzi E, Bosetti P. Notes on Numerical Methods for Solving Optimal Control Problems. IEEJ Journal of Industry Applications. 2015; 5:154-166
- Betts JT. A Survey of Numerical Methods for Trajectory Optimization. Control and Dynamics. 1998; 21:193-207.
- Rao AV. A survey of numerical methods for optimal control. Advances in the Astronautical Sciences. 2009; 135: 497-528.
- Joshi HR. Optimal control of an HIV immunology model. Optimal Control Applications and Methods. 2002; 23: 199-213.
- Zaman G, Han Kang Y, Jung IH. Stability analysis and optimal vaccination of an SIR epidemic model. BioSystems. 93: 240-249. 2008. PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18584947
- Pillis LG, Radunskaya AE. A mathematical model of immune response to tumor invasion. MIT. 2003; 1661-16668.
- De Pillis LG, W Gu, Fister KR, Head T, Maples K, et al. Chemotherapy for tumors: An analysis of the dynamics and a study of quadratic and linear optimal controls. Mathematical Biosciences. 2007; 209: 292-315.
- Bellman RE. Dynamic Programming. Courier Corporation. 2003.
- Pontryagin LS. Mathematical Theory of Optimal Processes. CRC Press. March 1987.
- Anita S, Arnautu V, Capasso V. An introduction to optimal control problems in life sciences and economics: from mathematical models to numerical simulation with MATLAB®. Modeling and simulation in science, engineering and technology. Birkhäuser. New York. 2011.
- Sweilam NH, AL-Mekhla M. On the Optimal Control for Fractional Multi-Strain TB Model. Optimal Control, Applications and Methods. 2016.
- Karush W. Minima of Functions of Several Variables with Inequalities as Side Constraints. Ph.D. Department of Mathematics. University of Chicago. Chicago. 1939
- H Kuhn, A Tucker. Nonlinear Programming. 1951; 481-492, California. University of California Press. Berkeley.
- Bryson AE, Ho YC. Applied optimal control. Hemisphere Publication Corporation. 1975.
- Aktas Z, Stetter HJ. A classification and survey of numerical methods for boundary value problems in ordinary differential equations. International journal for numerical methods in engineering. 1977; 11: 771-796.
- Shampine LF, Gladwell I, Thompson S. Solving ODEs with MATLAB. Cambridge University Press. 2003.20.
- Lenhart S and Workman JT. Optimal Control Applied to Biological Models. Chapman & Hall/CRC Mathematical and Computational Biology. CRC Press. Taylor & Francis Group. 2007.
- Mitter SK. The successive approximation method for the solution of optimal control problems. Automotica. 1996; 3:135-149.
- Hackbusch W. A numerical method for solving parabolic equations with opposite orientations. Computing. 1978; 20: 229-240.
- Victor VM. Practical Direct Collocation Methods for Computational Optimal Control. In Modeling and Optimization in Space Engineering. Volume 73 of Springer Optimization and Its Applications. Springer New York. 2013; 33-60.
- Chachuat B. Nonlinear and Dynamic Optimization: From Theory to Practice - IC-32: Spring Term. EPFL. 2009.
- Binder T, Blank L, Bock HG, Bulirsch R, Dahmen W, et al. Introduction to Model Based Optimization of Chemical Processes on Moving Horizons. In Introduction to Model Based Optimization of Chemical Processes on Moving Horizons. Springer Berlin Heidelberg. 2001; 295-339.
- Bock H, Plitt K. A multiple shooting algorithm for direct solution of optimal control problems. In 9th IFAC. Pergamon Press. 1984; 242-247.
- Diehl M, Findeisen R, Schwarzkopf S, Uslu I, Allgöwer F, et al. An Efficient Algorithm for Nonlinear Model Predictive Control of Large-Scale Systems Part I: Description of the Method. At-Automatisierungstechnik Methoden und Anwendungen der Steuerungs-, Regelungs-und Informationstechnik, 2002; 50: 557.
- Dickmanns ED, Well KH. Approximate solution of optimal control problems using third order hermite polynomial functions. In Marchuk GI, editor. Optimization Techniques IFIP Technical Conference Novosibirsk, number 27 in Lecture Notes in Computer Science, pages. Springer Berlin Heidelberg. 1974; 158-166.
- Törn A, Žilinskas A, Goos G, Hartmanis J, Barstow D, et al. Global Optimization, volume 350 of Lecture Notes in Computer Science. Springer Berlin Heidelberg. Berlin. Heidelberg. 1989.
- Biegler LT. Nonlinear programming: concepts, algorithms, and applications to chemical processes. Number 10 in MOS-SIAM series on optimization. SIAM. Philadelphia. 2010.
- Betts JT. Practical methods for optimal control and estimation using nonlinear programming. Advances in design and control. Society for Industrial and Applied Mathematics. Philadelphia. 2nd edition. 2010.
- Matthew PK. Transcription Methods for Trajectory Optimization A beginners tutorial. Technical report. Cornell University. 2015.
- E Hairer, Norsett SP, Wanner G. Solving Ordinary Differential Equations I Nonstiff Problems. Springer-Verlag Berlin Heidelberg, 2008.