Optimal Control of Anti-Angiogenesis and Radiation Treatments for Cancerous Tumor: Hybrid Indirect Solver

. Tis paper proposes a real-life volume reduction for cancer cells using optimal doses of radiation and an anti-angiogenic drug. A generalized dynamical system based on the difusion-consumption equation along with stimulation and inhibition factors is proposed. To achieve continuous and low-dose therapy, the related problem is simulated by an optimal regulator problem mathematically. By combining steepest descent, conjugate gradient, and Armijo techniques, a novel hybrid indirect iterative solver is designed. For accuracy and execution speed, the current solver is compared with an interior-point optimizer and sequential quadratic Hamiltonian methods. Cancer therapy under two diferent treatment strategies and 24 various versions of the general dynamical system is considered numerically. A comprehensive analysis of the corresponding outcomes is presented. Numerical results and related diagrams are provided.


Introduction
According to World Health Organization reports for the year 2019, cancer diseases are considered the main cause of death and the most important cause of life expectancy reduction [1].Malek and Abbasi [2][3][4][5] in years 2016, 2019, 2020, and 2022 introduced some mathematical models to simulate and solve optimal control problems for cancer therapy.Cancer cells form and reproduce abnormal cells irregularly, while the proliferation of adult normal cells is based on the division and replacement of dying cells.General cancer therapies have two categories: treatments that directly kill cancer cells and treatments that prevent cancer cell growth.Cancer therapy usually has negative side efects.Tus, treatment is acceptable if there is a reasonable balance between the outcome and its negative side efects.Anti-angiogenic treatment is a cancer therapy to prevent cancer cell proliferation.Radiation therapy leads to direct cancer cell killing and is frequently combined with antiangiogenic or other therapies to increase its efectiveness [6].Several mathematical models have been proposed to describe the evolution of tumor anti-angiogenesis as a dynamical system.Among them, the models developed by Hahnfeldt et al. [7] and Ergun et al. [8] are widely recognized as the most prominent.Hahnfeldt formulated a difusionconsumption equation to model the concentration of stimulators and inhibitors both inside and outside the tumor.His research revealed that the inhibitor would have a signifcant impact on the target endothelial cells in the tumor, ultimately leading to growth proportional to the tumor surface, while the efect of stimulators would be relatively independent of tumor or vascular size.Endothelial cells play a signifcant role in carrying blood, nutrients, and oxygen to vital organs, while also facilitating the removal of deoxygenated blood from these organs [9].In the Ergun model, the inhibition term is proportional to the tumor radius rather than its surface area, under the assumption that the system is in a steady state.Kienle et al. [10] considered an inconsistent convex combination of the Hahnfeldt (tumor surface) and Ergun (tumor radius) models; this combination would reduce the medical relevance.Tere are many studies that have separately used diferent dynamics of vascular carrying capacity.For example, Owolabi et al. [11] separately employed three diferent carrying capacity dynamics to investigate the singular arc for the optimal problems.In this article, we generalize the Hahnfeldt model by aiming to provide an improved description of the vasculature's carrying capacity while incorporating optimal strategies in the existence of anti-angiogenesis and radiotherapy.Here, we face both mathematical and clinical challenges; regarding the clinical aspect, treatment strategies for tumor reduction are encountered during treatment and fnal time, where the fnal time is fxed.From a mathematical point of view, we consider the quadratic term to control doses in the objective function.
Numerical methods for solving optimal control problems fall into two main categories: indirect and direct approaches.In the indirect method, the calculus of variations determines the frst-order optimality conditions for the optimal control problem [12][13][14].Tis approach yields a multipoint boundary-value problem to fnd optimal solutions.In direct methods, the state and control variables of the optimal control problem are approximated in an appropriate manner, transforming the problem into nonlinear programming (NLP) [15,16].Tese two approaches stem from diferent philosophies.Te indirect approach indirectly solves the problem (reason for naming indirect) by converting it into a boundary-value problem, solving a set of diferential equations for the optimal solution.In contrast, the direct method transcribes an infnite-dimensional optimization problem into a fnite-dimensional one, leading to the optimal solution.Iterative indirect methods, such as steepest descent, consist of two components.One part focuses on determining the direction of movement, while another part focuses on determining the length of the movement step required to reach the optimal solution [17,18].Hence, employing a hybrid indirect solver concept may be an appropriate approach.
Indirect methods utilize Pontryagin's maximum principle (PMP) to provide analytical insights into the optimal control problem, leading to a better understanding of the system's behavior.Moreover, indirect methods can yield highly accurate solutions when they converge.In contrast, while direct methods can handle a wide range of problems, the accuracy of their solutions depends on the discretization grid and initial guess [17,19].Furthermore, when dealing with NLP, we often fnd ourselves needing to employ black box software, with no access to its internal workings.Consequently, we are strongly motivated to adopt a hybrid indirect method for solving the optimal control problem, comprehensively analyzing the dynamical system under control.To achieve this, we prove the global convergence of the hybrid indirect solver.
Tis paper is organized as follows.In Section 2, we propose two diferent dynamical systems in the absence and presence of treatment with a novel generalized vasculature's carrying capacity.In Section 3, an optimal control therapy is presented that incorporates both anti-angiogenic and radiotherapy as two diferent control variables.Additionally, we propose a novel hybrid indirect method to solve the corresponding problem.Section 4 presents the numerical examples, while Section 5 discusses the results.Concluding remarks are provided in Section 6.

