Abstract

In this article, a modified gradient approach, based on the adjoint method, is introduced, to deal with optimal control problems involving inequality constraints. A common way to incorporate inequality constraints in the adjoint approach is to introduce additional penalty terms in the cost functional. However, this may distort the optimal control due to weighting factors required for these terms and raise serious concerns about the magnitude of the weighting factors. The method in this article avoids penalty functions and can be used for the iterative computation of optimal controls. In order to demonstrate the key idea, first, a static optimization problem in the Euclidean space is considered. Second, the presented approach is applied to the tumor anti-angiogenesis optimal control problem in medicine, which addresses an innovative cancer treatment approach that aims to inhibit the formation of the tumor blood supply. The tumor anti-angiogenesis optimal control problem with free final time involves inequality and final constraints for control and state variables and is solved by a modified adjoint gradient method introducing slack variables. It has to be emphasized that the novel formulation and the special use of slack variables in this article delivers high accurate solutions without distorting the optimal control.

1 Introduction

In the realm of complex systems and dynamic processes, the quest for efficiency and optimal performance has always been a driving force. Optimal control theory offers a powerful framework for tackling a wide range of problems. Whether it is a robot, a spacecraft, or an application in medicine, the ability to navigate and manipulate systems in the most effective way possible holds immense value.

The roots of optimal control theory can be traced back to the early 20th century, while significant advancements were made during the Apollo moon program in the 1960s [13]. Among the works mentioned, the works by Kelley, Bryson, and Ho are today fundamental for modern optimal control theory.

Despite the long history of optimal control, there is still research potential by combining principles from optimal control theory, optimization, and novel computational methods. Moreover, algorithms are becoming more robust to account for the growing complexity of problems. Various authors have already discussed on the complexity and computational strategies of optimal control theory [46].

Beside classical mechanics, optimal control is becoming increasingly important in modern medicine and can make a significant contribution to cancer treatment [7]. It involves using mathematical models to find the best ways to treat cancer. Several mathematical models have been proposed in the open literature to study tumor anti-angiogenesis [812].

One of the first biologically validated models of tumor growth dynamics was introduced by Hahnfeldt et al. [8]. This model serves as the basis for subsequent modifications and simplifications aimed at better understanding the dynamic properties of the underlying mechanisms. In the article of Ledzewicz and Schättler [12], the problem of minimizing tumor volume by administering predetermined amounts of anti-angiogenic and cytotoxic agents is examined. Optimal and suboptimal solutions are presented for four variations of the problem.

Another medical application can be found in the work of Gao and Huang [13], where an optimal control model of tuberculosis is presented to minimize both disease burden and intervention costs. More recently, the COVID-19 pandemic posed significant challenges to health systems in managing the disease. The work by Hametner et al. [14] discusses epidemiological models, allowing the prediction of infection dynamics using a flatness-based design.

Despite the great variety of mathematical models, the fundamental premise of all optimal control problems is to find the optimal input that steers a system toward desired goals, while accounting for various constraints. Handling constraints on the control variables, e.g., Graichen et al. [15] have systematically transformed an inequality constraint into a new equality constraint by means of saturation functions. In this article, we present a novel approach to address optimal control problems including inequality constraints by introducing a modified gradient method based on the adjoint technique. The adjoint approach, as, e.g., [1620], is meanwhile a well-established method for efficient gradient computation in optimization procedures.

Incorporating inequality constraints in the adjoint approach could be realized by the introduction of penalty terms in the cost functional. The novel method proposed in this article is an enhancement of Refs. [21] and [22] in order to circumvent the use of penalty functions. The basic idea of the proposed method traces back to a transformation of the inequality constraint to an equality constraint by the introduction of so-called slack variables [23 and 24]. The adjoint method coupled with the use of slack variables allows for an iterative computation of the optimal control, which ensures highly accurate solutions without biasing the optimal control. In order to show the robustness of the algorithm, the presented approach is applied to the tumor anti-angiogenesis optimal control problem. By avoiding the need for penalty terms and their associated weighting factors, the proposed method offers a promising avenue for addressing optimal control problems with inequality constraints effectively.

2 A Modified Gradient Approach for Optimization Problems With Inequality Constraints

Let us first consider the static optimization task
(1)
On the one hand, the function f(x) should be minimized and on the other hand, the variables x must fulfill the inequality constraint g˜(x)0 in each iteration step. In this article, we now introduce a handsome modification of the gradient approach, as presented in Ref. [22], for solving constrained optimization problems with inequality constraints. If we introduce the slack variable σR, tracing back to the Karush-Kuhn-Tucker conditions, the inequality constraint in Eq. (1)2 can be rewritten to an equality constraint. This approach is similarly employed in sources [23,24], where inequality constraints are also transformed into equalities. Hence, the equality constraint reads
(2)
or, respectively
(3)
Introducing y=(xT,σ)T, the variations of the function f in Eq. (1)1 and of the constraint g in Eq. (3) are given by
(4)
in which :=(/x1,/x2,,/xn,/σ)T. As a descent direction, approaching the minimum of f and g=0, we introduce a linear combination of the two gradients f and g
(5)
in which κ>0 is the update step size. For the determination of the multiplier μ, we claim that the descent direction should always point to the hypersurface g(y)=0 by setting
(6)
In this case, δg, i.e., the change of g, will be positive if g<0 and negative if g>0. Inserting Eq. (5) into Eq. (4)2 and equating the resulting term with Eq. (6) yields
(7)
From this, μ can be computed as
(8)
Hence, the combined search direction (cf. Eq. (5)) is given by
(9)

