ParaOpt for unstable systems

Djahou N Tognon

► To cite this version:

Djahou N Tognon. ParaOpt for unstable systems. 2025. hal-04966081

HAL Id: hal-04966081

https://inria.hal.science/hal-04966081v1

Preprint submitted on 25 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Djahou N. Tognon

Sorbonne Université, Université Paris Cité, CNRS, INRIA,
Laboratoire Jacques-Louis Lions, LJLL, EPC ANGE,
F-75005 Paris, France
djahou-norbert.tognon@inria.fr

ABSTRACT

ParaOpt is a two-level time-parallel method to solve the coupled forward/backward Euler-Lagrange system arising from partial differential equations (PDEs) constrained optimization. In this work, we present a convergence analysis of this algorithm in the case where the system under consideration is unstable. We complete this theoretical study with numerical experiments, where the properties of the algorithm are investigated on linear and nonlinear examples.

**Keywords** ParaOpt · Optimal Control · Unstable Systems

Introduction

Time parallel time integration for optimal control problems became in recent decades an active topic of research to solve very large-scale problems. The multicore architectures and graphics processing units that appeared due to the limited processing speed in current circuit technology pose great challenges to accelerate classical numerical algorithms. Moreover, optimal control problems which are typically addressed by solving their optimality system. Such systems generally comprise a combination of a direct equation and an adjoint equation, along with an equation that specifies the optimality of the control [24]. The discretizations in time and space of this type of problem lead to very large-scale problems. To solve such problems in parallel, various approaches have been proposed. All are based on different parallel in time algorithms [4], such as multiple shooting methods [19, 10] (and references therein), direct method [8], Schwarz method [7] and multigrid methods [2]. A simple way to combine time parallel and optimal control solvers consists of applying the former to the forward and backward problems appearing in the latter [11, 12, 21, 13]. In a more complex approach, one recasts the parallel solver methodology to solve the optimality system [1, 9, 18, 5]. This approach is followed in this paper. Finally, some properties of the system can be used to design specific methods. In the context of quantum control, the decomposition of the original control problem into smaller ones is considered in [20, 23]. Here, intermediate variables on time sub-intervals are defined at the sub-interval boundaries.

In ParaOpt [5], which is an algorithm based on a combined solving of forward and backward trajectories. This algorithm is based on the fundamental idea of the Parareal algorithm [19], hence it is related to the multiple shooting paradigm [22, 10], where the time interval is decomposed into a set of non-overlapping subintervals. Here, the solution procedure considers as unknowns the interface state and adjoint variables that are computed using an inexact Newton method so that trajectories are continuous on the boundaries. The main innovation of this approach consists in applying the Parareal methodology [19] to a two-point boundary value problem. In this way, a coarse solver is used to approximate the Jacobian in the Newton iteration.

Let us briefly introduce the optimal control framework that we are interested in in this work. We consider the linear optimal control problem associated with the cost functional

$$ \mathcal{J}(\nu) := \frac{1}{2} \|y(T) - y_{\text{target}}\|^2 + \frac{\alpha}{2} \int_0^T \|\nu(t)\|^2 dt, \quad (1) $$
where $\alpha$ is a fixed regularization parameter, $y_{\text{target}}$ is a target state. The evolution of the state function $y: [0, T] \to \mathbb{R}^r$ is described by the ordinary differential equation

$$ 
\begin{aligned}
\dot{y}(t) &= A y(t) + \nu(t), \\
y(0) &= y_{\text{init}}
\end{aligned}
\qquad (2) $$

where $\nu(t)$ is the control, and $A$ is real and symmetric operator independent on time variable. The spectral decomposition of $A$ allows to diagonalize the matrix and reduce the problem to a family of Dahlquist problems (see [14], Chap IV) of the form

$$ \dot{y}(t) = \sigma y(t) + \nu(t), \qquad (3) $$

where $y$ and $\nu$ are now scalars and $\sigma$ is an eigenvalue of $A$.

In [5], ParaOpt is applied to the discrete optimality system arising from (1)-(3) and takes the form of a linear iterative procedure where the coarse approximation of the Jacobian acts as an integration preconditioner. In that work, a spectral analysis of the iteration matrix is provided for the case where $\sigma$ is negative in (3), corresponding to dissipative or damped systems. This analysis is based on the argument principle of complex analysis. In this paper, we investigate the convergence properties of the algorithm in the case of positive $\sigma$, which typically occurs when considering unstable systems. Our analysis is based on Rouché's theorem (see [3]), where the spectral radius of the iteration matrix is bounded by an explicit constant depending on the parameters of the discrete optimization problem.

The paper is structured as follows. The discrete formulation of the optimal control problem is presented in Section 1. In Section 2, we present the linear formulation of the ParaOpt in our case. Section 3 is dedicated to the convergence analysis of the ParaOpt. We complete these results with some numerical experiments in Section 4.

# 1 Discrete Formulation

Let $\sigma$ be a real positive number. We consider the optimal control problem in dimension 1 which consists in minimizing the cost functional (1) subject to the constraint (3).

Instead of using the implicit Euler method as in [5] in the case of negative eigenvalue $\sigma$, here we use the explicit Euler method to discretize the Dahlquist test problem (3). The reason for this is detailed in Remark 2.1 below. In this way, given $M \in \mathbb{N}$, the time interval $[0, T]$ is discretized by a grid $t_0 = 0, \dots, t_M = T$, where $t_n = n\delta t$ and $\delta t = T/M$. The explicit Euler method applied to (2) then reads

$$ y_{n+1} = (1 + \sigma \delta t)y_n + \delta t \nu_n, \quad y_0 = y_{\text{init}}, $$

where $y_n$ and $\nu_n$ stand for the approximations of the state and the control at time $t_n$ respectively.

In this framework, we consider the discrete version of $\mathcal{J}$ (see (1)) defined by

$$ \mathcal{J}_{\delta t}(\nu) := \frac{1}{2} \|y_M - y_{\text{target}}\|^2 + \frac{\alpha}{2} \delta t \sum_{n=0}^{M-1} \|\nu_n\|^2. $$

The Lagrange operator (see [6, 15] for details) associated with this problem is given by

$$ \mathcal{L}_{\delta t}(y, \nu, \lambda) := \mathcal{J}_{\delta t}(\nu) - \sum_{n=0}^{M-1} \langle \lambda_{n+1}, y_{n+1} - (1 + \sigma \delta t)y_n - \delta t \nu_n \rangle, $$

where $\lambda$ is the adjoint-state function and $\lambda_n$ its approximations on $t_n$. Therefore, the Euler-Lagrange equations lead to the following optimality system.

$$ 
\begin{aligned}
y_0 &= y_{\text{init}}, \\
y_{n+1} &= (1 + \sigma \delta t)y_n + \delta t \nu_n, \quad n = 0, \dots, M-1, \\
\lambda_n &= (1 + \sigma \delta t)\lambda_{n+1}, \quad n = 1, \dots, M-1, \\
\alpha \nu_n &= -\lambda_{n+1}, \quad n = 0, \dots, M-1, \\
\lambda_M &= y_M - y_{\text{target}}.
\end{aligned}
 $$

Eliminating $\nu$ in this system gives

$$ y_{n+1} = (1 + \sigma \delta t)y_n - \frac{\delta t}{\alpha} \lambda_{n+1}, \quad n = 0, \dots, M-1, \qquad (4) $$

$$ \lambda_n = (1 + \sigma \delta t)\lambda_{n+1}, \quad n = 1, \dots, M-1, \qquad (5) $$

with $y_0 = y_{\text{init}}$ and $\lambda_M = y_M - y_{\text{target}}$.
## 2 Time-parallelization

