The Use of Slack Variables in the Adjoint Method Handling Inequality Constraints in Optimal Control and the Application to Tumor Drug Dosage

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 [1–3]. 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 [4–6].

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 [8–12].

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., [16–20], 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.

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)2−1$⁠. 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.

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$⁠.

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.

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.

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

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

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 $δΓ$⁠.

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.

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

   $x′(τ)=tff(x(τ),u(τ))$

   (48)

   with the initial condition $x(0)=x0$ in the interval $τ∈[0,1]$⁠.
2. Solve the adjoint differential equations

   $p′(τ)=−tf(hxT(x(τ),u(τ))+fxT(x(τ),u(τ))p(τ))$

   (49)

   with the final condition $p(1)=SxT(x(1),1)$ from $τ=1$ to $τ=0$⁠.
3. Solve the influence differential equations I

   $P′(τ)=−tffxT(x(τ),u(τ))P(τ)$

   (50)

   with the final condition $P(1)=ϕxT(x(1),1)$ from $τ=1$ to $τ=0$⁠.
4. Solve the influence differential equations II

   $Q′(τ)=−tf(gxT(x(τ),u(τ),σ(τ))+fxT(x(τ),u(τ))Q(τ))$

   (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

   $b=A(νμ.)$

   (52)
6. Select $κ$ and update the final time, controls and slack variables by

   $δtf=−κα(hf+S˙f+νTϕ˙f+μTgf)δu(τ)=−κ(huT+fuTp+fuTPν+guTμ+fuTQμ)δσ(τ)=−κβgσTμ$

   (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.

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.

in terms of the proposed algorithm.

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

   $v′(τ)=−tfζv(τ)ln(v(τ)w(τ))w′(τ)=tfw(τ)(b−μ−dv23(τ)−Gu(τ))y′(τ)=tfu(τ)$

   (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

   $p1′(τ)=2dw(τ)3v13(τ)tfp2(τ)+(1+ln(v(τ)w(τ)))ζtfp1(τ)p2′(τ)=−(b+μ+Gu(τ)+dv23(τ))tfp2(τ)−v(τ)w(τ)ζtfp1(τ)p3′(τ)=0$

   (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

   $P11′(τ)=2dw(τ)3v13(τ)tfP21(τ)+(1+ln(v(τ)w(τ)))ζtfP11(τ)P21′(τ)=−(b+μ+Gu(τ)+dv23(τ))tfP21(τ)−v(τ)w(τ)ζtfP11(τ)P31′(τ)=0$

   (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)₁ read

   $M1,11=∫t0tf(u(t)−σ12(t))(P31(t)−GP21(t)w(t))dt+α12u(tf)(σ12(tf)−u(tf))2(y(tf)−A)M1,12=∫t0tf(σ22(t)+u(t)−u¯)(P31(t)−GP21(t)w(t))dt+α12u(tf)(σ22(tf)+u(tf)−u¯)2(y(tf)−A)$

   (63)

   and $N1∈R$ can be determined by Eq. (41)₂, yielding

   $N1=∫t0tf((P31(t)−GP21(t)w(t))2)dt+αu2(tf)(A−y(tf))2$

   (64)

   Moreover, $b1∈R$ can be computed by Eq. (41)₃ and is given by

   $b1=∫t0tf(p3(t)−Gp2(t)w(t))(P31(t)−GP21(t)w(t))dt$

   (65)

   where the boundary term is zero, since $hf=0$⁠. The components of the matrix $M2=(M2,ij)∈R2×2$ in Eq. (45)₁ read

   $M2,11=∫t0tf(σ12(t)−u(t))2(1+4βσ12(t))dt+α14(σ12(tf)−u(tf))4M2,22=∫t0tf(σ22(t)+u(t)−u¯)2(1+4βσ22(t))dt+α14(σ22(t)+u(t)−u¯)4M2,12=∫t0tf(σ22+u(t)−u¯)(u(t)−σ12(t))dt+α14(σ22(t)+u(t)−u¯)2(σ12(t)−u(t))2$

   (66)

   where $M2,12=M2,21$⁠. Moreover, $N2∈R2$ can be determined by Eq. (45)₂

   $N2,1=∫t0tf(u(t)−σ12(t))(P31(t)−GP21(t)w(t))dt+α12u(tf)(u(tf)−σ12(tf))2(A−y(tf))N2,2=∫t0tf(σ22+u(t)−u¯)(P31(t)−GP21(t)w(t))dt+α12u(tf)(σ22+u(tf)−u¯)2(y(tf)−A)$

   (67)

   Inserting in Eq. (45)₃, $b2∈R2$ reads

   $b2,1=∫t0tf(u(t)−σ12(t))(p3(t)−Gp2(t)w(t))dtb2,2=∫t0tf(σ22(t)+u(t)−u¯)(p3(t)−Gp2(t)w(t))dt$

   (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

   $δtf=−κα12[μ1(σ12−u(tf))2+μ2(σ22+u(tf)−u¯)2+2νu(tf)(y(tf)−A)−2γv(tf)ζln(v(tf)w(tf))]δu(τ)=−κ[p3(τ)+12μ1(σ12(τ)−u(τ))2+12μ2(σ22(τ)+u(τ)−u¯)2−Gp2(τ)w(τ)+ν1(P31(τ)−GP21(τ)w(τ))]$

   (69)

   and the updates for the slack variables read

   $δσ1(τ)=−2κβμ1σ1(τ)(σ12(τ)−u(τ))δσ2(τ)=−2κβμ2σ2(τ)(σ22(τ)+u(τ)−u¯)$

   (70)