As a simple example, we consider the objective function f(x,y)=2x2+y2 and the inequality constraint g˜(x,y)(x+a)2+(y+b)21. Choosing ε=0.9 and κ=0.42, an update of x, y, and σ, which reduces f and incorporates g˜0, can be computed from Eq. (9). In order to show the robustness of the proposed method, we consider two different cases summarized in Figs. 1 and 2. The contour lines of g˜ are shown in gray, while the contour lines of f are depicted in black. Moreover, the contour line g˜=0 indicates the border of the restricted area and is shown as a dashed line. Note that the update becomes very small near the optimum; thus, the gradient has been scaled in both figures.

Fig. 1
Visualization of the optimization procedure with an inequality constraint along the search direction
Fig. 1
Visualization of the optimization procedure with an inequality constraint along the search direction
Close modal
Fig. 2
Visualization of the optimization procedure with a minimum inside the restricted area
Fig. 2
Visualization of the optimization procedure with a minimum inside the restricted area
Close modal

In the first case, the constraint g˜ is defined by a=2 and b=1 and the optimization procedure is initiated with x0=4, y0=2, and σ0=2. Here, the inequality condition would be violated if the steepest descent of f to the minimum were followed. As depicted in Fig. 1, the method converges to the global minimum (x*,y*)=(0,0) and σ*=2, circumventing the forbidden area g˜<0. In the second case, the constraint g˜ is given by a=1/2 and b=1/2 and the optimization procedure is started again with x0=4, y0=2, and σ0=2. Thus, if the global minimum of f in (0,0) lies in the region of the inequality constraint, the method converges to the constrained optimum as depicted in Fig. 2. At the optimal point, given by (x*,y*)=(0.113,0.290) and σ*=0, the contour line of f(x,y) is tangent to the boundary of the constraint g˜(x,y). In summary, the method converges to the optimal solution in both scenarios, i.e., for constraints along the search direction and for constrained minima.

In the third case, the constraint g˜ is given by a=1/2 and b=3/2 and the initial values read x0=4, y0=2, and σ0=2. Note that in this case the global minimum lies exactly on the boundary of the inequality constraint g˜=0. As shown in Fig. 3, the method converges to the global minimum (x*,y*)=(0,0) and σ*=0.

Fig. 3
Visualization of the optimization procedure with a minimum on the boundary of the restricted area
Fig. 3
Visualization of the optimization procedure with a minimum on the boundary of the restricted area
Close modal

3 Optimal Control Problem

In this section, we focus on the optimization of an optimal control problem minimizing a cost functional, which includes a scrap function. Moreover, we investigate two types of constraints imposing limitations on the system behavior, encompassing final conditions and inequality constraints.

3.1 Problem Formulation.

If we are dealing with an optimal control problem, we face a similar task as presented above. In optimal control, a dynamic system is mostly described by a set of first order differential equations
(10)
with a state vector x(t)Rn, a control vector u(t)Rm, and a given initial state vector x0. The goal is to find controls u(t) that minimize an objective function measuring, e.g., energy consumption or the operation period of a process. For free final time problems, we introduce a cost functional of the form
(11)
where S(x(tf),tf) denotes a final cost term and h(x(t),u(t)) describes the cost per time. In order to impose final conditions on the states x, we require the system to satisfy
(12)
in which r is the number of final conditions. So far, control constraints are commonly considered by penalty terms in the cost functional when using a gradient based method. However, using the presented modified gradient approach, penalty terms in Eq. (11) can be avoided. Therefore, we consider a set of inequality constraints
(13)
If we introduce a vector of slack variables σ(t)Rl, the inequality constraints can also be formulated as equation constraints
(14)
or after rearranging
(15)
Instead of Eq. (15), we may now consider
(16)
and integrate Eq. (16) over the time interval. Hence, we obtain
(17)
Note that Eq. (17) is equivalent to Eqs. (14)(16), because g(x(t),u(t),σ(t))0 and, hence, Γ(x(t),u(t),σ(t)) can only become zero if g(x(t),u(t),σ(t))=0 at every instant of time. The use of slack variables in dynamic optimization problems is also considered in Ref. [23,24], in which path constraints are transformed into integral constraints. Hence, in order to consider system limitations, the function
(18)