Problem Formulation
Te simplifed geometry for cancer treatment includes three types of cells: healthy cells, cancer cells, and endothelial cells.In this paper, a fact on the mathematical model is on the rate of nutrient consumption and the difusion of nutrients and oxygen to various parts of the tumor and vascular blood vessels.More precisely, we focused on the interaction between cancer cells and endothelial cells.To describe their interaction, a 2-compartment, cell population-based model is employed.Te model considers the cancer cell volume, denoted as p, and the endothelial cell volume, denoted as q, as variables.In the following subsections, we will discuss tumor dynamics, propose a general carrying capacity dynamics, and present its general dynamical systems in the absence and presence of treatments in full detail.

Tumor Dynamics.
To describe tumor growth dynamics, there are three growth models: exponential growth, logistic growth, and Gompertzian growth [9].In 1825, Benjamin Gompertz introduced and applied his growth and mortality law, which is often used to describe the growth of animals and plants, as well as the quantity or volume of bacteria and cancer cells over time [20].Te Gompertzian growth model ofers several advantages compared to the exponential and growth models.Exponential and logistic growth models failed to accurately describe the experimental data, whereas the Gompertzian model exhibited remarkable descriptive capability [21].Above all, the Gompertzian growth model efectively captures how tumor cells rely on nutrients, oxygen, and space in a continuous process.As the tumor grows, the availability of these essential resources gradually diminishes, causing the growth rate to slow down until the tumor reaches its maximum size [22].Hence, we consider the Gompertzian growth model in the form of an ordinary diferential equation (ODE): where ξ denotes the tumor growth parameter.Equation (1) implies that a tumor shrinks when p(t)/q(t) > 1.However, when p(t) tends to zero, the Gompertzian growth model will sufer from singularity.Tus, the Gompertzian model is not an adequate description of very small tumor volumes [23].
Note that when tumor grows, we have p(t)/q(t) < 1 and in medical studies [24,25], the Gompertzian growth model is responsible.

General Carrying Capacity
Dynamics.Te dynamical model for endothelial cells is represented by a balance equation between the stimulatory efect S(p, q) and the inhibitory efect I(p, q), given in the form: Hahnfeldt proposed a difusion-consumption equation with quasi-steady state to describe the concentration n of stimulators and inhibitors inside the tumor as 2 Journal of Mathematics where D 2 is the difusion coefcient, s 0 is the cells' secreting rate, and c is the cells' clearing rate.Te difusionconsumption equation ( 3) can be represented in both Cartesian and spherical coordinate systems as follows.
(i) In a Cartesian coordinate system, the inhibitor or stimulator concentration at the point (x, y, z) within the tumor at time t will be denoted as n(x, y, z, t).
Te concentration n can be determined by solving the following diferential equation: (ii) In a spherical coordinate system ( ρ,  θ,  φ) [26], Hahnfeldt considered the concentration as an explicit function of the tumor radius, assuming a symmetric tumor for simplifcation.Consequently, the concentration n at a specifc point within the tumor, located at a distance  r from the center, is introduced by solving the following diferential equation: After solving equation ( 5), Hahnfeldt determined the interior concentrations of the inhibitor n I ( r) and stimulator n S ( r) at a specifc point within the tumor as follows: where  r 0 represents the initial tumor radius.By considering equation (6), it is concluded that the inhibitor would have a signifcant impact on the target endothelial cells within the tumor, leading to growth proportional to the tumor surface or (volume) (2/3) .In contrast, the efect of stimulators would remain relatively independent of the tumor or vasculature size [7].However, by utilizing integration in spherical coordinate system [26], we drive the total stimulation and inhibition concentration inside the tumor as follows: Based on equation (7), we can conclude that the total stimulating efect is in accordance with either  r 3 0 or (volume) (3/3) .Tus, the general form of the total stimulating efect can be proposed as follows: where birth rate b > 0 and α s + β s � (3/3).Similarly, according to equation (8), the fnal expression of the total inhibition efect is in accordance with either  r 5 0 or (volume) (5/3) .Tus, we can propose the general form for the total inhibition efect as where α i + β i � (5/3) and death rate d > 0. Hahnfeldt made the assumption that the inhibition term prevents tumor cell production, which has an impact on endothelial cells.According to the total concentration of inhibition equation (8), there can also be other cases in which the inhibition term, in addition to preventing tumor cell production, also inhibits the endothelial cells production, thereby impacting endothelial cells once again.Tus, without losing generality, we propose the general carrying capacity dynamic as From equation (11), we understand that working with only the tumor radius means working with the smallest geometrical dimension.Consequently, there are 4 choices for the stimulating efect and 6 choices for the inhibiting efect as follows: According to the multiplication principle, the carrying capacity model can be written in 24 diferent ways, with each model being consistent with the actual tumor characteristics.For example, Hahnfeldt gives α s � 1, β s � 0, α i � (2/3), and Ergun et al. [8] selected questionable values: α s � 0, β s � (2/3), α i � 0, and β i � (4/3), which are inconsistent with the actual tumor characteristics.Te reason is that the inhibitor factor must be proportional to the tumor surface, while Ergun considered it to be proportional to the tumor radius (see equations ( 12) and ( 13)).Unfortunately, Kienle et al. [10], without providing any justifcation, made two mistakes: (i) using the above questionable values of the Ergun model and (ii) using an inconsistent convex combination of the Hahnfeldt model (proportional to tumor surface) and the Ergun model (proportional to tumor radius).From a mathematical point of view, with knowledge of the authors, the only possible selections for α s , β s , α i , and β i are the 24 cases given in this paper, and any selection except these 24 cases is questionable.

