Implicit Third-Order Peer Triplets with Variable Stepsizes for Gradient-Based Solutions in Large-Scale ODE-Constrained Optimal Control

TL;DR

This work advances gradient-based ODE-constrained optimal control by introducing two third-order implicit Peer two-step triplets, AP4o33vgi and AP4o33vsi, that preserve high-order accuracy on variable time grids. The methods feature symmetry properties, boundary-step designs, and triangular boundary iterations to enable efficient solutions for large-scale problems, supported by a posteriori error estimators that drive adaptive equi-distribution of global errors. Theoretical results establish global error bounds and super-convergence under appropriate grid conditions, while numerical experiments on a 1D heat boundary-control problem and a 2D prostate cancer growth model demonstrate substantial gains in accuracy and computational efficiency. The approach significantly enhances the practical applicability of high-order Peer methods to complex PDE-constrained control tasks, including multi-dose drug protocol design in oncology.

Abstract

This paper is concerned with the theory, construction and application of variable-stepsize implicit Peer two-step methods that are super-convergent for variable stepsizes, i.e., preserve their classical order achieved for uniform stepsizes when applied in a gradient-based solution algorithm to solve ODE-constrained optimal control problems in a first-discretize-then-optimize setting. Gradients of the objective function can be computed most efficiently using approximate adjoint variables. High accuracy with moderate computational effort can be achieved through time integration methods that satisfy a sufficiently large number of adjoint order conditions for variable stepsizes and provide gradients with higher-order consistency. In this paper, we enhance our previously developed variable implicit two-step Peer triplets constructed in [J. Comput. Appl. Math. 460, 2025] to get ready for large-scale dynamical systems with varying time scales without losing efficiency. A key advantage of Peer methods is their use of multiple stages with the same high stage order, which prevents order reduction - an issue commonly encountered in semi discretized PDE problems with boundary control. Two third-order methods with four stages, good stability properties, small error constants, and a grid adaptation by equi-distributing global errors are constructed and tested for a 1D boundary heat control problem and an optimal control of cytotoxic therapies in the treatment of prostate cancer.

Figures (3)

• Figure 1: Dirichlet heat problem with $m=250$ spatial points. Exemplary mesh density function $\psi(t)$ for AP4o33vgi, $N=63$ and adapted time grids for $N=15,31,63,123$ (top left). Convergence of the maximal control errors $\|U_{ni}-u(t_{ni})\|_\infty$ (top right), state errors $\|y(T)$$-$$y_h(T)\|_\infty$ (bottom left), and adjoint errors $\|p(0)$$-$$p_h(0)\|_\infty$ (bottom right) for uniform and adaptive time grids.
• Figure 2: AP4o33vgi for PCa problem: single drug cycle. The target optimal cytotoxic 1-dose drug $U_{d1}(t)$ (bottom left) yields a significant smaller tumor volume $V_{\phi,d1}$ (top right) and serum PSA $P_{s,d1}$ (bottom right), compared to $V_{\phi,0}$ and $P_{s,0}$ obtained for the standard docetaxel protocol $U_0(t)$. Refining the time grid at both ends of the interval using the equi-distribution principle reduces the mesh density function there (top left) and leads to a reduction of the estimated maximum global error by nearly 40%. The measures of the adaptive mesh are $\sigma_n\in [0.75,1.26]$ and $|\eta_n|\le 2$.
• Figure 3: AP4o33vgi for PCa problem: 3-dose drug cycle. The numerical target optimal cytotoxic drug $U_\star(t)$ (bottom left) yields a reduced tumor volume $V_{\phi,\star}$ (top right) and serum PSA $P_{s,\star}$ (bottom right), compared to $V_{\phi,d3}$ and $P_{s,d3}$ obtained for the 3-dose docetaxel protocol $U_{d3}(t)$. Refining the time grid at both ends of the interval and at the drug delivery time points $t_c$ using the equi-distribution principle reduces the mesh density function in critical regions (top left) and leads to a reduction of the estimated maximum global error by 38%. The measures of the adaptive mesh are $\sigma_n\in [0.75,1.20]$ and $|\eta_n|\le 2$.

Theorems & Definitions (6)

• Lemma 5.1
• Remark 5.1
• Remark 5.2
• Theorem 6.1
• Lemma 7.1
• Remark 8.1