can now be used to put bounds on x(t) and u(t). Hence, the optimal control problem is formulated as the problem of finding controls u(t) and slack variables σ(t) in the time interval t[t0,tf] which minimize the functional in Eq. (11) and satisfy at the same time Eqs. (12) and (18). In order to apply a solution strategy to our optimal control problem, we first compute the variation of the cost functional in Eq. (11) with respect to the control signal. Moreover, we derive the variation of the final conditions in Eq. (12) and the variation of the constraints in Eq. (18), both with respect to the control signal. Finally, we consider how the variations can be combined in a reasonable way to obtain an update formula for applying a descent optimization algorithm.

3.2 Adjoint Equations.

As the first step, to derive the adjoint equations and the gradient of the cost functional, we expand J in Eq. (11) with the state equations, yielding
(19)
in which p(t) denotes a vector of arbitrary adjoint variables. It has to be mentioned that in the notation, the time dependencies have been omitted for better readability. For free final time, the first variation of the extended cost functional in Eq. (19) is given by
(20)
Note that the variations δtf and δxf are not independent, confer Ref. [25, Sec. 4.3]. For free final time problems, the variation δx(tf) describes the variation of x(t) at fixed end time tf, while δxf is the variation of x(t), including also the final time variation. Using the relation
(21)
and applying an integration by parts of the last term of the integral pT(t)δx˙(t) yields
(22)
where δx(t0)=0 as the initial conditions are fixed. Herein, hx, hu, fx, and fu denote the partial derivatives of h and f with respect to x and u. We have not imposed any restrictions on the adjoint variable vector p(t); thus far, we now define p(t) through the linear time-varying final value problem
(23)
If we introduce the total time derivative S˙(t):=Sx(x(t),t)x˙(t)+St(x(t),t), the first variation of the cost functional in Eq. (22) is given by
(24)

which shows the influence of δu(t) and δtf on δJ.

3.3 Influence Equations for the Final Conditions.

Analogous to the cost functional above, we now derive the first-order variation of the final conditions. Therefore, Eq. (12) is now extended by the state equations reading
(25)
in which one more set of adjoint variables arranged in the matrix P(t)Rn×r is introduced. Note that the integral term on the right side is zero for any matrix P(t), if the state equations are satisfied. As a first step, we compute the variation of ϕ, yielding
(26)
in which ϕx and ϕt denote partial derivatives of the function ϕ(x(t),t) with respect to x and t. Utilizing the relation from Eq. (21) and applying an integration by parts of PT(t)δx˙(t) in Eq. (26), we obtain
(27)
Now, we define P(t) by the matrix differential equations
(28)
which can be integrated backward in time. Therefore, we have to solve one set of n ordinary differential equations for each component of the final condition ϕ(x(tf),tf)=0. With P(t) from Eq. (28), and with
(29)
as the total time derivative of ϕ(x(t),t) at t=tf, Eq. (27) reduces to
(30)

showing the direct influence of δu(t) and δtf on δϕ.

3.4 Influence Equations for the Inequality Constraints.

In order to account for inequality constraints in the optimal control problem, we now deduce the variation of Eq. (17). For an arbitrary adjoint matrix Q(t)Rn×l, the function in Eq. (17) is still zero if we expand it as in Eq. (25) before, which gives
(31)
Moreover, the infinitesimal change of the function Γ is given by
(32)
After an integration by parts of QT(t)δx˙(t) and collecting all variations, we get
(33)
Now, we may choose the matrix Q(t) such that
(34)
where we get one column in Q(t) for each inequality constraint. Finally, the first variation δΓ reduces to
(35)

which directly relates the variation of the inputs δu(t), the variation of the slack variables δσ(t) and the variation of the free final time δtf with the variation of the constraints δΓ.

3.5 Gradient Combination.

Similar to Eq. (4) for a static optimization problem, Eqs. (24), (30), and (35) describe the variations of the cost functional and two sets of auxiliary conditions in the case of an optimal control problem with varying end times. Since we are dealing with first-order variations, we can combine Eqs. (24), (30), and (35) by a linear combination
(36)
in which we introduced the multipliers νRr and μRl. After inserting Eqs. (24), (30), and (35) for δJ, δϕ, and δΓ, we receive
(37)
Finally, the update formulas of the final time tf, the controls u(t), and the slack variables σ(t) are obtained by choosing
(38)
where the numbers α and β serve for conditioning the problem and κ governs the step size in the optimization procedure, requiring appropriate selection. Moreover, the multipliers ν and μ account for weighted associations between δJ, δϕ, and δΓ and can be determined similar to Eqs. (6)(8). First, we aim to improve the approximation of the condition ϕ(x(tf),tf)=0; hence, we set
(39)
where 0<ε1<2 (confer [22, Sec. 2.3] for more information). In order to compute ν, we equate Eq. (30) with Eq. (39) and insert Eq. (38) for δu and δtf
(40)
or after rearranging and introducing the abbreviations
(41)
we obtain
(42)
In a next step, we consider the inequality constraints and how we can achieve Γ(x(t),u(t),σ(t))=0. Since g(x(t),u(t))0 for an arbitrary δu(t) with t[t0,tf], we deduce from Eq. (6) that the constraints in Eq. (17) can be approached if we set
(43)
In order to compute μ, we set Eq. (35) equal to Eq. (43) and insert Eq. (38) for δtf, δu(t), and δσ(t)
(44)
in which gf=g(x(tf),u(tf),σ(tf)). Using the abbreviations
(45)
in Eq. (44), we finally receive
(46)
Hence, Eqs. (42) and (46) may be summarized in the following matrix notation for ν and μ:
(47)