General Dynamical System in the Absence of Treatment.
Here, we propose and provide a description of the dynamical system under consideration in the absence of treatment for cancer cells and endothelial cells: where α s + β s � (3/3) and α i + β i � (5/3).In the rest of this subsection, we discuss the existence and uniqueness solution for dynamical system (15).

Existence and Uniqueness Solution for
General Dynamical System.Te purpose of this subsection is to establish the existence and uniqueness of the solution for the general dynamical system (15).To achieve this goal, we need to demonstrate measurability, locally Lipschitz, and locally integrable.Subsequently, we will apply the existence and uniqueness theorem as presented and proven by Sontag [27].
Lemma 1 (measurable).Assume that dynamical system (15) is in the form _ x(t) � f(t, x(t)) such that f: I × X ⟶ R 2 , I is an interval from t 0 initial time to t f fnal time, and X � (p, q)|p, q ∈ R + is an open subset in R 2 .If the following properties hold, then for all values of α s , β s , α i , and β i that satisfy equations (12) and (13), f(t, x(t)) is measurable.
Proof.(i) Consider an arbitrary set x 1 � (p 1 , q 1 ) in X and substitute in the dynamical system (15); we have vector function as follows: We know that constant real-value functions are measurable.Tus, the vector function f(., x) is measurable (see [26]; Section 3 Teorem 6).

□
Proof.(ii) By substituting an arbitrary time t 1 ∈ I in the dynamical system (15), the vector function f(t 1 , .) is similar to f(., x 1 ).We know that constant real-value functions and their linear combinations are continuous.Tus, the vector function f(t, .) is continuous.

□
Lemma 2 (locally Lipschitz).Suppose Lemma 1 holds.Ten for all values of α s , β s , α i , and β i that satisfy equations ( 12) and (13), f is locally Lipschitz on X. Proof.To prove locally Lipschitz, we must show that there are for every x 0 ∈ X a real number δ > 0 and a locally integrable function α: References for * and * * are [7,8].Te absorbed energy per unit mass of tissue is typically measured in units of gray (Gy). 4 Journal of Mathematics where B δ (x 0 ) is a ball with radius δ and center x 0 contained in X.For this purpose, let x: � (p 1 , q 1 ), y: � (p 2 , q 2 ), and f i be the ith element of vector function f for i ∈ 1, 2 { }.We have According to the Lemma 1, vector function f(t, .) is continuous and B δ (x 0 ) is bounded and is not contained (0, 0), so the following corresponding bounds can be considered.For i � 1, we have: If the elements of x and y do not tend towards zero, then in this case, we can defne , and similarly for i � 2, we have: ), we have: Tus, f is locally Lipschitz.