In this section, we present the formulation of ParaOpt to numerically solve the optimality system (4)-(5). We consider a partition of $[0, T]$ into $L$ subintervals of uniform length $\Delta T$ with $T_\ell = \ell \Delta T$ such that $\Delta T = N \delta t$ and $M = NL$ (see Figure 1). We also consider two sets of intermediate states $(Y_\ell)_{\ell=0,...,L}$ and $(\Lambda_\ell)_{\ell=1,...,L}$ corresponding to the approximations of state $y$ and adjoint state $\lambda$ at time $T_\ell$, $\ell = 0, ..., L$ and $\ell = 1, ..., L$, respectively. We then express the discrete optimality system (4)-(5) in terms of the variables $Y_\ell$ and $\Lambda_\ell$ located in $T_\ell$. In other words, we eliminate the unknowns related to the interior of each subinterval $[T_\ell, T_{\ell+1}]$. To this end, let $n_1$ and $n_2$ be two integers such that $0 \le n_1 \le n_2 \le M$. From (4), we get

$$y_{n_2} = (1 + \sigma \delta t)^{n_2 - n_1} y_{n_1} - \frac{\delta t}{\alpha} \sum_{j=0}^{n_2 - n_1 - 1} (1 + \sigma \delta t)^j \lambda_{n_2 - j}. \quad (6)$$

Also, (5) implies that

$$\lambda_{n_2-n_1} = (1 + \sigma \delta t)^{n_2-n_1} \lambda_{n_2}. \quad (7)$$

Now, let $j < n_2$ and set $n_1 = n_2 - j$. Because of (7), we have $\lambda_{n_2-j} = (1 + \sigma \delta t)^j \lambda_{n_2}$, so that substituting $\lambda_{n_2-j}$ into (6), we obtain

$$y_{n_2} = (1 + \sigma \delta t)^{n_2 - n_1} y_{n_1} - \frac{\delta t}{\alpha} \left[ \sum_{j=0}^{n_2 - n_1 - 1} (1 + \sigma \delta t)^{2j} \right] \lambda_{n_2}. \quad (8)$$

Setting $n_1 = (\ell - 1)N$ and $n_2 = \ell N$ and using (8) and (7), the reduced optimality system (4)-(5) reads

$$\begin{aligned} Y_{\ell} - \beta_{\delta t} Y_{\ell-1} + \frac{\gamma_{\delta t}}{\alpha} \Lambda_{\ell} &= 0, & \ell &= 1, \dots, L; \\ \Lambda_{\ell-1} - \beta_{\delta t} \Lambda_{\ell} &= 0, & \ell &= 2, \dots, L, \end{aligned} \quad (9)$$

with $Y_0 = y_{init}$ and $\Lambda_L - Y_L = -y_{target}$, where

$$\beta_{\delta t} := (1 + \sigma \delta t)^{\frac{\Delta T}{\delta t}}, \quad (10)$$

$$\gamma_{\delta t} := \delta t \left[ \sum_{j=0}^{N-1} (1 + \sigma \delta t)^{2j} \right] = \frac{\beta_{\delta t}^2 - 1}{\sigma(2 + \sigma \delta t)}. \quad (11)$$

Then collecting all the unknowns $Y_\ell$ and $\Lambda_\ell$ in the vector $X^T = (Y_0, Y_1, \dots, Y_L | \Lambda_1, \dots, \Lambda_L)$, (9) yields

$$\mathcal{F}_{\delta t}(X) := \mathcal{M}_{\delta t}X - \mathfrak{b} = 0, \quad (12)$$

where $\mathfrak{b}^T = (y_{init}, 0, \dots, 0, \dots, 0, -y_{target})$ and

$$\mathcal{M}_{\delta t} := \left( \begin{array}{c|c} \begin{array}{cccc} 1 & & & \\ & -\beta_{\delta t} & 1 & \\ & & \ddots & \\ & & & -\beta_{\delta t} \end{array} & \begin{array}{c} 0 \\ \gamma_{\delta t}/\alpha \\ \vdots \\ 0 \end{array} \\ \hline \begin{array}{c} \vdots \\ -1 \end{array} & \begin{array}{c} 1 & -\beta_{\delta t} \\ & \ddots & \ddots \\ & & 1 & -\beta_{\delta t} \\ & & & 1 \end{array} \end{array} \right).$$

We consider another time step $\Delta t$ such that $\delta t \le \Delta t \le \Delta T$, see Figure 1. In what follows, $\Delta t$ and $\delta t$ play the role of coarse and fine time steps, that are usually considered in the Parareal algorithm [19]. The values of these discretization parameters are the same in each subinterval $[T_\ell, T_{\ell+1}]$. To solve numerically (12), ParaOpt (see [5]) consists of the following Newton-type iteration

$$\mathcal{F}'_{\Delta t}(X^i)(X^{i+1} - X^i) = -\mathcal{F}_{\delta t}(X^i),$$

where $\mathcal{F}'$ denotes the derivative of $\mathcal{F}$ with respect to its variable $X$. In the linear framework we are considering, this equation is equivalent to

$$\mathcal{M}_{\Delta t}(X^{i+1} - X^i) = -(\mathcal{M}_{\delta t}X^i - \mathfrak{b}). \quad (13)$$

In the next section, we provide the convergence analysis of (13) which is based on the eigenvalues of the iteration matrix $(I - \mathcal{M}_{\Delta t}^{-1}\mathcal{M}_{\delta t})$.
<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic representation of the decomposition of $[0, T]$ into sub-intervals of the form $[T_\ell, T_{\ell+1}]$. Each subinterval $[T_\ell, T_{\ell+1}]$ is discretized using a coarse and a fine grid. In this example, we have taken $\Delta t = 2\delta t$.</figcaption></figure>

**Remark 2.1.** Similar analysis is provided in [5], where an implicit Euler method is used to discretize (3). In our case for large values of $\sigma$, $\beta_{\delta t} = (1 - \sigma\delta t)^{\frac{\Delta T}{\delta t}}$ which makes sense only if $\Delta T < 1/\sigma$, or equivalently $L > \sigma T$ when $\sigma$ is a real positive number. As a consequence, $L$ must be scaled to $\sigma$, making this approach numerically inefficient in the case of large $\sigma$.

### 3 Eigenvalues Analysis

We present in this section the convergence analysis of ParaOpt given by the linear iteration (13). In this way, we denote by $\mu$ an eigenvalue of the iteration matrix $I - \mathcal{M}_{\Delta t}^{-1}\mathcal{M}_{\delta t}$. Let $\hat{x}$ be the eigenvector associated with the eigenvalue $\mu$. Then we have

$$ (\mathcal{M}_{\Delta t} - \mathcal{M}_{\delta t}) \hat{x} = \mu \mathcal{M}_{\Delta t} \hat{x}. \quad (14) $$

To fit with the bloc structure of $I - \mathcal{M}_{\Delta t}^{-1}\mathcal{M}_{\delta t}$, we write $\hat{x}$ in the form

$$ \hat{x} = (v_0, v_1, \dots, v_L | w_1, \dots, w_L)^T, $$

so that $v_\ell$ and $w_\ell$ are the intermediate state and adjoint state associated with the sub-interval $[T_\ell, T_{\ell+1}]$. Since $\mathcal{M}_{\Delta t} - \mathcal{M}_{\delta t}$ has two zero rows, the eigenvalue $\mu = 0$ must have a multiplicity of at least two. Now, let $\mu$ be a non-zero eigenvalue. Recall that if no such eigenvalue exists, the iteration matrix is nilpotent, and the iteration converges in a finite number of steps. Writing (14) componentwise yields