The linear system of equations can be readily solved for ν and μ, which can then be inserted into the update formulas in Eq. (38) for approaching the final conditions and inequality constraints.

4 Algorithm

Finally, we provide a summary of the steps involved in the modified gradient approach to compute optimal controls with final conditions and inequality constraints. To facilitate the process, a transformation into a unit interval is beneficial due to the free final time tf. Thus, we introduce a new time coordinate, τ[0,1], related with t=tfτ to the original time coordinate. Derivatives with respect to τ are denoted as (·) and are defined as (·)=tfd(·)/dt. The procedure for solving optimal control problems can be summarized as follows:

Before applying an optimization scheme, select a final time guess tf, an initial control signal u(τ), and a slack variable signal σ(τ), where τ[0,1]. Choose the weighting factors α and β and the parameters ε1 and ε2. For the optimization problem conditioning, the scaling parameters α and β are chosen by balancing the magnitude of the update formulas δtf, δu(τ), and δσ(τ). Moreover, the values of the parameters ε1 and ε2 determine the convergence toward the final conditions and the fulfillment of the rewritten inequality constraints. The range of values for ε1 and ε2 is discussed in Ref. [22], but experience has shown that the setting of the parameters only has an influence on the convergence speed and not on the convergence success.

  1. Solve the state equations
    (48)
    with the initial condition x(0)=x0 in the interval τ[0,1].
  2. Solve the adjoint differential equations
    (49)
    with the final condition p(1)=SxT(x(1),1) from τ=1 to τ=0.
  3. Solve the influence differential equations I
    (50)
    with the final condition P(1)=ϕxT(x(1),1) from τ=1 to τ=0.
  4. Solve the influence differential equations II
    (51)
    with the trivial final condition Q(1)=0 from τ=1 to τ=0.
  5. Compute the multipliers ν and μ by solving the linear system of equations
    (52)
  6. Select κ and update the final time, controls and slack variables by
    (53)
    A constant step size κ can be chosen, but for improved convergence, utilizing a merit function is beneficial as described in Ref. [26].
  7. Check the gradient of J and the values of the final conditions ϕ(x(1),1) and inequality constraints Γ(τ).

Repeat steps 1 through 6 until the gradient of the cost functional J is sufficiently small and the final conditions and inequality constraints in step 7 are satisfied. Note that if a merrit function, as described in Ref. [27, Sec. 15.5], is used to determine the update step size κ, the method is open for standard numerical gradient methods to minimize a function. In order to improve the algorithm convergence, the nonlinear conjugate gradient method can be applied, incorporating information from a prior update step.

5 Tumor Anti-Angiogenesis Optimal Control Problem

To demonstrate the application of the modified gradient approach, we consider the tumor anti-angiogenesis problem presented in Ref. [12]. The tumor anti-angiogenesis optimal control problem was selected as an example to test the algorithm with a variety of conditions, including terminal conditions, control inequality constraints, and free final time. Although state inequality constraints are not included in this example, in the authors' experience such problems show a similar convergence for this type of constraints.

The anti-angiogenic therapy is an innovative cancer treatment approach that aims to inhibit the formation of the blood supply system necessary for tumor growth.

Most cancer chemotherapy treatments fail due to a combination of intrinsic and acquired drug resistance. While cancer cells become increasingly resistant, the drugs continue to kill healthy cells, leading to therapy failure.

Therefore, finding cancer treatment methods that overcome drug resistance is crucial. One possible therapeutic approach offers anti-angiogenesis for the treatment of tumors. Above a certain tumor size, tumors form their own blood supply by developing a vascular system through a process called angiogenesis.

Angiogenesis is a signaling process between tumor cells and endothelial cells, i.e. cells that cover the inside of blood vessels. Tumor cells release a vascular endothelial growth factor to stimulate endothelial cell growth in blood vessels. In response, blood vessels provide vital nutrients to the tumor, supporting its growth.

However, endothelial cells possess specific receptors that make them responsive to inhibitors of angiogenesis, such as the naturally occurring inhibitor protein endostatin. In pharmacological therapies, the primary objective is often to target the growth factor vascular endothelial growth factor, with the aim of hindering the formation of new blood vessels and capillaries. By impeding the development of these vital structures, the growth of the tumor can be effectively blocked.

5.1 Problem Formulation.