□ Lemma 3 (locally integrable). Suppose Lemma 1 holds. Ten f is locally integrable on t.
Proof.To prove that f is locally integrable, we must show that the following inequality holds: where β is a locally integrable function.Consider an arbitrary set x 0 � (p 0 , q 0 ) in X and f i as the ith element of vector function Tus, by introducing β(t) � C + e t , the lemma is proved.
Te presented Lemmas 1, 2, and 3 are the assumptions for Teorem 54 in [27], and by applying the aforementioned theorem, the existence and uniqueness of the solution for the general dynamical system (15) are concluded.
Te following matters must be considered if one wants to use the proposed general dynamical system (15).
(i) Scale.Te Hahnfeldt model is based on serious asymptotic analysis of the difusion-consumption equation underlying tumor anti-angiogenesis and makes the reasonable assumption that the inhibitory factor in the dynamic is proportional to the tumor surface, while the Ergun model questionably scales down these to the tumor radius.Te proposed general dynamical system (15) works with all the possible scaling features (see equations (12) and Journal of Mathematics ( 13)).However, the authors in Section 5 show that there are better selected values for α s , β s , α i , and β i compared with the Hahnfeldt model [7].(ii) Medical Relevance.In the Hahnfeldt model and our proposed general dynamical system, the tumorinduced inhibition parameter d has mm − 2 day − 1 dimension, while in the Ergun model, d has mm − 1 day − 1 dimension.Indeed, it is not feasible to consider two diferent physical concepts for a given biological phenomenon when one is measured in terms of death rates per unit length and the other in terms of death rates per unit area.
□ 2.4.General Dynamical System Modeling with Anti-Angiogenic Treatment.Te process of stimulating the formation of new blood vessels and capillaries in order to supply the tumor with the necessary nutrients is called tumor angiogenesis [11].In contrast, anti-angiogenic treatment prevents the cancer cell nutrition by eradicating the current tumor blood vessels while restricting the formation of new blood vessels [28,29].Initially, it was assumed that anti-angiogenic therapy might not exhibit toxicity compared to other chemotherapeutic agents due to the genetic stability and quiescence of endothelial cells under normal physiological conditions, as well as the selectivity of targeted drugs.However, this assumption proved to be a miscalculation [30].In Figure 1, it is shown that tumors initiate angiogenesis through the activation of vascular endothelial growth factor (VEGF).Over 60 anti-angiogenic agents (AAs) have existed in clinical trials in the US since 2006.For example, endostatin inhibits VEGF, disrupting the endothelial cell growth that forms the lining of the newly developed blood vessels [31][32][33].
We consider a control variable u that represents the AA dose and is taken as a Lebesgue-measurable function in a compact interval [0, u max ] with u max denoting the highest dose.Tus, from equation ( 2), the dynamical model for endothelial cells volume under anti-angiogenic treatment yields [9] _ q(t) � S(p, q) − I(p, q) − cq(t)u(t), (24) where c denotes the anti-angiogenic elimination parameter.Based on the above discussion and the proposed general dynamical system (15), the ODE system for anti-angiogenic treatment as monotherapy is given as follows: Numerous medical studies confrm that once monotherapy is halted, the tumor will grow back [34][35][36].Tus, anti-angiogenic therapy is not efcient as a stand-alone treatment.Still, in combination with other therapies that destroy cancer cells, such as chemotherapy or radiotherapy, it can enhance their efect and lead to synergistic benefts [6,8].

General Dynamical System Modeling with Radiotherapy.
For all types of tumors, radiotherapy is a common therapeutic approach.It is carried out in 60-70% of newly diagnosed cancer patients, either as monotherapy or as a combination treatment with surgery, chemotherapy, and anti-angiogenic therapy, which are useful [37,38].We consider a control variable w that represents the amount of the radiotherapy agent (RA) rate in a compact interval [0, w max ] with w max denoting the highest dose rate.Wein, in the year 2000, modeled cancer cells radiation damage by [39] − p(t) α + β where α and β are cancer cell radiosensitivity parameters, ρ is the tumor repair rate, and the fnal time t f is fxed in every instance within this paper.Note that is the solution to the linear ODE given by Tus, the term that quantifes the radiation damage to the cancer cells can be written in the form Notice that the radiation has a destructive efect on endothelial cells, so the damage to endothelial cells is given by where η and δ are endothelial cell radiosensitivity parameters.In the literature, often the efects of radiotherapy on healthy cells, cancer cells, and endothelial cells are modeled by a system of ODEs [40].Based on the preceding discussion and the proposed general dynamical system (15), we write: with initial conditions p(t 0 ) � p 0 , q(t 0 ) � q 0 , and r(t 0 ) � r 0 .
6 Journal of Mathematics

General Dynamical System with Combined Terapy.
Te use of diversifed methods in anti-cancer therapy ofers a broader range of options for clinical treatment and enables the formation of robust partnerships [30].For this purpose, the general dynamical system with control variables u(t) and w(t) for cancer therapy using the combination of antiangiogenic treatment and radiotherapy is in the following form: where the parameter values and notation descriptions can be found in Table 1.Te classifcation of the general dynamical system (32) using equations ( 12) and ( 13) leads to 24 distinct types of dynamical systems (see second column in Table 2).