$$ \begin{aligned} \mu v_0 &= 0, \\ \mu(-\beta_{\Delta t} v_{\ell-1} + v_\ell + \frac{1}{\alpha}\gamma_{\Delta t} w_\ell) &= (\beta_{\delta t} - \beta_{\Delta t}) v_{\ell-1} + \frac{1}{\alpha}(\gamma_{\Delta t} - \gamma_{\delta t}) w_\ell, \quad \ell = 1, \dots, L \\ \mu(w_{\ell-1} - \beta_{\Delta t} w_\ell) &= (\beta_{\delta t} - \beta_{\Delta t}) w_\ell, \quad \ell = 2, \dots, L \\ \mu(w_L - v_L) &= 0. \end{aligned} $$

To simplify the notation, we introduce $\beta = \beta_{\Delta t}$, $\gamma = \gamma_{\Delta t}$, $\delta\beta = \beta - \beta_{\delta t}$ and $\delta\gamma = \gamma - \gamma_{\delta t}$, which gives

$$ \begin{aligned} v_0 &= 0 \\ \mu(-\beta v_{\ell-1} + v_\ell + \frac{1}{\alpha}\gamma w_\ell) &= -\delta\beta v_{\ell-1} + \frac{1}{\alpha}\delta\gamma w_\ell, \quad \ell = 1, \dots, L \end{aligned} \quad (15) $$

$$ \mu(w_{\ell-1} - \beta w_\ell) = -\delta\beta w_\ell, \quad \ell = 2, \dots, L \quad (16) $$

$$ w_L - v_L = 0. \quad (17) $$

Using (15) and (16), we get

$$ \begin{aligned} v_\ell &= av_{\ell-1} + bw_\ell, \quad \ell = 1, \dots, L \\ w_{\ell-1} &= aw_\ell, \quad \ell = 2, \dots, L, \end{aligned} $$

where

$$ a = \beta - \frac{\delta\beta}{\mu}, \quad b = \frac{1}{\alpha} \left[ \frac{\delta\gamma}{\mu} - \gamma \right], \quad (18) $$

so that

$$ w_\ell = a^{L-\ell} w_L \quad \text{and} \quad v_\ell = \sum_{j=0}^{\ell} a^j b w_{\ell-j}. \quad (19) $$
Combining the two expressions, we get

$$v_L = b \left[ \sum_{\ell=0}^{L-1} a^{2\ell} \right] w_L. \qquad (20)$$

Using (17) and (20), we obtain

$$\left( 1 - b \sum_{\ell=0}^{L-1} a^{2\ell} \right) w_L = 0.$$

The component $w_L$ of the eigenvector $\hat{x}$ cannot vanish. In fact, suppose that $w_L = 0$. According to (19) and (20), we find that $\hat{x} = 0$, which contradicts the fact that $\hat{x}$ is an eigenvector. Therefore,

$$1 - b \sum_{\ell=0}^{L-1} a^{2\ell} = 0. \qquad (21)$$

Using the definition of $b$ in (18), we obtain $\mu = \delta\gamma/(\alpha b + \gamma)$. Introducing

$$\tau := \beta - \gamma \frac{\delta\beta}{\delta\gamma}, \qquad (22)$$

and combining it with (18), we get $b = (\tau - a)\delta\gamma/\alpha\delta\beta$. Therefore, (21) becomes

$$\mathcal{P}(a) := \alpha \frac{\delta\beta}{\delta\gamma} - (\tau - a) \sum_{\ell=0}^{L-1} a^{2\ell} = 0. \qquad (23)$$

Each eigenvalue $\mu$ of $I - M_{\Delta t}^{-1}M_{\delta t}$ is associated with a root of $\mathcal{P}$ through the relation.

$$\mu = \frac{\delta\beta}{\beta - a}. \qquad (24)$$

Each root $a$ depends on the parameters $\beta$, $\gamma$, $\beta_{\delta t}$, $\gamma_{\delta t}$ and $\alpha$. We now analyze the behavior of the roots of $\mathcal{P}$ to get conditions for the convergence of the iterate process associated with $I - M_{\Delta t}^{-1}M_{\delta t}$. We start with some properties of those parameters.

**Lemma 3.1.** Let $\sigma > 0$. Then we have $\delta\beta < 0$, $\delta\gamma < 0$, $(1 + \sigma\Delta T) \le \beta < \exp(\sigma\Delta T)$ and $\gamma > 0$.

*Proof.* Let $h : x \mapsto (1 + x)^{\frac{k}{x}}$ and $g : x \mapsto (\frac{1}{x(1+x)} - \frac{1}{x^2}\ln(1+x))$, for $x \in (0, k)$, with $k := \sigma\Delta T$. We have $h'(x) = kg(x)h(x)$ and $g$ is negative, then $h$ is decreasing on $(0, k)$. A direct calculation gives

$$\delta\beta = h(\sigma\Delta t) - h(\sigma\delta t) < 0.$$

Form $(1 + k) \le h(x) < \exp(k)$, we get $1 + \sigma\Delta T \le \beta < \exp(\sigma\Delta T)$.

Next, $\gamma > 0$ follows from (11) and $\beta - 1 \ge k > 1$. Moreover, we consider the function

$$p : x \mapsto \frac{h^2(x) - 1}{\sigma(2 + x)}, \quad x \in (0, k),$$

whose derivative is given by

$$x \mapsto -\frac{h^2(x) - 1}{\sigma(2 + x)^2} + \frac{2kg(x)h^2(x)}{\sigma(2 + x)},$$

which is negative, meaning that $p$ decreases on $(0, k)$. Hence

$$\delta\gamma = p(\sigma\Delta t) - p(\sigma\delta t) < 0,$$

and the result follows.

**Lemma 3.2.** Let $\sigma > 0$. Then we have $1 < \tau < \beta$.

*Proof.* Let us consider $q : x \mapsto \frac{h(x)+1}{2+x}$ on $(0, k)$, $k = \sigma\Delta T$ whose derivative is given by

$$x \mapsto -\frac{h(x) + 1}{(2 + x)^2} + \frac{kg(x)h(x)}{2 + x},$$
which is negative, meaning that $q$ is strictly decreasing on $(0, k)$. Thus, $q(\sigma\Delta t) < q(\sigma\delta t)$, i.e.,

$$ \frac{\beta + 1}{2 + \sigma\Delta t} < \frac{\beta_{\delta t} + 1}{2 + \sigma\delta t}. $$

It follows that

$$ \frac{\beta^2 - 1}{(2 + \sigma\Delta t)(\beta - 1)} < \frac{\beta_{\delta t}^2 - 1}{(2 + \sigma\delta t)(\beta_{\delta t} - 1)}, \quad \text{then} \quad \frac{\gamma}{\beta - 1} < \frac{\gamma_{\delta t}}{\beta_{\delta t} - 1}. $$

This inequality leads to $(\gamma\beta_{\delta t} - \gamma) < (\gamma_{\delta t}\beta - \gamma_{\delta t})$. Adding $\gamma\beta$ to each sides of the latter inequality, we get

$$ (\gamma\beta_{\delta t} - \gamma) + \gamma\beta < (\gamma_{\delta t}\beta - \gamma_{\delta t}) + \gamma\beta, \quad \text{i.e.,} \quad \beta\delta\gamma - \gamma\delta\beta < \delta\gamma. $$

Thus, dividing each member by $\delta\gamma$, we obtain $\beta - \gamma\delta\beta/\delta\gamma > 1$, i.e., $\tau > 1$. From the Lemma 3.1, we have $\gamma \frac{\delta\beta}{\delta\gamma} > 0$ so that $\tau = \beta - \gamma \frac{\delta\beta}{\delta\gamma} < \beta$, which concludes the proof. $\square$

In order to locate the roots of $\mathcal{P}$, we need the following lemma.

**Lemma 3.3.** Let $Q_{L,\beta}(z) := (z^{2L-2}-1)(z-1)-1+\frac{1}{\beta}$. If $L \ge 3$ and $\tau \ge \frac{3}{2}$ then $Q_{L,\beta}$ has only one root in $(1, \frac{3}{2})$.

