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 specified 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 })=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₁,Bsub>1 are positive constants representing the weights of the terms. The first 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)$ reflects 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 justifies the quadratic terms in the functional. So the functional given in Eq. (2.1) should be minimized.

The dynamical system is defined 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}T\left(0\right)\ge T_0,I_H\left(0\right)\ge 0,I_R\left(0\right)\ge 0,D\left(0\right)\ge 0.$

Where

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

• II_H (t) is the active CTL cells (hunting CTL cells).

• II_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.

$\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, ${\displaystyle {\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 {ℝ}^{n_x}\times {ℝ}^{n_u}\times {ℝ}^{n_x}\to ℝ,$ 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 $ℒ:[t_I,t_F]\times {ℝ}^{n_x}\times {ℝ}^{n_u}\times {ℝ}^{n_q}\times {ℝ}^{n_x}\to ℝ$ 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{adjoint equations,}\hfill \\ \partial \mathscr{H}\hfill & =\hfill & 0\hfill & ,\text{control equations,}\hfill \\ \lambda (t_F)\hfill & =\hfill & {\left[\frac{\partial \phi }{\partial x}\right]}_{t=t_F}\hfill & ,transversalityconditions,\hfill \\ \lambda (t_I)\hfill & =\hfill & 0\hfill & ,transversalityconditions.\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₁, ..., 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{initial value),}\hfill \\ \hfill & u\stackrel{´}{}(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)),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,b_j={\displaystyle {\int }_0^1{\ell }_j}(t)dt,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 defined 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₀, 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}fort\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),x\left(t_0\right)=x_0,t\in \left[t_0,t_f\right],$

which is solved to yield the state vector $x\left(t,q\right)$ in the time interval [t₀, 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}minJ\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¹, ..., 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ⁱ, tⁱ⁺¹] 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).

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ⁱ, and has a block sparse structure.

$\begin{array}{llll}\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{subject to}\hfill & \hfill & \hfill \\ \hfill & s_0-x_0\hfill & =\hfill & 0,(\text{initial value),}\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,j=1,\mathrm{...},K,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,t\in [t_i,t_{i+1}],\text{ (3}\text{.4)}$

suppose further that the coefficients (a⁰,..., 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),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{}_{i+1}=x\left({\tau }_{i+1}\right))$ at Eq.(3.6) is defined as:

let ${\overline{\tau }}_j=\left({\tau }_j+{\tau }_{j+1}\right)/2,and\overline{u}({\overline{\tau }}_j)=\left(u\left({\tau }_j\right)+u\left({\tau }_{+1}\right)\right)/2,$ $\begin{array}{llll}\text{<mover accent="true"> x <mo>˜</mo> </mover>}(\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{<mover accent="true"> x <mo>˜</mo> </mover>}({\overline{\tau }}_{\text{j}}),\text{<mover accent="true"> u <mo>¯</mo> </mover> (}{\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ˆ_i which represent the states at the collocation points in each subinterval.

$\begin{array}{llll}\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{subject to}\hfill & \hfill & \hfill \\ \hfill & s_0-x_0\hfill & =0,\hfill & (\text{initial value),}\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{discretized 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 & =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{¯}{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{¯}{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}\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{¯}{x}}_i)),\hfill \\ \hfill & \text{subject to}\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¹ = 2M, 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 βⁱ(t) form a basis for C₀ or C₁ cubic B-splines with n¹ = 2M, 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,for\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,fork=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,,for\text{}k=1,\mathrm{...},M,$

the state variable is C¹ cubic function whereas the control variable is C⁰ 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{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}{\max }\left\{a_i\left|{\tilde{y}}_i\left(t_k+h_k\right)-{\stackrel{⌢}{y}}_i\left(t_k+h_k\right)\right|\right\},$

For $i=1,2,\mathrm{...},n,$ where the weights aⁱ 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

${ϵ}_k\approx \underset{i}{\max }\frac{{\eta }_{i,k}}{\left(w_i+1\right)},$

where the scale weight $w_i=\underset{k=1}{\overset{M}{\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).