We consider a mathematical cell population model of a tumor, described in Ref. [12], including differential equations for the tumor volume and its vascular carrying capacity, which allow investigation of anti-angiogenic cancer treatment. The state equations governing tumor growth are given by
(54)
in which the variables v(t) and w(t) represent volumes. In detail, the variable v(t) denotes the volume of the primary tumor and w(t) designates the vascular carrying capacity, which is associated with the effective vascular support provided to the tumor. The variable u(t) denotes the control input of the system and describes the angiogenic dose rate, e.g., endostatin. Moreover, the variable ζ denotes a tumor growth parameter, b is the birth rate and d is the death rate. The parameter μ is used to consider the baseline loss of vascular support, while the term μw(t) describes the loss of endothelial cells due to natural causes. The loss of endothelial cells due to additional outside inhibition is modeled by the term Gu(t)w(t), where G is the anti-angiogenic killing parameter. In Ref. [8], the proposed model was subjected to biological validation, and parameters were estimated from medical data obtained from animal experiments in mice with implanted Lewis lung carcinoma.The tumor anti-angiogenesis optimal control problem is stated as follows: For a free terminal time tf, we are looking for the control angiogenic dose rate u(t) to minimize the functional
(55)
in which γ=0.01 is a scaling parameter and v(tf) is a scrap function to measure the final tumor size. The boundary conditions for t0=0 are given by v(0)=v0 and w(0)=w0. In order to limit the maximal cumulative dose of all angiogenic agents administered, we introduce the integral constraint
(56)
which can be rewritten as y(tf)=A, if we add the state equation y˙=u(t) to Eq. (54) with y(0)=0. Hence, the final constraint reads
(57)
Moreover, the control inequality path constraint is given by
(58)
or, respectively
(59)

in terms of the proposed algorithm.

5.2 Solution Strategy.

To demonstrate the proposed method and its algorithmic implementation, we follow the step-by-step procedure in Sec. 4 to show how the method can effectively solve the problem at hand.

  1. In a first step, according to Eq. (48), we transform the state equations of the tumor anti-angiogenesis problem into the τ-domain yielding
    (60)
    Note that Eq. (60) can be rewritten by Eq. (48) if we introduce x=(v,w,y)T. With the initial conditions x(0)=(v0,w0,0)T, Eq. (60) can be solved in the interval τ[0,1].
  2. The adjoint differential equations in Eq. (49) for p=(p1,p2,p3)T are given by
    (61)
    and can be solved with the final conditions p(1)=(γ,0,0)T from τ=1 to τ=0.
  3. Note that we are dealing with only one final constraint, given by Eq. (57), hence, the first set of influence differential equations for P=(Pij)R3×1 from Eq. (50) are given by
    (62)
    which may be also solved backward in time by starting at τ=1 with the final conditions P(1)=(0,0,y(tf)A)T. Since no objective function h is defined in the present problem, h=0 here. Hence, the right-hand side of the adjoint equations in Eq. (61) correspond to the right-hand side of the influence equations in Eq. (62).
  4. Moreover, we have no state inequality constraints, hence, the final condition for Q(τ) is given by Q(1)=0 and Q(τ) in Eq. (51) is zero for the entire τ-domain and have therefore no effect on the gradient formula.

  5. Before we can compute a combined update for the system inputs, the multipliers νR and μR2 must be determined by Eq. (52). Thus, the components of the matrix M1=(M1,ij)R1×2 in Eq. (41)1 read
    (63)
    and N1R can be determined by Eq. (41)2, yielding
    (64)
    Moreover, b1R can be computed by Eq. (41)3 and is given by
    (65)
    where the boundary term is zero, since hf=0. The components of the matrix M2=(M2,ij)R2×2 in Eq. (45)1 read
    (66)
    where M2,12=M2,21. Moreover, N2R2 can be determined by Eq. (45)2
    (67)
    Inserting in Eq. (45)3, b2R2 reads
    (68)
    Arranging Eqs. (63)(68) in matrix notation, the multipliers ν and μ can be computed by solving Eq. (52).
  6. Finally, the gradient formulas are obtained from Eq. (53), where the update for the final time tf and the control signal u(τ) are given by
    (69)
    and the updates for the slack variables read
    (70)

By adding the updates to the previous estimates, we receive an improvement by: tf+δtf as new final time, u(τ)+δu(τ) as new control signal, σ1(τ)+δσ1(τ) and σ2(τ)+δσ2(τ) as new slack variables.

5.3 Numerical Results.

For numerical computations we use the parameter set from Hahnfeldt et al. [8]. In summary, the parameters read ζ=0.084day1, b=5.85day1, d=0.00873mm2day1, G=0.15kg/mg, and μ=0.02day1. Moreover, the initial conditions are given by v0=8600mm3 and w0=4500mm3, and the final state of y(t) is prescribed by yf=A=15mg/kg. For the initial guess, we choose a final time tf=2days, a control history u(τ)=0, and a slack variable history σ1=σ2=1 for τ[0,1].