*Proof.* The derivative of $Q_{L,\beta}$ with respect to $z$ is positive in $(1, \infty)$ which means that $Q_{L,\beta}$ increases in $(1, \infty)$ then in $(1, \frac{3}{2})$. Moreover, $Q_{L,\beta}(\frac{3}{2}) = \frac{3^{2L-2}-1}{2^{2L-1}} - 1 + \frac{1}{\beta}$ and $L \mapsto Q_{L,\beta}(\frac{3}{2})$ increases with $Q_{3,\beta}(\frac{3}{2}) = \frac{3}{2} + \frac{1}{\beta} > 0$. Therefore, $Q_{L,\beta}(\frac{3}{2}) > 0$ for $L \ge 3$. Furthermore, $Q_{L,\beta}(1) < 0$ since $\beta > 1$, hence $Q_{L,\beta}$ has only one root in $(1, \frac{3}{2})$. $\square$

Note that the explicit definition of $\tau$ in (22) enables to check a priori whether the assumption $\tau \ge \frac{3}{2}$ holds. We can, moreover, exhibit some choices of discretization parameters for which this assumption holds a priori. An example of a sufficient condition to obtain this inequality is presented in Appendix A; see, in particular, Lemma A.1 which shows that for $\delta t$ small enough, $\Delta t$ close enough to $\delta t$ and $\sigma, T, L$ given such that $3 \le L \le \frac{\sigma T}{k_0}$, with $k_0 := 0.681$, we have $\tau \ge \frac{3}{2}$.

**Lemma 3.4.** Let $L \ge 3, \tau \ge \frac{3}{2}$ and $f$ be defined by

$$ f(z) := -(\tau - z) \sum_{\ell=0}^{L-1} z^{2\ell}. $$

For $z_\theta := \beta + (\beta - \tau_0) \exp(i\theta)$, $\theta \in [0, \pi]$, we get $|f(z_\theta)| \ge (\tau - \tau_0)L$.

*Proof.* We first observe that

$$ |f(z_\theta)| = |\tau - z_\theta| \left| \sum_{\ell=0}^{L-1} z_\theta^{2\ell} \right|. $$

Let us now estimate each term $|\tau - z_\theta|$ and $\left|\sum_{\ell=0}^{L-1} z_\theta^{2\ell}\right|$. By definition we have

$$ |\tau - z_\theta| = |\tau - \beta - (\beta - \tau_0)e^{i\theta}|, $$

so that

$$ |\tau - z_\theta| \ge ||\tau - \beta| - |\beta - \tau_0|| = \tau - \tau_0. \quad (25) $$

Next, the square of the modulus of $z_\theta$ for $\theta \in [0, \pi]$ is given by $|z_\theta|^2 = (\beta + (\beta - \tau_0) \cos \theta)^2 + (\beta - \tau_0)^2 \sin^2 \theta$, which implies that $|z_\theta|^2 = \tau_0^2 + 4\beta(\beta - \tau_0) \cos^2 \frac{\theta}{2}$. Considering the change of variable $x = 4(\beta - \tau_0) \cos^2 \frac{\theta}{2}$ for $\theta \in [0, \pi]$, we get $|z_\theta|^2 = \tau_0^2 + \beta x$, with $x \in [0, 4(\beta - \tau_0)]$. We proceed analogously to get $|z_\theta - 1|^2 = (\tau_0 - 1)^2 + (\beta - 1)x$, and $|z_\theta + 1|^2 = (\tau_0 - 1)^2 + (\beta + 1)x$. In this way, let us consider the function $\psi$ defined by

$$ \psi(x) := (\tau_0 + \beta x)^L - L \left[ (\tau_0 - 1)^2 + (\beta - 1)x \right]^{\frac{1}{2}} \left[ (\tau_0 + 1)^2 + (\beta + 1)x \right]^{\frac{1}{2}} - 1, $$

for $x \in [0, 4(\beta - \tau_0)]$, whose derivative is given by

$$ \psi'(x) = L\varphi(x) + \frac{L}{2}(\beta + 1) \left( 1 - \frac{1}{\vartheta(x)} \right), $$
where $\varphi(x) = \beta(\tau_0^2 + \beta x)^{L-1} - \frac{1}{2}(\beta + 1) - \frac{1}{2}(\beta - 1)\vartheta(x)$ and

$$ \vartheta(x) = \left[ \frac{(\tau_0 + 1)^2 + (\beta + 1)x}{(\tau_0 - 1)^2 + (\beta - 1)x} \right]^{\frac{1}{2}}. $$

It is easily seen that $\vartheta(x) > 1$ on $[0, 4(\beta - \tau_0)]$. Then we get

$$ 1 - \frac{1}{\vartheta(x)} > 0. \qquad (26) $$

Furthermore, the functions $x \mapsto \beta(\tau_0^2 + \beta x)^{L-1} - \frac{1}{2}(\beta + 1)$ and $x \mapsto -(\beta - 1)\vartheta(x)$ are increasing in $[0, 4(\beta - \tau_0)]$. Then $\varphi$ increases and $\varphi(0) = \beta Q_{L,\beta}(\tau_0)/(\tau_0 - 1)$. Since $Q_{L,\beta}(\tau_0) = 0, \varphi(0) = 0$ we get

$$ \varphi(x) \ge 0, \text{ on } [0, 4(\beta - \tau)]. \qquad (27) $$

Combining (26) and (27), we get $\psi'(x) > 0$ on $[0, 4(\beta - \tau_0)]$. Hence $\psi$ is increasing and $\psi(0) = (\tau_0^{2L} - 1) - L(\tau_0^2 - 1)$. Recalling that $\tau_0 > 1$, we get

$$ \frac{\tau_0^{2L} - 1}{\tau_0^2 - 1} = \sum_{\ell=0}^{L-1} \tau_0^{2\ell} > L, $$

which implies that $\psi(0) > 0$. Thus $\psi(x) > 0$ on $[0, 4(\beta - \tau_0)]$, which is equivalent to

$$ (|z_{\theta}|^{2L} - 1) - L|z_{\theta} - 1| |z_{\theta} + 1| > 0, \theta \in [0, \pi]. $$

As a consequence,

$$ \left| \sum_{\ell=0}^{L-1} z_{\theta}^{2\ell} \right| \ge \frac{|z_{\theta}|^{2L} - 1}{|z_{\theta}^2 - 1|} > L. \qquad (28) $$

Combining (25) and (28), we get $|f(z_{\theta})| > (\tau - \tau_0)L$, for $\theta \in [0, \pi]$, and the proof is complete. $\square$

We now need to estimate $|f(z)|$ in order to apply Rouché's theorem. Numerically, we observe $|f(z_{\theta})| \ge (\tau - 1)L$ on $\{z_{\theta} = \beta + (\beta - e^{i\theta}), \theta \in [0, \pi]\}$. However, we encountered difficulties in proving this result analytically. On the other hand, we have established in Lemma 3.4 that if $\tau \ge \frac{3}{2}$, $|f(z_{\theta})| \ge (\tau - \tau_0)L$ on contracted half-circle $\{z_{\theta} = \beta + (\beta - \tau_0)e^{i\theta}, \theta \in [0, \pi]\}$ (the existence of $\tau_0$ is guaranteed in this case, see Lemma 3.3). This inequality gives rise to the next result.

**Theorem 3.1.** Let $\tau \ge \frac{3}{2}$ and $L_0$ be defined by

$$ L_0 := \frac{\beta - \tau}{\gamma(\tau - \tau_0)}. \qquad (29) $$

Then if $L > \max\{3, \alpha L_0\}$ the spectrum of $I - M_{\Delta t}^{-1}M_{\delta t}$ is contained in $D_{\sigma} \cup \{\mu^*\}$, where $D_{\sigma} = \{\mu \in \mathbb{C} : |\mu| < |\beta|/(\beta - \tau_0)\}$ and $\mu^*$ is a real negative number.