Optimal Control Problem for a General Dynamical System
Tis paper aims to determine the optimal dosage scheme of AA and RA that minimizes tumor volume while causing the least harm to healthy cells.According to Food and Drug Administration recommendations and clinical discussions of the efectiveness of some AAs, it is reasonable to use the lowest possible dose [6,41,42].On the other hand, in radiotherapy, the healthy cells around the cancer cells are almost destroyed, and the side efects depend on accuracy and the amount of radiation dose [43].Hence, the importance of optimal control methods in this context is signifcant.

Optimal Control Formulation.
Based on the general dynamical system discussed in Section 2, one faces both mathematical and clinical challenges.Regarding the clinical aspect, treatment strategies for tumor reduction are encountered during treatment and the fnal time, where the fnal time is fxed.From a mathematical point of view, both linear and quadratic terms can be employed to control doses in the objective function.Numerous studies have demonstrated that administering continuous and low doses of therapeutic drugs leads to improved outcomes [44].Ledzewicz et al. [45] investigated linear controls and determined that optimal control is a singular arc, while bang-bang control fails to be optimal or desirable due to the need for dose continuity [12].Terefore, we consider the quadratic term in the objective function to control doses, aiming to achieve continuous and low doses, as demonstrated in a previous study [12].To ensure that the tumor volume remains nonnegative without introducing additional constraints, we express the tumor volume in a quadratic form within the objective function.Te following optimal control problem (OCP) is considered as an optimal regulator problem [46]: Here, according to control theory, X � (p, q, r) is a state space, and U � (u, w) is a control space.

Solving Optimal Control Problems: Direct and Indirect
Methods.Tere are two main numerical approaches to solve optimal control problems, depending on whether the problem is frst optimized then discretized or vice versa.Te purpose of the indirect method (optimize-then-discretize) is to fnd a solution to the ODE system resulting from Pontryagin's minimum principle.Te direct method, a more recent approach, involves discretizing the optimal control problem and then fnding the optimal solution to this discrete problem, which transforms into a fnitedimensional problem of nonlinear programming.Each of these two approaches has advantages and disadvantages [47,48].In recent years, numerous studies on cancer therapy with combinations of anti-angiogenic and radiotherapy treatments are done by direct methods.Using the direct method yields an NLP, which is solved by the interior-point optimizer (IPOPT) [49][50][51][52] via the Applied Mathematical Programming Language (AMPL) [53].On the other hand, Kienle et al. [10] applied the sequential quadratic Hamiltonian (SQH) method as the indirect method, which was frst introduced in [48,54], to solve the related optimal control problem.Here, in this article, we propose a novel hybrid indirect method to solve the corresponding problem.

Optimality Conditions for the Optimal Control Problem.
We propose an optimization algorithm based on the indirect method for OCP (33).For this purpose, we apply the optimality necessary condition by using PMP.Te general form for OCP (33) is as follows: with initial conditions X t 0  � p 0 , q 0 , r 0 , where U ad is admissible control space and is given by We defne the Hamiltonian function, which involves the integrand of the performance index and the costate variables with the right-hand side of the ODEs (34) in the form [48] H � L(X(t), U(t), t) + λ T f(X(t), U(t), t), (36) where λ represents the costate space.Specifcally, the Hamiltonian function for OCP (33) is as follows: Journal of Mathematics where λ p , λ q , and λ r are the corresponding costate variables for the state variables p, q, and r, respectively.Canonical Hamiltonian equations [46] for the OCP (33) using PMP yield with transversality conditions λ p t f   � τ p p t f  , λ q t f   � 0, and λ r t f   � 0.
Te general dynamical system (32) and canonical Hamiltonian equation (38) form a two-point boundaryvalue problem (TPBVP).If the boundary conditions are all given at either t 0 or t f , one would numerically integrate the reduced diferential equations to derive X * (t) and λ * (t), t ∈ [t 0 , t f ] [48].Unfortunately, the transversality conditions are dependent on state variables.Moreover, initial and transversality conditions are split, so the TPBVP cannot be solved in a standard manner [48].Terefore, it is motivated to apply an iterative numerical technique to overcome these difculties.By the hybrid technique here, we apply the steepest descent as an indirect method [48], conjugate gradient for search direction [55], and Armijo technique for the step length [56].Although each of the three mentioned techniques exists in the literature, by integrating and implementing them, we have designed and developed a hybrid method to solve the optimal control problem.Tis hybrid method exhibits remarkable performance when compared to direct and indirect algorithms, which will be discussed in detail in the numerical results section.