The update formulas of the final time and the slack variables in Eq. (53) are scaled by α=0.001 and β=0.005. Moreover, we set the update parameters to ε1=0.2 and ε2=0.8, respectively. The presented algorithm is implemented in the software tool Matlab. The solution of the state and adjoint variables in Eqs. (48)(51) are obtained by the ode45 routine using the default solver settings for accuracy and variable step size control. The control input is discretized with 200 data points and due to the step size control, the numerical grids in the forward and backward solutions are not equivalent. Hence, the state and control trajectories are interpolated by piecewise polynomials using the modified Akima algorithm makima, which is suited to handle rapid changes between flat regions that occur during bang-bang controls. The solution to the tumor anti-angiogenesis problem is shown in the figures below. The convergence history is depicted in Fig. 4 and shows that after 150 iterations, the cost functional value is reduced from 8600mm3 to 7550mm3 and the final constraint converges to zero.

Fig. 4
Convergence of the optimization: Decreasing scrap function value in ratio to γ (top) and decreasing constraint value (bottom)
Fig. 4
Convergence of the optimization: Decreasing scrap function value in ratio to γ (top) and decreasing constraint value (bottom)
Close modal

The identified system inputs are presented in Fig. 5. Based on Pontryagin's minimum principle [28], the optimal control for bounded inputs, which appear linear in the Hamiltonian, is either a bang-bang or a singular control. In the tumor anti-angiogenesis problem, the angiogenic dose rate is of the bang-bang type. The optimal control is depicted in Fig. 5(a) and shows that the control input is either at its maximum or minimum value. The slack variables σ1(t) and σ2(t) associated with the control, which are used to account for the control limitations, are shown in Fig. 5(b). Here, σ1(t) and σ2(t) are associated with the upper and lower limits of the control.

Fig. 5
Control input and slack variables: (a) optimal angiogenesis inhibition and (b) slack variables
Fig. 5
Control input and slack variables: (a) optimal angiogenesis inhibition and (b) slack variables
Close modal

At the beginning, the control input u(t) remains in the upper constraint bound, resulting in a positive value of σ1(t). As the lower bound is not active, σ2(t) remains zero. When the control input switches from the upper to the lower constraint bound, we observe that σ1(t) changes to zero, while σ2(t)>0.

In Fig. 6(a), the optimal dosage decreases the volume of the primary tumor v(t) from 8600mm3 to approximately 7550mm3 after tf=1.2days. Figure 6(b) shows that the carrying capacity of the vasculature w(t) decreases during the administration of angiogenic agents and increases after discontinuation of therapy.

Fig. 6
State histories for the tumor anti-angiogenesis: (a) primary tumor volume, (b) carrying capacity of the vasculature, and (c) total amount of agents
Fig. 6
State histories for the tumor anti-angiogenesis: (a) primary tumor volume, (b) carrying capacity of the vasculature, and (c) total amount of agents
Close modal

Furthermore, in Fig. 6(c), the total amount of angiogenic agents administered over the whole period is monitored. While the maximum dose of drug is administered, the total amount of angiogenic agents increases linearly and remains constant after drug therapy is suspended.

The solution in Fig. 5(a) and in Figs. 6(a)6(c) is compared with the solution obtained by the software tool GPOPS-II presented in Ref. [29, Sec. 4.3.16]. GPOPS-II is a method used for solving nonlinear, continuous, and multiphase optimal control problems. It utilizes the Gauss pseudo-spectral method to transform the optimal control problem into a sparse, finite nonlinear programing problem. The results of the proposed gradient method shows good agreement for the states v, w and the control u with the presented results of GPOPS-II [29, Fig. 60]. Note that during the control change a discontinuity occurs both in the control and in the state derivative. Hence, in Ref. [29] the tumor anti-angiogenesis problem served as a test example to analyze convergence properties of various mesh refinement methods. The article describes that the GPOPS-II software tool is facing a large number of mesh points around the discontinuity. In contrast to the latter discussion in Ref. [29], a rapid change of control is not a problem in the presented method, as the number of control points can be chosen as required. The discontinuity of the carrying capacity w in Fig. 6(b) is handled by the automatic step size control of ode45.

6 Conclusion

In summary, the proposed method has several advantages over classical solution strategies using a penalty approach. It allows the use of initial controls that do not a priori satisfy the end conditions and inequality constraints. Moreover, no penalty terms are required, which provides more accurate solutions and makes it more robust in solving the underlying two-point boundary value problem. Despite the trivial choice of the initial control, the method converges to the optimal solution. Since the method enables an automatable process to find optimal solutions, its applicability covers even more challenging problems in optimal control. However, there are a few disadvantages to consider. Typical for gradient-based procedures, the method requires a high number of iterations near the optimum to achieve accurate results. Furthermore, the choice of process parameters, i.e., update and scaling parameters, demands experience and expertise.

Overall, the advantages of the proposed method, including flexibility in initial control, robustness, convergence, and applicability to complex systems, outweigh its drawbacks, thus making it a promising approach for solving optimal control problems in various domains, e.g., as shown here for optimal dose control in drug therapies. In many cases, it is difficult to derive the state equations of complex systems using a set of minimal independent coordinates. Therefore, many multibody systems are modeled by commercial software tools dealing with index-three differential algebraic equations (DAEs). A modified gradient approach tailored for multibody systems, however without applicability to inequality constraints, has been investigated in a previous work [30]. The prework is also based on the adjoint method including a backward differentiation formula scheme to solve the resulting adjoint DAE system. Building on the prework along with the knowledge acquired in this article, the theory can be extended for DAE systems to handle control problems with inequality constraints in multibody dynamics.