Assuming that $\alpha L_0 \ge 3$, this theorem requires the constraints $\tau \ge \frac{3}{2}$ and $L > \alpha L_0$. The former constraint has been discussed above. Let us assume that this constraint holds and focus on the latter constraint. From Lemma 3.2, we get $\beta > \tau$. Since $\tau \ge \frac{3}{2}$, $\beta > \frac{3}{2}$. Using Lemma A.2, we obtain in the case $\Delta t \to \delta t$ and $\delta t \to 0$

$$ \frac{L_0}{\sigma} = \frac{\beta - \tau}{\sigma\gamma(\tau - \tau_0)} = \frac{\delta\beta}{\sigma\delta\gamma(\tau - \tau_0)} < q_1(k) := \frac{8ke^k}{(2k+1)e^{2k}-1} \le \kappa := 1.31. $$

Therefore, $\alpha L_0 < \alpha\sigma q(k) \le \alpha\sigma\kappa$. It follows that $L > \alpha L_0$ if $L > \alpha\sigma\kappa$.

We conclude that, in the case $\Delta t \to \delta t$ and $\delta t \to 0$, the two conditions $\tau \ge \frac{3}{2}$ and $L > \alpha L_0$ hold simultaneously when $L$ satisfies

$$ \alpha\kappa \le \frac{L}{\sigma} < \frac{T}{k_0}. $$

*Proof.* We first locate the roots of $\mathcal{P}$. In this view, we consider the contour

$$ \mathcal{C} = \{z_{\theta} = \beta + (\beta - \tau_0)e^{i\theta} : \theta \in [-\pi, \pi]\}. \qquad (30) $$

Using (22), we observe that

$$ |\mathcal{P}(z_{\theta}) - f(z_{\theta})| = \alpha \frac{\beta - \tau}{\gamma}, $$
holds for $z_θ \in C$. Substituting (29) into $L > \alpha L_0$ gives rise to $(\tau - \tau_0)L > \alpha(\beta - \tau)/\gamma$. It follows that

$$|\mathcal{P}(z_θ) - f(z_θ)| < (\tau - \tau_0)L.$$

From Lemma 3.4, $|f(z_θ)| > (\tau - \tau_0)L$ for $z_θ = \beta + (\beta - \tau_0)e^{i\theta}$, $\theta \in [0, \pi]$. Moreover, for $\bar{z}_θ := z(-\theta) = \beta + (\beta - \tau_0)e^{-i\theta}$, $\theta \in [0, \pi]$, (28) holds, and

$$|\tau - \bar{z}_θ| = \tau - \beta - (\beta - \tau_0)e^{-i\theta}| \geq |\tau - \beta| - |\beta - \tau_0| = \tau - \tau_0.$$

Hence, $|f(z_θ)| > (\tau - \tau_0)L$ holds for $\theta \in [-\pi, 0]$. In follows that, $|f(z_θ)| > (\tau - \tau_0)L$ holds on $C$, so that

$$|\mathcal{P}(z_θ) - f(z_θ)| < |f(z_θ)|.$$

By Rouché's theorem (see [3]) we conclude that $\mathcal{P}$ and $f$ have the same number of roots in the domain delimited by $C$. The roots of $f$ are $z_0 = \tau$ and $z_j = e^{i\frac{j\pi}{L}}$ with $j = 1, \dots, L-1, L+1, \dots, 2L-1$. The roots $z_j \neq \tau$ are on the unit circle and only $z_0 = \tau \in C$, therefore $\mathcal{P}$ has only one root $a^* \in C$. The polynomial $\mathcal{P}$ has an even number of complex roots, since its coefficients are real numbers. As a consequence, $a^*$ must be real. On the other hand, given $x > \tau$, we have

$$\mathcal{P}(x) = \alpha \frac{\delta\beta}{\delta\gamma} + (x - \tau) \sum_{\ell=0}^{L-1} x^{2\ell} > 0,$$

meaning that $a^* < \tau$. Using (24), we get

$$\mu^* := \delta\beta/(\beta - a^*) < 0. \qquad (31)$$

Furthermore, the eigenvalues associated with the roots outside $C$ satisfy $|\mu| < |\beta|/(\beta - \tau_0)$, which means that they belong to $D_\sigma$, and the result follows. $\square$

**Lemma 3.5.** Let $\sigma > 0$, $T$ be given and $L$, $\Delta t$ and $\delta t$ be fixed. Then

$$\frac{|\delta\gamma|}{\gamma} \leq \sigma(\Delta t - \delta t) \left[ \frac{1}{2} + \left( \frac{\sigma}{2}(\Delta t - \delta t) + 1 \right) e^{2\sigma\Delta T} \right].$$

*Proof.* We have

$$\begin{aligned} \frac{|\delta\gamma|}{\gamma} &= \frac{(2 + \sigma\delta t)(\beta_{\delta t}^2 - \beta^2) + \sigma(\Delta t - \delta t)(\beta_{\delta t}^2 - 1)}{(2 + \sigma\delta t)(\beta^2 - 1)} \\ &= \frac{\beta_{\delta t}^2 - \beta^2}{\beta^2 - 1} + \frac{\sigma(\Delta t - \delta t)}{2 + \sigma\delta t} \left[ \frac{\beta_{\delta t}^2 - \beta^2}{\beta^2 - 1} + 1 \right]. \end{aligned}$$

It follows that

$$\frac{|\delta\gamma|}{\gamma} \leq \frac{\sigma}{2}(\Delta t - \delta t) + \left[ \frac{\sigma}{2}(\Delta t - \delta t) + 1 \right] \frac{\beta_{\delta t}^2 - \beta^2}{\beta^2 - 1}.$$

Let us focus on the term $(\beta_{\delta t}^2 - \beta^2)/(\beta^2 - 1)$. Define $h_0(x) = (1 + x)^{\frac{2k}{x}}$, $x \in (0, k)$, so that $\beta_{\delta t}^2 = h_0(\sigma\delta t)$ and $\beta^2 = h_0(\sigma\Delta t)$. The mean value theorem leads to

$$\beta_{\delta t}^2 - \beta^2 = (\Delta t - \delta t)|h_0'(x^*)|,$$

for some $x^* \in (\delta t, \Delta t)$, with $h_0'(x) = 2kg(x)h_0(x)$. Using the fact that $\ln(1 + x) \leq x$ for $x \geq 0$, we deduce that

$$|h_0'(x)| \leq \frac{2\sigma k}{1 + x} h_0(x) \leq 2\sigma k e^{2k},$$

which gives $\beta_{\delta t}^2 - \beta^2 \leq 2\sigma^2\Delta T e^{2\sigma\Delta T}$. Using Lemma 3.1, we obtain $\beta^2 - 1 \geq \sigma\Delta T(2 + \sigma\Delta T)$. Therefore

$$\frac{\beta_{\delta t}^2 - \beta^2}{\beta^2 - 1} \leq (\Delta t - \delta t) \frac{2\sigma e^{2\sigma\Delta T}}{2 + \sigma\Delta T} \leq \sigma(\Delta t - \delta t)e^{2\sigma\Delta T},$$

which is the desired result.

**Theorem 3.2.** Let $\tau \geq \frac{3}{2}$ and $L > \alpha L_0$. Then the spectral radius $\rho$ of the iteration matrix $I - M_{\Delta t}^{-1}M_{\delta t}$ satisfies

$$\rho < \sigma(\Delta t - \delta t) \left[ \frac{1}{2} + \left( \frac{\sigma}{2}(\Delta t - \delta t) + 1 \right) e^{2\sigma\Delta T} \right]. \qquad (32)$$
Proof. Theorem 3.1 implies that the spectral radius satisfies $\rho = |\mu^*|$, where $\mu^*$ is defined in (31).