Armijo Conjugate Steepest Descent Solver.
Te steepest descent method is an iterative numerical technique used for determining the minimum of a diferentiable function and is employed as an indirect method for solving the general dynamical system (32) and canonical Hamiltonian equation (38).Te main idea for solving them is to utilize the framework of the steepest descent method [48].To search the direction, we use the conjugate gradient method as an exact line search [55].Moreover, in order to determine distinct step lengths for the two control variables, we employ the Armijo technique [56] as an inexact line search method.Te proposed hybrid method is called the Armijo conjugate steepest descent (ACSD) method.Here, the conjugate gradient direction is similar to the path opposite to ∇J with respect to the control variables u and w.For simplicity, in each step, the general dynamical system (32) and canonical Hamiltonian equation (38) are considered as the following form.TPBVP: and TPBVP (39) satisfes the boundary conditions X i t 0  � p t 0 , q t 0 , r t 0  , Nominal control history in each step will update by direction ψ i and step length ϵ i : 10 Journal of Mathematics For better readability, in the following, we defne g i : � (zH/zU)(X i (t), λ i (t), U i (t), t).Te value of ψ i is determined through the conjugate gradient method using the following equation [55]: where ζ i is a scalar in the following form: Furthermore, the computation of distinct step lengths ϵ i � (ϵ u i , ϵ w i ) for the two control variables involves an inexact line search using the Armijo technique, as outlined in [44,56]: where κ ∈ (0, 1) is known as Armijo backtracking parameter.Restrictions on AAs and RAs for nonnegativity and admissibility force us to put the following step to Algorithm 1.
where U max � (u max , w max ).Finally, the iteration is stopped when the criterion ‖ψ i ‖ ≤ ϑ is satisfed, where ϑ is a preselected positive constant and we have: Hence, we summarize the Armijo steepest descent method for the general dynamical system (32) and canonical Hamiltonian equation (38) in Algorithm 1.

Convergence Analysis of the Armijo Conjugate Steepest
Descent Method.In this subsection, we present the convergence analysis of the Armijo conjugate steepest descent method.For this purpose, we will demonstrate that the conjugate gradient coefcient satisfes a sufcient descent condition, and we will prove global convergence.Here, we defne an additional condition by the Hamiltonian gradient and direction to ensure that the algorithm makes meaningful progress, as follows: It is important to mention that equation ( 47) is more general than the curvature condition in the strong Wolfe conditions (see equation (3.7b) in [44]).44) and (47), then for ‖g i ‖ ≠ 0, we have