Acknowledgment

This research was funded in whole or in part by the Austrian Science Fund (FWF): P 37110-NBL (No. 10.55776/P37110) Karin Nachbagauer acknowledges also support from the Technical University of Munich – Institute for Advanced Study.

Conflict of Interest

The authors have no competing interests to declare that are relevant to the content of this article.

References

1.
Johnson
,
C.
, and
Gibson
,
J.
,
1963
, “
Singular Solutions in Problems of Optimal Control
,”
IEEE Trans. Autom. Control
,
8
(
1
), pp.
4
15
.10.1109/TAC.1963.1105505
2.
Kelley
,
H. J.
,
1962
, “
Method of Gradients – Optimization Techniques With Applications to Aerospace Systems
,”
Mathematics in Science and Engineering
,
Academic Press
, New York, Vol.
5
, pp.
206
254
.10.1016/S0076 5392(08)620949
3.
Bryson
,
A.
, and
Ho
,
Y.
,
1975
,
Applied Optimal Control: Optimization, Estimation and Control
,
Hemisphere
,
Washington, DC
.
4.
Bertolazzi
,
E.
,
Biral
,
F.
, and
Lio
,
M. D.
,
2005
, “
Symbolic-Numeric Indirect Method for Solving Optimal Control Problems for Large Multibody Systems
,”
Multibody Syst. Dyn.
,
13
(
2
), pp.
233
252
.10.1007/s11044-005-3987-4
5.
Bottasso
,
C. L.
,
Croce
,
A.
,
Ghezzi
,
L.
, and
Faure
,
P.
,
2004
, “
On the Solution of Inverse Dynamics and Trajectory Optimization Problems for Multibody Systems
,”
Multibody Syst. Dyn.
,
11
(
1
), pp.
1
22
.10.1023/B:MUBO.0000014875.66058.74
6.
Sargent
,
R. W. H.
,
2000
, “
Optimal Control
,”
J. Comput. Appl. Math.
,
124
(
1–2
), pp.
361
371
.10.1016/S0377-0427(00)00418-0
7.
Jarrett
,
A. M.
,
Faghihi
,
D.
,
Hormuth
,
D. A.
,
Lima
,
E. A.
,
Virostko
,
J.
,
Biros
,
G.
,
Patt
,
D.
, and
Yankeelov
,
T. E.
,
2020
, “
Optimal Control Theory for Personalized Therapeutic Regimens in Oncology: Background, History, Challenges, and Opportunities
,”
J. Clin. Med.
,
9
(
5
), p.
1314
.10.3390/jcm9051314
8.
Hahnfeldt
,
P.
,
Panigrahy
,
D.
,
Folkman
,
J.
, and
Hlatky
,
L.
,
1999
, “
Tumor Development Under Angiogenic Signaling: A Dynamical Theory of Tumor Growth, Treatment Response, and Postvascular Dormancy
,”
Cancer Res.
,
59
(
19
), pp.
4770
4775
.http://epbiradivot.cwru.edu/EPBI473/files/wk13tumorTherapy/philCanRes99model.pdf
9.
d'Onofrio
,
A.
, and
Gandolfi
,
A.
,
2004
, “
Tumour Eradication by Antiangiogenic Therapy: Analysis and Extensions of the Model by Hahnfeldt et al. (1999)
,”
Math. Biosci.
,
191
(
2
), pp.
159
184
.10.1016/j.mbs.2004.06.003
10.
Ergun
,
A.
,
Camphausen
,
K.
, and
Wein
,
L. M.
,
2003
, “
Optimal Scheduling of Radiotherapy and Angiogenic Inhibitors
,”
Bull. Math. Biol.
,
65
(
3
), pp.
407
424
.10.1016/S0092-8240(03)00006-5
11.
Glick
,
A. E.
, and
Mastroberardino
,
A.
,
2017
, “
An Optimal Control Approach for the Treatment of Solid Tumors With Angiogenesis Inhibitors
,”
Mathematics
,
5
(
4
), p.
49
.10.3390/math5040049
12.
Ledzewicz
,
U.
, and
Schättler
,
H.
,
2008
, “
Analysis of Optimal Controls for a Mathematical Model of Tumour Anti-Angiogenesis
,”
Optim. Control Appl. Methods
,
29
(
1
), pp.
41
57
.10.1002/oca.814
13.
Gao
,
D.
, and
Huang
,
N.
,
2018
, “
Optimal Control Analysis of a Tuberculosis Model
,”
Appl. Math. Modell.
,
58
, pp.
47
64
.10.1016/j.apm.2017.12.027
14.
Hametner
,
C.
,
Kozek
,
M.
,
Bohler
,
L.
,
Wasserburger
,
A.
,
Du
,
Z. P.
,
Kolbl
,
R.
,
Bergmann
,
M.
,
Bachleitner-Hofmann
,
T.
, and
Jakubek
,
S.
,
2021
, “
Estimation of Exogenous Drivers to Predict Covid-19 Pandemic Using a Method From Nonlinear Control Theory
,”
Nonlinear Dyn.
,
106
(
1
), pp.
1111
1125
.10.1007/s11071-021-06811-7
15.
Graichen
,
K.
,
Kugi
,
A.
,
Petit
,
N.
, and
Chaplais
,
F.
,
2010
, “
Handling Constraints in Optimal Control With Saturation Functions and System Extension
,”
Syst. Control Lett.
,
59
(
11
), pp.
671
679
.10.1016/j.sysconle.2010.08.003
16.
Haug
,
E. J.
,
Wehage
,
R. A.
, and
Mani
,
N. K.
,
1984
, “
Design Sensitivity Analysis of Large-Scale Constrained Dynamic Mechanical Systems
,”
ASME J. Mech., Trans., Autom. Des.
,
106
(
2
), pp.
156
162
.10.1115/1.3258573
17.
Bestle
,
D.
, and
Eberhard
,
P.
,
1992
, “
Analyzing and Optimizing Multibody Systems
,”
Mech. Struct. Mach.
,
20
(
1
), pp.
67
92
.10.1080/08905459208905161
18.
Cao
,
Y.
,
Li
,
S.
,
Petzold
,
L.
, and
Serban
,
R.
,
2003
, “
Adjoint Sensitivity Analysis for Differential-Algebraic Equations: The Adjoint DAE System and Its Numerical Solution
,”
SIAM J. Sci. Comput.
,
24
(
3
), pp.
1076
1089
.10.1137/S1064827501380630
19.
Steiner
,
W.
, and
Reichl
,
S.
,
2012
, “
The Optimal Control Approach to Dynamical Inverse Problems
,”
ASME J. Dyn. Syst., Meas., Control
,
134
(
2
), p.
021010
.10.1115/1.4005365
20.
Nachbagauer
,
K.
,
Oberpeilsteiner
,
S.
,
Sherif
,
K.
, and
Steiner
,
W.
,
2015
, “
The Use of the Adjoint Method for Solving Typical Optimization Problems in Multibody Dynamics
,”
ASME J. Comput. Nonlinear Dyn.
,
10
(
6
), p.
061011
.10.1115/1.4028417
21.
Eichmeir
,
P.
,
Nachbagauer
,
K.
, and
Steiner
,
W.
,
2020
, “
The Adjoint Gradient Method for Time-Optimal Control of a Moon Landing – Ascent, Descent, and Abort
,”
ASME
Paper No. DETC2020-22034.10.1115/DETC2020 22034
22.
Eichmeir
,
P.
,
Nachbagauer
,
K.
,
Lauß
,
T.
,
Sherif
,
K.
, and
Steiner
,
W.
,
2021
, “
Time-Optimal Control of Dynamic Systems Regarding Final Constraints
,”
ASME J. Comput. Nonlinear Dyn.
,
16
(
3
), p.
031003
.10.1115/1.4049334
23.
Jacobson
,
D.
, and
Lele
,
M.
,
1969
, “
A Transformation Technique for Optimal Control Problems With a State Variable Inequality Constraint
,”
IEEE Trans. Autom. Control
,
14
(
5
), pp.
457
464
.10.1109/TAC.1969.1099283
24.
Feehery
,
W. F.
,
1998
, “
Dynamic Optimization With Path Constraints
,”
Ph.D. thesis
,
Massachusetts Institute of Technology
, Cambridge, MA.http://hdl.handle.net/1721.1/29879
25.
Kirk
,
D.
,
2004
,
Optimal Control Theory: An Introduction
,
Dover Publications
,
New York
.
26.
Bhatti
,
M. A.
,
2012
,
Practical Optimization Methods: With Mathematica® Applications
,
Springer Science & Business Media
,
New York
.
27.
Luenberger
,
D. G.
, and
Ye
,
Y.
,
2016
,
Linear and Nonlinear Programming
,
Springer
,
New York
.
28.
Pontryagin
,
L.
,
Boltyanskii
,
V.
,
Gamkrelidze
,
R.
, and
Mischchenko
,
E.
,
1962
, “
The Mathematical Theory of Optimal Processes
,”
Wiley Interscience
,
New York
.
29.
Rao
,
A.
, and
Hager
,
W.
,
2018
, “
Solution of Optimal Control Problem for High Speed Ascent and Reentry Vehicles
,”
Florida University
,
Gainesville
, Report No. AFRL-RW-EG-TR-2018-058.https://apps.dtic.mil/sti/trecms/pdf/AD1063085.pdf
30.
Eichmeir
,
P.
,
2022
, “
The Adjoint Gradient Method for Time-Optimal Control of Multibody Systems
,”
Ph.D. thesis
,
Vienna University of Technology
,
Vienna, Austria
.10.34726/hss.2022.98567