Moreover, $\beta - a^* > \beta - \tau$ then, $1/(\beta - a^*) < 1/(\beta - \tau)$. Recalling (24), it follows that

$$|\mu^*| = \frac{|\delta\beta|}{\beta - a^*} < \frac{|\delta\beta|}{\beta - \tau} = \frac{|\delta\gamma|}{\gamma}. \qquad (33)$$

Inequality (32) then follows from Lemma 3.5. □

We show that the spectral radius of the iteration matrix in (13) can be bounded as in (32) under the conditions of Theorem 3.2. Thus, if the upper bound in (32) is less than 1, those conditions are sufficient to ensure convergence of the ParaOpt. This estimate is made under conditions $\tau \ge \frac{3}{2}$ and $L > \alpha L_0$. Note that for $L \le \alpha L_0$ we observe numerically that for some cases, the spectral radius remains bounded by 1. However, the proof of this result remains an open problem, which will be investigated in future work.

# 4 Numerical Results

In this section, we present some numerical experiments to test the accuracy of the bounds (32) and (33). We complete these experiments with a non-linear control example.

## 4.1 Linear scalar ODE: Sensitivity with respect to the discretization parameters

In this part, we set $\alpha = 0.5$, $T = 2$, $L = 10$ and investigate the dependence of the spectral radius of $I - \mathcal{M}_{\Delta t}^{-1}\mathcal{M}_{\delta t}$ on the discretization parameters $\delta t$ and $\Delta t$. In our test, $L$ satisfies $L > \alpha L_0$.

In our first test, we compute the true spectral radius of $I - \mathcal{M}_{\Delta t}^{-1}\mathcal{M}_{\delta t}$ and the estimates (32) and (33) in three different cases. In the first case, we fix $\sigma = 1$ and $\Delta t = \Delta T/4$ for various $\delta t$, that is, $\delta t = \frac{\Delta t}{2^k}$, $k = 1, \dots, 20$. Next, we fix $\sigma = 3.408$, $\delta t = 5.10^{-2} \cdot 2^{-20}$ and let $\Delta t$ vary according to $\Delta t = \frac{\Delta T}{2^k}$, $k = 1, \dots, 19$. Finally, we vary $\Delta t = \frac{\Delta T}{2^k}$, $k = 1, \dots, 19$ and fix $\frac{\Delta t}{\delta t} = 2$. The results are shown in Figure 2. In the tree cases, we observe that the

<figure><img src="image_3.png" /><figcaption>Figure 2: Spectral radius and its estimates. Top left: varying $\delta t$ (with fixed $\Delta t$); top right: varying $\Delta t$ (fixed $\delta t$); bottom: varying $\Delta t$ and $\delta t$ (with fixed ratio $\frac{\Delta t}{\delta t}$).</figcaption></figure>
upper bound in (33) and the exact spectral radius are very close. Note that the highest admissible value of $\delta t$ is $\Delta t$. In

<figure><img src="image_4.png" /><figcaption>Figure 3: Spectral radius and the estimate (32) for various values of $\alpha$.</figcaption></figure>

the case where $\Delta t \to \delta t$, both the spectral radius and the upper bounds tend to 0, which is consistent with the fact that $I - \mathcal{M}_{\Delta t}^{-1}\mathcal{M}_{\delta t}$ tends to the null matrix.

Moreover, we observe that the spectral radius and the upper bounds (33) and (32) increase as $\Delta t$ increases (see Figures (2) top right and bottom). On the other hand, when $\Delta t$ is fixed and $\delta t \ll \Delta t$, the spectral radius and the upper bounds remain constant (see Figure 2 top left). This behavior shows that the approximation of the Jacobian of $\mathcal{F}_{\delta t}(X)$ becomes coarser when $\Delta t$ increases, which affects the convergence rate (see [17]).

We also performed several tests where the constraint $L > \alpha L_0$ in Theorem 3.2 is not satisfied, that is $L \le \alpha L_0$. For example, keeping $T = 2$ and $L = 10$ and setting $\Delta t = \Delta T/50$, $\delta t = \Delta T/200$, $\sigma = 2$, we have $\tau < \frac{3}{2}$. Moreover, for $\alpha$ from $10^1$ to $10^4$, $L \le \alpha L_0$, however, (32) holds; see Figure 3. Consequently, further investigations are required to refine the analysis and extend the application range of (32).

## 4.2 A nonlinear example

We finally consider a nonlinear example which consists of minimizing the cost functional

$$ \mathcal{J}(\nu) = \frac{1}{2}|y(T) - y_{\text{target}}|^2 + \frac{\alpha}{2} \int_0^T |\nu(t)|^2 dt, \quad (34) $$

subject to the state dynamic of an inverted pendulum described by

$$ \begin{aligned} \dot{y}_1(t) &= y_2(t) \\ \dot{y}_2(t) &= \frac{g}{l} \sin(y_1(t)) + \nu(t) \cos(y_1(t)), \end{aligned} \quad (35) $$

with $t \in [0, T]$, $y(0) = y_{\text{init}}$, $g$ the gravity acceleration, $l$ the length of the pendulum and $\nu$ the horizontal acceleration (see [16]). The first component $y_1$ is the angular position of the pendulum in $t$ and its angular velocity $y_2$. Contrary to the case of the linear control, the control intervenes here in a multiplicative way, i.e., in the factor of the state. This means that the ParaOpt approach cannot be applied directly, and we need to use implicit function-type techniques to eliminate the control and get the reduced optimality system. In this way, we rewrite (35) as $\dot{y} = f(y, \nu)$ where

$$ f(y, \nu) := \left( \begin{array}{c} y_2(t) \\ \int_0^t \sin(y_1(t)) + \nu(t) \cos(y_1(t)) \end{array} \right). $$

The first optimality system reads to

$$ \dot{y} = f(y, \nu), \quad \dot{\lambda} = -\nabla_y f^T(y, \nu) \cdot \lambda, \quad (36) $$

and the algebraic equation

$$ \alpha\nu = -\nabla_{\nu} f^T(y, \nu) \cdot \lambda, \quad (37) $$

with the initial and final conditions $y(0) = y_{\text{init}}$ and $\lambda(T) = y(T) - y_{\text{target}}$ respectively. From (37), we get $\nu(t) = \phi(y(t), \lambda(t)) := \frac{1}{\alpha}\lambda_2(t) \cos y_1(t)$. Substituting $\phi$ by the latter expression into (36), we obtain the reduced optimality system

$$ \begin{aligned} \dot{y} &= f(y, \phi(y, \lambda)) \\ \dot{\lambda} &= -\nabla_y f^T(y, \phi(y, \lambda)) \cdot \lambda \end{aligned} \quad (38) $$
<figure><img src="image_5.png" /><figcaption>Figure 4: $L^2$-norm of the convergence errors for various values of the ratio $r = \frac{\delta t}{\Delta t}$.</figcaption></figure>

<table><caption>Table 1: Equivalent sequential time of ParaOpt in the case of the forward Euler method solver for inverted pendulum.</caption><thead><tr><td>L</td><td># Iters</td><td>L.T<sub>CPU</sub>(L)</td></tr></thead><tbody><tr><td>8</td><td>8</td><td>80.258993</td></tr><tr><td>16</td><td>8</td><td>72.285331</td></tr><tr><td>32</td><td>9</td><td>73.279832</td></tr><tr><td>64</td><td>9</td><td>74.521707</td></tr></tbody></table>

where

$$f(y, \phi(y, \lambda)) = \left( \frac{g}{l} \sin y_1 - \frac{1}{\alpha} \lambda_2 \cos^2 y_1 \right) \quad \text{and} \quad \nabla_y f^T(y, \phi(y, \lambda)).\lambda := \left( \frac{g}{l} \lambda_2 \cos y_1 + \frac{1}{2\alpha} \lambda_2^2 \sin 2y_1 \right),$$