Lemma 4. If we have Armijo conjugate steepest descent with conditions equations (
Proof.From equation ( 27) and the Cauchy-Schwarz properties, we have Dividing both sides by ‖ψ i ‖ yields Te right-hand side of the above inequality is always negative.Tus, the proof is complete.44) and ( 47), then we have
Proof.From equation (42), for i � 1, it is clear that ‖g 1 ‖/‖ψ 1 ‖ � 1, and for i ≥ 2 by multiplying g T i+1 , we have Taking the absolute value and applying the Cauchy-Schwarz properties, Terefore, we have and by dividing both sides of the above inequality by ‖g i+1 ‖‖ψ i+1 ‖, we have Journal of Mathematics Tis suggests that From equation (47), Tus, the proof is complete.
□ Theorem 6.Assuming that ψ i is generated using equation ( 42) and ζ i is obtained from equation ( 43), a sufcient condition is satisfed when where M > 0 holds.
Proof.Multiplying both sides of equation ( 42) by g T i+1 and using ζ i+1 yields and we have For the right-hand side of the above inequality, we can take the absolute value and apply the Cauchy-Schwarz properties: By employing the condition stated in equation ( 47), we obtain: Te proof is complete.Terefore, for ‖g i+1 ‖ ≠ 0, this implies that Te following theorem provides proof of the Zoutendijk conditions [57], establishing the global convergence of Armijo conjugate steepest descent with the new condition (47).□ Theorem 7. Consider the Armijo conjugate steepest descent method, where ϵ i is obtained using the Armijo inexact line search rule equations ( 44) and (47).Assuming Lemma 4, Lemma 5, and Teorem 6 hold true, then either Proof.We employ a proof by contradiction argument.If Teorem 7 is not true, then there exists a > 0, such that From equation (42), Upon squaring both sides of the above equation, we have (1) Subdividing the interval [t 0 , t f ] into N equal subintervals and considering piecewise-constant controls: Applying the assumed control U i to integrate the state equations from t 0 to t f with initial condition X i (t 0 ) and store the state matrix X i .(3) Calculate λ i (t f ) from equation (40), use U i and X i to integrate the ODE costate from t f to t 0 , and store the costate matrix λ i .(4) Evaluate ψ i from equation ( 42) and store as a vector.
(5) if‖ψ i ‖ ≤ ϑ then (6) terminate the iterative procedure and output U i and X i .( 7) else (8) fnd the ϵ i from Armijo technique for the two control variables such that then replace U i by U i+1 and return to step 2. (11) end if ALGORITHM 1: ACSD method.12 Journal of Mathematics (67) Dividing both sides by (g T i+1 ψ i+1 ) 2 yields (68) Applying equation ( 43) yields From equation ( 57) of Lemma 5, we have and based on equation (48) of Lemma 4, we have On the other hand, from equation ( 59), we have: For the right-hand side of the above inequality, we can take the absolute value and apply the Cauchy-Schwarz properties: By employing the condition stated in equation ( 47), we obtain: We know that (2 − 1/ω) > 1, and thus Substitute equation (75) in equation ( 71), and we get Dividing both sides by (2 − (1/ω)) yields By applying equation ( 57), we have Te right-hand side of the above inequality is always positive; then Tis contradicts equation (63) in Teorem 6. Terefore, the proof is completed.
Te questions to be asked at this stage are these: (i) Can we determine values for α s , β s , α i , and β i from the values of second column in  2. Note that dynamical system number 21 in Table 2 is known as the Hahnfeldt dynamical system.Figure 2 depicts the comparison among the four selected dynamical systems in Example 1, including tumor volume, vascular carrying capacity, and the amounts of anti-angiogenic and radiation agents.In Figure 2(b), the optimal vascular carrying capacity of dynamical system 21 initially increased slightly and then decreased.Te reason for this event is the role of vascular carrying capacity in the inhibiting factor.Strictly speaking, by reducing the efect of vascular carrying capacity, the response to anti-angiogenic treatment also decreases.In Figure 2(c), considering the daily dosage limit for the anti-angiogenic drug, dynamical system 21 exhibited the highest consumption of the drug between days 5 and 8.However, its performance in terms of the fnal tumor volume was not as efective.In contrast, dynamical system 19 had lower anti-angiogenic drug usage throughout the treatment period and achieved the smallest fnal tumor volume.In Figure 2(d), in accordance with the limitation of the daily radiation therapy dose, as expected, the maximum allowable daily dosage has been administered.
Among the 24 possible dynamical systems listed in Table 2, the four dynamical systems (13, 19, 20, and 21) with the lower costs have been selected.Example 2. Tis example aims to compare results obtained from both direct and indirect methods.Here, the treatment strategy focuses on reducing tumor volume during treatment using the lowest possible anti-angiogenic and radiation agents as optimal regulator problem in Lagrange form.As in Example 1, we employ the four selected dynamical systems in order to make a regular comparison between two treatment strategies while serving the purpose of abbreviation.Other parameter values in performance index are given by t f � 10, τ p � 0, θ p � 1, θ u � 1, and θ w � 1.
Figure 3 depicts the comparison of tumor volume, vascular carrying capacity, and the amounts of antiangiogenic and radiation therapy for the four selected dynamical systems in Example 2. According to Figure 3(c), an anti-angiogenic drug is proposed with a bang-bang property for each selected dynamical system (13, 19, 20, and 21).Given the constraint on the daily permissible dose of radiation therapy within the acceptable range [58], it was as expected that the optimal radiation therapy dose would be selected as the highest possible value throughout the treatment (shown in Figures 2(d

. Discussion of Results
Based on the numerical results reported in Section 4 regarding the dynamical system, solver method, and treatment strategy, we have arrived at the following fndings.

Dynamical System.
Based on Example 1 and its numerical results, the discussion for the 24 distinct versions of a general dynamical system are as follows: For stimulatory efect, one can say the following."By increasing α s (which increases the tumor volume infuence) or decreasing β s (which decreases the carrying capacity infuence), all values of J, p(t f ), and q(t f ) experience a signifcant decrease.Tis implies that if the stimulating efect of a tumor exceeds its volume in comparison to the vascular carrying capacity, there will be a substantial reduction in p(t f ), q(t f ), and the overall treatment cost.Consequently, low-cost treatments can be employed where tumor volume plays a crucial role in the stimulation efect." For the inhibitory efect, we have the following."Trough the α i increase (which enhances the tumor volume infuence) or the β i decrease (which diminishes the carrying capacity impact), almost signifcant increases are observed in the values of J, p(t f ), q(t f ), and D AA .Tis indicates that if the inhibitory efect of a tumor exceeds its volume in relation to the vascular carrying capacity, there will be substantial increases in tumor volume, vascular carrying capacity at the end of the treatment, and the overall treatment cost.Terefore, high-cost treatments can be justifably applied where tumor volume plays a critical role in the inhibition efect.

Journal of Mathematics
For the dynamical systems, the following aspects can be considered.
"Te Hahnfeldt dynamical system, among the four above selected dynamical systems, has consumed the highest amount of anti-angiogenic drugs.It can be concluded that the Hahnfeldt dynamical system compared to the three other dynamical systems (13, 19, and 20) has more resistance to treatment.Tus, in order to respond to the arguments in [40], we propose to the physician to use dynamical system number 19 that has less cost, tumor and carrying capacity volumes at fnal time, and total dose of anti-angiogenic agent among other dynamical systems."16 Journal of Mathematics zero, then the IPOPT-AMPL encounters "cannot compute 0/0" error.Tis error is due to checking all grid points for the optimal value of NLP.Another disadvantage of the AMPL-IPOPT method is that it encounters a "slack too small" error when applied to dynamical systems 4, 5, 6, 11, and 12 with objective function in Example 2. Tus, there are limitations in the AMPL-IPOPT implementation, making it unsuitable for solving many optimal control problems, which has been confrmed and further analyzed in [59].(iii) Te Armijo conjugate steepest descent method almost provides a better solution for optimal control problems compared to the other two methods.

Treatment Strategy.
Upon comparing the numerical results presented in Examples 1 and 2 to investigate the treatment strategies, one can conclude the following (see Tables 2 and 3)."In the treatment strategy for reducing the cancer cell volume during treatment, a large amount of anti-angiogenic drug is used compared to the other treatment strategy.However, this strategy does not show a signifcant change in the reduction of tumor and vascular carrying capacity volumes.Additionally, the use of an anti-angiogenic drug with a bang-bang property contradicts the low-dose and continuous treatment approach for each of the selected dynamical systems (13, 19, 20, and 21) (see Figure 3(c)).Consequently, the administration of high doses of antiangiogenic drug comes with its own side efects [60], which can be clinically harmful.However, if it is necessary to reduce the vascular carrying capacity volume, one can consider adopting the treatment strategy of decreasing tumor volume during treatment by using a higher dosage of anti-angiogenic drugs, while also accepting its associated side efects.By implementing this strategy, the probability of tumor growth after the treatment signifcantly decreases.On the other hand, for optimal control problems with dynamical system 19, the tumor volume at the fnal time reaches levels between 19 and 22, while the carrying capacity volume at the fnal time reaches levels between 235 and 240, indicating strong medical sense."

Conclusion
In this paper, cancer therapy with anti-angiogenic and radiotherapy treatments is discussed.Te general dynamical system (15) is proposed.Te hybrid indirect solver is proposed which successfully solves the related optimal control problems (33).Mathematically, the existence and uniqueness of the solution, as well as the convergence of the hybrid indirect method, were proven.Te fndings of this paper are as follows: (a) Te low-cost treatments can be employed, where tumor volume plays a crucial role in the stimulation efect.
(b) Te application of high-cost treatments can be justifed when tumor volume signifcantly infuences the inhibitory efect.(c) Ofer to physicians and scientists: among the available options, there is at least one dynamical system that demonstrates lower costs, reduced tumor and carrying capacity volumes at the fnal time of treatment, and a lower total dose of the antiangiogenic agent compared to the other dynamical systems.(d) Te Armijo conjugate steepest descent method almost ofers a better solution for optimal control problems compared to the IPOPT-AMPL and SQH methods.(e) Comparison of two treatment strategies (see Examples 1 and 2) makes it clear that the medical experts can use diferent strategies based on the situation and weights that they put on the way of the treatment.
In the next future work, we will mathematically present the reasons behind the resistance observed in reducing the carrying capacity of blood vessels when increasing the use of anti-angiogenic drugs.

Figure 1 :
Figure 1: Angiogenic process and intense growth of cancer cells and endothelial cells.

Figure 2 :Figure 3 :
Figure 2: Optimal regulator problem in Bolza form: treatment strategy to reduce tumor volume at the fnal time based on data in Example 1 for the selected dynamical systems (13, 19, 20, and 21) using the Armijo conjugate steepest descent method.(a) Optimal cancer cell volume.(b) Optimal endothelial cell volume.(c) Optimal anti-angiogenic agent.(d) Optimal radiation agent.

Table 1 :
Notation and numerical values for model parameters and state and control variables.

Table 2 :
(33)rical results for 24 distinct versions of a general dynamical system for the optimal control problem(33), with a treatment strategy in Example 1.

Table 2
[58] and total RA doses (computed as L 1 -integrals) denoted as D RA utilized in the treatment.Goedegebuure et al.[58]confrmed that fractionated low-dose radiotherapy, i.e., daily fractions of up to 2 Gy, has a positive efect on cancer cells.Tus, in the following two diferent examples, we bound the w max to be less than or equal to 2 Gy.Te results have been conducted on an Intel Core i5 CPU 6th Generation at 2.4 GHz, with 4 GB RAM.� 10, τ p � 1, θ p � 0, θ u � 1, θ w � 1, and others in Table1.Our proposed Armijo conjugate steepest descent method is used as the solver, and the results are reported in Table Table 3 for the analysis of the solver methods, we have reached the following conclusions: (ii) Te execution speed of the IPOPT-AMPL method is very slow and depends on the lower bound values of state variables.If the lower bounds of the state variables are chosen to be greater than or equal to

Table 3 :
Numerical results obtained by three diferent solvers (ACSD, SQH, and IPOPT-AMPL) applied to selected dynamical systems with a treatment strategy in Example 2.