with initial and final conditions $y(0) = y_{init}$ and $\lambda(T) = y(T) - y_{target}$ respectively. We now numerically solve the optimal control problem of our example applying the ParaOpt to (38). In this test, the time discretization of the ParaOpt method is performed with the forward Euler method. The error is defined as the $L^2$-norm of the difference between the state and adjoint state values obtained from a convergence fine-grid solution and the state and adjoint state values at each iteration of ParaOpt. We set $T = 1$, $\alpha = 0.01$, $L = 8$, $\delta t = 5e-6$, $g = 9.81m/s^2$, $l = 0.5m$, $y_{init} = (\frac{\pi}{4}, \frac{\pi}{6})$, and $y_{target} = (0, 0)$. The results are shown in Figure 4 where the error is represented as a function of the ratio $r = \frac{\delta t}{\Delta t}$.

For $r = 1$, the ParaOpt actually corresponds to the true Newton method, which explains the observed quadratic convergence. On the other hand, for values of $r$ smaller than 1, $\Delta t$ increases and the approximation of the Jacobian becomes a coarse approximation of the true Jacobian. This leads to a decrease in the rate of convergence (see [17]).

We observe that for $L = 1$ the algorithm does not converge, which means that without time domain decomposition, the iterative solver diverges even when working with the true Jacobian. Moreover, for $L > 1$, if $T/L \le 1/5$ the algorithm converges and diverges otherwise.

In this test, we consider the sequential time $L.T_{CPU}(L)$ where $T_{CPU}(L)$ is the required time to convergence to a tolerance of $10^{-13}$ when $L$ processors are used for a parallel implementation of the algorithm. Then we fix $\delta t = 2.5e-6$ to extend the choice of the number of processors $L$, $\delta t/\Delta t = 10^{-1}$ and we keep the other parameters as before. Since the algorithm does not converge for $L < 5$, the corresponded sequential times $L.T_{CPU}(L)$ are not defined. We then present the results where $L \ge 5$ in Table 1. We observe that the ParaOpt has the weak scalability property in the sense that the computing time remains roughly constant as the number of processors increases. Note that the weak scalability of ParaOpt has also been established in a case of linear control (i.e. the control enters as a source term in the state equation) of nonlinear dynamics, namely Lotka-volterra system in [5].

## References

[1] Andrew T. Barker and Martin Stoll. Domain decomposition in time for pde-constrained optimization. *Comput. Phys. Commun.*, 197:136–143, 2015.

[2] M. Emmett and M. L. Minion. Toward an Efficient Parallel in Time Method for Partial Differential Equations. *Communications in Applied Mathematics and Computational Science*, 7:105–132, 2012.

[3] T. W. Gamelin. *Complex Analysis*. Springer-Verlag New York, Inc, 2001.
[4] M. J. Gander. 50 years of time parallel time integration, in multiple shooting and time domain decomposition methods. Springer, New York, pages 69–113, 2015.

[5] M. J. Gander, F. Kwok, and J. Salomon. ParaOpt: A parareal algorithm for optimality systems. SIAM Journal on Scientific Computing, 42(5):A2773–A2802, 2020.

[6] M. J. Gander, F. Kwok, and G. Wanner. Constrained optimization: from lagrangian mechanics to optimal control and pde constraints. in Optimization with PDE Constraints, Springer, pages 151–202, 2014.

[7] Martin J Gander. Optimized schwarz methods. SIAM Journal on Numerical Analysis, 44(2):699–731, 2006.

[8] Martin J Gander and Stefan Güttel. Paraexp: A parallel integrator for linear initial-value problems. SIAM Journal on Scientific Computing, 35(2):C123–C142, 2013.

[9] Martin J Gander and Felix Kwok. Schwarz methods for the time-parallel solution of parabolic control problems. In Domain decomposition methods in science and engineering XXII, pages 207–216. Springer, 2016.

[10] Martin J Gander and Stefan Vandewalle. Analysis of the parareal time-parallel time-integration method. SIAM Journal on Scientific Computing, 29(2):556–578, 2007.

[11] S. Götschel and M.L. Minion. Parallel-in-time for parabolic optimal control problems using pfasst, in domain decomposition methods in science and engineering xxiv. Springer, pages 363–371, 2019.

[12] SEBASTIAN Gotschel and Michael L Minion. An efficient parallel-in-time method for optimization with parabolic pdes. SIAM Journal on Scientific Computing, 41(6):C603–C626, 2019.

[13] Stefanie Günther, Nicolas R Gauger, and Jacob B Schroder. A non-intrusive parallel-in-time approach for simultaneous optimization with unsteady pdes. Optimization Methods and Software, 34(6):1306–1321, 2019.

[14] E. Hairer and G. Wanner. Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems. Springer-Verlag Berlin Heidelberg, 1996.

[15] K. Ito and K. Kunisch. Lagrange Multiplier Approach to Variational Problems and Applications. Society for Industrial and Applied Mathematics, 2008.

[16] H. P. Kalmus. The inverted pendulum. American journal of Physics, 38(7):874–778, 1970.

[17] C. T. Kelley. Solving Nonlinear Equations with Newton’s Method. Society for Industrial and Applied Mathematics, 2003.

[18] Felix Kwok and Djahou N Tognon. A Parallel in Time Algorithm Based on ParaExp for Optimal Control Problems. In 2024 IEEE 63rd Conference on Decision and Control (CDC), MiCo, Milan, Italy, December 2024.

[19] J.-L. Lions and G. Maday, Y.and Turinici. Résolution d’edp par un schéma en temps ’pararéel’. Comptes Rendus de l’Académie des Sciences-Series I-Mathematics, 332(7):661–668, 2001.

[20] Yvon Maday, Julien Salomon, and Gabriel Turinici. Monotonic parareal control for quantum systems. SIAM JOURNAL ON NUMERICAL ANALYSIS, 45(6):2468–2482, 2007.

[21] Tarek P Mathew, Marcus Sarkis, and Christian E Schaerer. Analysis of block parareal preconditioners for parabolic optimal control problems. SIAM Journal on Scientific Computing, 32(3):1180–1200, 2010.

[22] David D Morrison, James D Riley, and John F Zancanaro. Multiple shooting method for two-point boundary value problems. Communications of the ACM, 5(12):613–614, 1962.

[23] M. K. Riahi, J. Salomon, S. J. Glaser, and D. Sugny. Fully efficient time-parallelized quantum optimal control algorithm. PHYSICAL REVIEW A, 93(4), APR 13 2016.

[24] Emmanuel Trélat. Contrôle optimal: théorie & applications. Vuibert, Paris, 2nd edition, 2008.

A Additional Results

**Lemma A.1.** Let $h(x) = (1 + x)^{\frac{k}{x}}$, $p(x) = \frac{h^2(x)-1}{2+x}$ and for $0 < x < y$,

$$ \chi(x, y) = h(y) - p(y) \left( \frac{h(y) - h(x)}{p(y) - p(x)} \right). $$

If $k \ge k_0 = 0.681$ then

$$ \lim_{x \to 0} \lim_{y \to x} \chi(x, y) > \frac{3}{2}. $$
Proof. We first have the following equality

$$ \lim_{y \to x} \frac{h(y) - h(x)}{p(y) - p(x)} = \lim_{y \to x} \frac{\frac{h(y)-h(x)}{y-x}}{\frac{p(y)-p(x)}{y-x}} = \frac{h'(x)}{p'(x)}. \quad (39) $$

Then,

$$ \lim_{x \to 0} \lim_{y \to x} \chi(x, y) = \lim_{x \to 0} \left( h(x) - p(x) \frac{h'(x)}{p'(x)} \right) = e^k - \frac{(e^{2k} - 1)}{2} \left( \lim_{x \to 0} \frac{h'(x)}{p'(x)} \right). $$

Next, $h'(x) = g(x)h(x)$ where $g(x) = \frac{k}{x(1+x)} - \frac{k}{x^2} \ln(1+x)$. Since

$$ \ln(1+x) = x - \frac{x^2}{2} + x \left( \frac{x^2}{2} - \frac{x^3}{4} + \dots \right), $$

the Taylor expansion of $g$ is given by

$$ \begin{aligned} g(x) &= k \left[ \frac{1}{x(1+x)} - \frac{1}{x} + \frac{1}{2} - x \left( \frac{1}{3} - \frac{x^2}{4} + \dots \right) \right], \\ &= -\frac{k}{1+x} + \frac{k}{2} - kx \left( \frac{1}{3} - \frac{x^2}{4} + \dots \right), \end{aligned} $$

so that, $\lim_{x \to 0} g(x) = -k/2$ and $\lim_{x \to 0} h'(x) = -\frac{k}{2}e^k$. Besides, $p'(x) = -\frac{h^2(x)-1}{(2+x)^2} + \frac{g(x)h(x)}{2+x}$ and thus

$$ \lim_{x \to 0} p'(x) = -\frac{e^{2k} - 1}{4} - \frac{ke^{2k}}{2}. $$

Therefore,

$$ \lim_{x \to 0} \frac{h'(x)}{p'(x)} = \frac{2ke^k}{(2k+1)e^{2k}-1}, \quad (40) $$

and

$$ \lim_{x \to 0} \lim_{y \to x} \chi(x, y) = \frac{(k+1)e^{3k} + (k-1)e^k}{(2k+1)e^{2k}-1}. $$

Now, let us consider

$$ \frac{(k+1)e^{3k} + (k-1)e^k}{(2k+1)e^{2k}-1} - \frac{3}{2} = \frac{\chi_0(k)}{2((2k+1)e^{2k}-1)}, $$

where $\chi_0(k) := 2(k+1)e^{3k} - 3(2k+1)e^{2k} + 2(k-1)e^k + 3$. The function $\chi_0$ is differentiable and $\chi'_0(k) = \phi(k)e^k$, $\phi(k) = (6k+8)e^{2k} - 12(k+1)e^k + 2k$. Also, $\phi'(k) = (12k+14)e^{2k} - 12(k+2)e^k + 2$ and $\phi^{(2)}(k) = 4(6k+14)e^{2k} - 12(k+3)e^k$.

Since $\phi^{(2)}(k) > 0$ for $k > 0$, $\phi'$ increases in $(0, \infty)$ with $\phi'(0.38) < 0 < \phi'(0.39)$. Thus, $\phi'(k) = 0$ has only one solution $\kappa_0$ that satisfies $0.38 < \kappa_0 < 0.39$. It follows that $\phi'(k) < 0$ in $(0, \kappa_0)$ and $\phi'(k) > 0$ in $(\kappa_0, \infty)$. As a consequence, $\phi$ decreases in $(0, \kappa_0)$ and increases in $(\kappa_0, \infty)$. As $\phi(0) < 0$, $\phi(k) > 0$ in $(0, \kappa_0)$. In addition, as $\phi(0.45) \cdot \phi(0.46) < 0$, the sign of $\phi$ changes in $(\kappa_0, \infty)$ and therefore $\phi(k) = 0$ has only one solution $\kappa_1$ in $(0.45, 0.46)$. Then, $\phi(k) < 0$ in $(0, \kappa_1)$ and $\phi(k) > 0$ in $(\kappa_1, \infty)$, which implies that $\chi'_0(k) < 0$ in $(0, \kappa_1)$ and $\chi'_0(k) > 0$ in $(\kappa_1, \infty)$. It follows that $\chi_0$ decreases in $(0, \kappa_1)$ and increases in $(\kappa_1, \infty)$. Furthermore, $\chi_0(0.680) \cdot \chi_0(0.681) < 0$, therefore, $\chi_0(k) = 0$ has only one solution $\kappa_2$ in $(0.68, 0.681) \subset (\alpha_1, \infty)$. As a consequence, $\chi_0(k) > 0$ for $k \ge \kappa_0 = 0.681 > \kappa_2$, so that,

$$ \frac{\chi_0(k)}{2(2k-1)e^{2k}-2} > 0, $$

which achieves the proof.

<figure><img src="image_6.png" /></figure>

**Lemma A.2.** Let $P_{\ell, \kappa}(z) = (z^{2\ell-2}-1)(z-1)-1+\frac{1}{\kappa}$ where $\kappa > \frac{3}{2}$ and $\ell \ge 3$. We have

$$ \lim_{x \to 0} \lim_{y \to x} \frac{h(y) - h(x)}{(p(y) - p(x))(\chi(x, y) - z_0)} < \kappa_3, $$

where $\kappa_3 = 1.31, k \ge k_0$ and $z_0$ the real root of $P_{\ell, \kappa}$ in $(1, \frac{3}{2})$.
Proof. We first estimate $\lim_{x \to 0} \lim_{y \to x} \chi(x, y) - z_0$. Since $\kappa > \frac{3}{2}$, $\frac{1}{\kappa} < \frac{2}{3}$ and then,

$$P_{\ell}(z) := (x^{2\ell-2}-1)(x-1)-1+\frac{2}{3} = (x^{2\ell-2}-1)(x-1)-\frac{1}{3} > P_{\ell,\kappa}(z).$$

The polynomial $P_\ell$ increases in $(1, 2)$. Moreover, $P_\ell(\frac{5}{4}) = \frac{5^{2\ell-2}-1}{4^{2\ell-1}} - \frac{1}{3}$ increases as function of $\ell$ with $P_3(\frac{5}{4}) = \frac{5^4-1}{4^5}-\frac{1}{3} > 0$. Thus, $P_\ell(\frac{5}{4}) > 0$ for $\ell \ge 3$. As $P_\ell(1) < 0$, we conclude that $P_\ell$ has a real root $z_1$ in $(1, \frac{5}{4})$. Now, because $P_\ell(z) > P_{\ell,\kappa}(z)$ we get $z_1 > z_0$. Using Lemma A.1 we obtain $\lim_{x \to 0} \lim_{y \to x} \chi(x, y) - z_0 > \frac{3}{2} - z_0 > \frac{3}{2} - z_1$. It follows that

$$\lim_{x \to 0} \lim_{y \to x} \chi(x, y) - z_0 \ge \frac{3}{2} - \frac{5}{4} = \frac{1}{4},$$

since $z_1 < \frac{5}{4}$. Hence,

$$\lim_{x \to 0} \lim_{y \to x} \frac{1}{\chi(x, y) - z_0} < 4. \qquad (41)$$

Next, using (39) and (40), we get

$$\lim_{x \to 0} \lim_{y \to x} \frac{h(y) - h(x)}{p(y) - p(x)} = \frac{2ke^k}{(2k+1)e^{2k}-1},$$

and together with (41), we obtain

$$\lim_{x \to 0} \lim_{y \to x} \frac{h(y) - h(x)}{(p(y) - p(x))(\chi(x, y) - z_0)} < q_1(k), \quad q_1(k) := \frac{8ke^k}{(2k+1)e^{2k}-1}.$$

Furthermore,

$$q_1'(k) = 8 \frac{-(2k^2 + k - 1)e^{3k} - (k + 1)e^k}{((2k + 1)e^{2k} - 1)^2} = - \frac{8(k + 1)e^k ((2k - 1)e^{2k} + 1)}{((2k + 1)e^{2k} - 1)^2} < 0 \text{ on } k \ge 0.$$

Then $q_1$ decreases in $(0, \infty)$, so that for $k > 0$, $q_1(k) < 2$ and for $k \ge k_0$, $q_1(k) \le q_1(k_0) = \kappa_3 = 1.31$, which achieves the proof. $\square$