# Decompose-then-optimize: a new approach to design domain decomposition methods for optimal control problems

Pierre-Henri Cocquet, Martin Gander, Alexandre Vieira

► To cite this version:

Pierre-Henri Cocquet, Martin Gander, Alexandre Vieira. Decompose-then-optimize: a new approach to design domain decomposition methods for optimal control problems. Numerical Algorithms, 2025, 10.1007/s11075-025-02046-4. http://dx.doi.org/10.1007/s11075-025-02046-4

HAL Id: hal-04999634

https://hal.science/hal-04999634v1

Submitted on 20 Mar 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
Decompose-then-optimize: a new approach to
design domain decomposition methods for
optimal control problems.

Pierre-Henri Cocquet<sup>1a,1b*</sup>, Martin J. Gander<sup>2*</sup>
and Alexandre Vieira<sup>3*</sup>

<sup>1a*</sup>Université de Pau et des Pays de l'Adour, E2S-UPPA, SIAME, Pau,
64000, France.

<sup>1b*</sup> IRL CRM-CNRS - International Research Lab.

<sup>2*</sup>UNIGE, Genève, Switzerland.

<sup>3*</sup>Université Côte d'Azur, INRIA, LJAD, Nice, 06100, France.

*Corresponding author(s). E-mail(s): pierre-henri.cocquet@univ-pau.fr; martin.gander@unige.ch; alexandre.vieira@inria.fr;

**Abstract**

For optimal control problems there is a classical discussion of whether one should first optimize the problem and then discretize it, or the other way round. We are interested in exploring a similar question related to domain decomposition methods for optimal control problems which have received substantial attention over the past two decades, but new methods were mostly developed using the optimize-then-decompose approach. After a detailed introduction to this subject, we present and analyze a new domain decomposition method for optimal control problems that comes from the decompose-then-optimize strategy which is less common. We use as our model problem a linear quadratic optimal control problem which we decompose and then solve using an augmented Lagrangian optimization technique. This leads to a new domain decomposition algorithm for such problems that has very good scalability properties. We prove that, when the algorithm converges, it necessarily converges to an optimal point of the original, non-decomposed problem. We illustrate the efficiency of our new domain decomposition method with numerical examples from which we obtain very desirable properties for domain decomposition methods, namely that the convergence is independent of the meshsize and the number of subdomain.

**Keywords:** Domain decomposition, Optimal control, Augmented Lagrangian.
<figure><img src="image_2.png" /><figcaption>Fig. 1: Discretize then optimize versus optimize then discretize (left part) and decompose then optimize versus optimize then decompose (right part).</figcaption></figure>

# 1 Introduction

In optimal control and PDE constrained (equivalently constraint is commonly used) optimization problems of the form

$$ \min_{y,u} \mathcal{J}(y, u) \quad \text{s.t.} \quad g(y, u) = 0, \qquad (1) $$

a subject of discussion is the question whether one should first optimize and then discretize the problem, i.e. compute the Lagrangian first order optimality conditions

$$ \nabla \mathcal{L}(y, u, \lambda) = 0, \qquad (2) $$

and then discretize them (optimize-then-discretize or indirect method, path going from the middle down and then to the left in Figure 1), or if it is better to first discretize and then optimize the problem (discretize-then-optimize or direct method, path going left and then down in Figure 1), see [51] for a presentation of these methods. There are advantages to both approaches: for discretize-then-optimize, one always obtains the true gradient of the discrete problem, even when the discretization is coarse, and symmetric formulations remain symmetric. Also, for time dependent problems, due to the Pontryagin maximum principle (see [26] for a historical introduction), the first order optimality conditions are of the form of a Hamiltonian differential equation, and if one discretizes first and then optimizes, this geometric structure can automatically be preserved [14], see also [7, 29] and [30, Chapter VI, exercises 14,15,16]. This is rather elegant, and can be important in low dimensions close to a critical value of the Hamiltonian, but not in general since the first order optimality condition is a boundary value problem and the time interval is not very large in general. In optimize-then-discretize, none of the above properties hold, but discretizations are much more flexible, one can adapt locally independently in the forward and backward problem which can be an advantage, also in hyperbolic problems where CFL conditions need to be met. In some cases, the two approaches also commute.

We are interested here in a new, analogous question for optimal control problems and PDE constrained optimization when designing and analyzing parallel algorithms for their solution using domain decomposition (DD), namely the decompose-then-optimize and optimize-then-decompose approaches, see the right part of Figure 1. The
situation here is more involved, and much less explored in the literature. Historically, the FETI method [19] was developed in this spirit, albeit for unconstrained optimization, and using the equivalence of the Laplace problem and the minimization of the Dirichlet integral, see the blue arrows in Figure 1: for FETI, one first considers the PDE problem as an equivalent minimization problem (blue arrow going up), then decomposes the minimization problem at the variational level (red arrow going to the right and blue arrow going down) which makes Neumann traces match automatically, and thus only the Dirichlet traces need to be imposed to match explicitly for the DD solution to be a solution of the underlying PDE which is done in FETI using Lagrange multipliers. FETI is also historically not written as an iteration, like the equivalent dual Schur complement method at the discrete level, and the Conjugate Gradient method is used to solve the decomposed system. Using however a stationary residual correction method for the FETI system shows that it is a classical Dirichlet-Dirichlet method (blue arrow going down), the dual of the Neumann-Neumann method [9] or primal Schur complement method at the discrete level, see e.g. [15, Section 4.8].

The *optimize-then-decompose* approach (black arrow going down in the middle and then right green arrow in Figure 1) has received quite some attention in the literature, since one can directly apply standard DD methods to the optimality system (2), and then study their convergence, which leads to interesting new results for DD methods, see for example [21–23, 25], and also [45]. In contrast to the *optimize-then-discretize* approach, in the *optimize-then-decompose* approach, it is also often shown that the resulting DD iterations can be interpreted as solving optimal control problems on the subdomains during the DD iteration, i.e. the green arrow going up on the right in Figure 1, see e.g. [16, 17, 21, 25, 27, 43, 44, 53]. DD methods can also be interpreted as an optimization problem, see for example [18], where optimization techniques are used to minimize the jumps in interface traces. This approach is also related to the *indirect shooting method*, used for numerically solving optimal control problems of ODEs where the optimality conditions are first derived, and then a parallel-in-time method is applied to the derived system; see [31] for more details.

We focus here on the rather new *decompose-then-optimize* approach indicated by the red arrows going right and down in Figure 1. The idea here is to add some *continuity constraints* in the feasible set associated to the decomposed PDEs. These constraints must then be handled by different optimization techniques. This approach is related to the *direct shooting method* for solving optimal control problems of ODEs where the continuity of the solution in time also appears as constraints in the optimization problem [31]. Due to the many optimization techniques, one can discover new DD methods doing this, and even try to interpret their meaning as DD methods for the first order optimality system, as indicated by the left going red arrow in Figure 1. A very fruitful source for such new methods is the augmented Lagrangian technique introduced by Hestenes in [32], see also Powell [50]. In this approach one has, for (some of) the constraints, a penalization function composed of a quadratic penalty term and a scalar product involving a Lagrange multiplier. There is then a precise iterative algorithm on how to update the Lagrange multiplier and the penalty parameter in order to converge to an optimum, the augmented Lagrangian method. Without
DD, this method has been successfully applied to many PDE-constrained optimization problems over the past decades, see [1, 3–5, 11–13, 33, 36–38], and the method is suitable for the analysis in infinite dimensional spaces, see for instance [39–41].

Our goal is to show how one can design and analyze new DD algorithms based on the principle of decompose-and-optimize, and we follow Emile Picard's recommendation [49] to do so: "Les méthodes d'approximation dont nous faisons usage sont théoriquement susceptibles de s'appliquer à toute équation, mais elles ne deviennent vraiment intéressantes pour l'étude des propriétés des fonctions définies par les équations différentielles que si l'on ne reste pas dans les généralités et si l'on envisage certaines classes d'équation<fnref n="1" />." We will thus choose a specific PDE-constrained optimization problem and use the augmented Lagrangian method to solve it to discover a new domain decomposition method based on the decompose-then-optimize approach.

Our paper is organized as follows: in Section 2, we present the decompose-then-optimize method based on the augmented Lagrangian algorithm for a linear quadratic optimization problem. In Section 3 we study the convergence of the new DD method for 2 subdomains, and explain how the results can be extended to more general decompositions; we also use Fourier techniques to get more insight into the convergence behavior. Afterwards, we show numerically the excellent scalability properties of our algorithm in Section 4, and we draw some conclusions in Section 5.

*Notations.* We denote by $\nabla q$ the gradient of a real-valued function. Assuming that we have a Hilbert space $\mathcal{H}$ and a subspace $\mathcal{X}$ such that $\mathcal{X} \subset \mathcal{H} \subset \mathcal{X}^*$ is a Gelfand triple, the directional derivative of a function $F: x \in \mathcal{X} \mapsto F(x) \in \mathbb{R}$ is

$$ \partial_x F(x)[\delta x] = \lim_{t \to 0} \frac{F(x + t \delta x) - F(x)}{t} = \langle \partial_x F(x), \delta x \rangle_{\mathcal{X}^*, \mathcal{X}} $$

where $\partial_x F(x)$ is the gradient of $F$. The notation $A \lesssim B$ means that there exists a constant $C$ which can depend only on the domain $\Omega$, $C = C(\Omega)$, such that $A \le C(\Omega)B$. For $\Gamma \subset \partial\Omega$, we denote by $H^s(\Gamma)$ the restriction to $\Gamma$ of distributions in $H^s(\partial\Omega)$, and by $H_{00}^{1/2}(\Gamma)$ the set of distributions defined on $\Gamma$ such that their extension with 0 on $\partial\Omega$ is in $H^{1/2}(\partial\Omega)$. We refer to [46] for more information about these trace spaces.

## 2 A decompose-then-optimize method based on the augmented Lagrangian algorithm

As indicated above, we use as our model problem the linear quadratic PDE constrained optimization problem<fnref n="2" />

$$ \begin{array}{ll} \min & \tilde{\mathcal{J}}(y, u) := \frac{1}{2}\|y - y_{\text{target}}\|_{L^2(\Omega)}^2 + \frac{\alpha}{2}\|u\|_{L^2(\Omega)}^2 \\ \text{s.t.} & \begin{cases} -\Delta y = f + u & \text{on } \Omega, \\ y|_{\partial\Omega} = 0, \end{cases} \end{array} \quad (3) $$

<footnote n="1">¹The successive approximation methods we use can in theory be applied to any equation, but they only become really interesting for the study of properties of functions defined by differential equations if one does not remain in generalities, and considers certain classes of equations.</footnote>

<footnote n="2">²We use $\tilde{\mathcal{J}}$ here to have $\mathcal{J}$ for the decomposed optimization problem.</footnote>
where $\Omega$ is a bounded open subset of $\mathbb{R}^d$, $d \ge 1$, with Lipschitz boundary, $u \in L^2(\Omega)$ is our control and $f \in L^2(\Omega)$ is an imposed source term.

For FETI, an ad hoc augmented Lagrangian formulation penalizing the Dirichlet jumps has already been proposed in [47] which led to Robin transmission conditions when interpreted as a DD method, see also [2] for an augmented Lagrangian formulation with inequality constraints for FETI. For the PDE constrained optimization problem (3) however, we are only aware of [28, 35], where a quadratic penalization approach was studied with the PDE constraint decomposed into subdomains, and the trace jumps at the interfaces were then penalized in the augmented cost formulation. In [42], the authors also use an augmented Lagrangian method but to penalize the entire PDE, in addition to a decomposition technique, with no proof on convergence.

We decompose the domain of the PDE constraint in (3) into $J$ non-overlapping subdomains, $\bar{\Omega} = \cup_{j=1}^J \bar{\Omega}_j$, $\bar{\Omega}_i \cap \bar{\Omega}_j = \emptyset$, and then add continuity constraints along the interfaces $\Gamma_{ij} := \partial\Omega_i \cap \partial\Omega_j = \Gamma_{ji}$. We denote the normal on $\Gamma_{ij}$ outward from $\Omega_i$ into $\Omega_j$ by $n_{ij}$. The decomposed optimization problem then reads

$$ \begin{array}{ll} \min_{y,u,g} & \mathcal{J}(y, u, g) = \sum_{j=1}^{J} \frac{1}{2} \|y_j - y_{\text{target}}\|_{L^2(\Omega_j)}^2 + \frac{\alpha}{2} \|u_j\|_{L^2(\Omega_j)}^2 \\ \text{s.t.} & \left\{ \begin{array}{ll} -\Delta y_j = f + u_j & \text{in } \Omega_j, \\ y_j = 0 & \text{on } \partial\Omega \cap \Omega_j, \\ \partial_{n_{ji}} y_j |_{\Gamma_{ji}} = g_{ji}, \quad g_{ji} = -g_{ij} & \text{on } \Gamma_{ij} \neq \emptyset, \\ y_j = y_i & \text{on } \Gamma_{ij} \neq \emptyset, \\ (u_j, g_j) \in \mathcal{U}. \end{array} \right. \\ & \text{(4)} \end{array} $$

<figure><img src="image_3.png" /><figcaption>**Fig. 2:** Example decomposition of $\Omega$ into four subdomains $\Omega_j$.</figcaption></figure>

The controls $u = (u_1, u_2, \dots, u_J)$ and the states $y = (y_1, y_2, \dots, y_J)$ are now split into subdomain quantities, and $\mathcal{U}$ is the set of admissible controls which is defined below. Note that we added additional unknowns function $g := (g_{ji})_{ji}$ on the interfaces which represent the normal derivative there that must match between the two subdomains, and also equations, because the matching conditions $\partial_{n_{ji}} y_j |_{\Gamma_{ji}} = \partial_{n_{ji}} y_i |_{\Gamma_{ji}}$ are replaced in (4) by the two equations stating that each normal trace must equal $g_{ji}$. Hence, the new unknowns $g_{ji}$, like the unknown solutions $y_j$, must now appear in the minimization problem formulation in the first line in (4) on the left.

The set of admissible controls $\mathcal{U}$ is defined as

$$ \mathcal{U} := \{ (u_j, g_{ji}) \mid u_j \in L^2(\Omega_j), g_{ji} \in H^{-1/2}(\Gamma_{ji}) \text{ and } g_{ji} = -g_{ij} \}. $$

We can prove that the problems (3) and (4) have the same feasible sets and are therefore equivalent. Note that we could also have added an unknown for the Dirichlet traces instead, or even both Dirichlet and Neumann which would lead to further new DD algorithms with the approach we describe below.
We next define a reduced problem associated to (4) by eliminating the subdomain solutions $y_j$ (like static condensation),

$$ \min_{u,g} \hat{\mathcal{J}}(u,g) := \sum_{j=1}^{J} \frac{1}{2} \|y_j(u_j, g_{ji}) - y_{\text{target}}\|_{L^2(\Omega_j)}^2 + \frac{\alpha}{2} \|u_j\|_{L^2(\Omega_j)}^2 \quad (5) $$

$$ \text{s.t. } y_j|_{\Gamma_{ji}} - y_i|_{\Gamma_{ji}} = 0, \quad (u,g) \in \mathcal{U}, $$

where the statically condensed solutions $y_j(u_j, g)$ satisfy

$$ \left\{ \begin{array}{ll} -\Delta y_j = u_j + f & \text{in } \Omega_j, \\ y_j = 0 & \text{on } \partial\Omega \cap \partial\Omega_j, \\ \partial_{n_{ji}} y_j|_{\Gamma_{ji}} = g_{ji} & \text{on each part of } \partial\Omega_j \setminus \partial\Omega = \cup_i \Gamma_{ji}. \end{array} \right. \quad (6) $$

Our new DD method is obtained when applying an augmented Lagrangian algorithm to the reduced optimization problem (10). In the augmented Lagrangian algorithm, one does not only add the Lagrange multiplier term to the cost functional to form the Lagrangian, but in addition also penalizes (some of) the constraints in the Lagrangian. The augmented Lagrangian algorithm is then given by an iteration, in which at each iteration step an unconstrained optimization problem is solved for fixed values of the Lagrange multipliers and penalization parameters, and based on the result, these multipliers and parameters are updated, see [32]. In the context of (5), the continuity constraints $y_j|_{\Gamma_{ji}} - y_i|_{\Gamma_{ji}} = 0$ are those we penalize. The augmented cost functional we use for the reduced optimization problem (5) is

$$ \hat{\mathcal{J}}^{+}(u, g; \lambda, \rho) = \hat{\mathcal{J}}(u, g) + \sum_{(i,j)} \int_{\Gamma_{ij}} (y_i(u_i, g) - y_j(u_j, g))\lambda_{ij} + \sum_{(i,j)} \frac{\rho_{ij}}{2} \|y_i(u_i, g) - y_j(u_j, g)\|_{L^2(\Gamma_{ij})}^2 \quad (7) $$

where the sum is taken so that each interior interface is considered only once (there is a slight abuse of notation in (7) since, if $\lambda_{ij} \in H^{-1/2}(\Gamma_{ji})$, the integral $\int_{\Gamma_{ij}}$ must be understood as a duality bracket). The augmented Lagrangian algorithm then solves approximately at each iteration for given values $\lambda = \lambda^k$ and $\rho = \rho^k$ the reduced minimization problem

$$ \min \hat{\mathcal{J}}^{+}(u, g; \lambda^{k}, \rho^{k}) \text{ s.t. } (u, g) \in \mathcal{U}. \quad (8) $$

This iteration defines our new DD method for the optimal control problem (3) decomposed into subdomains as in (5). This approach is summarized in Algorithm 2.1. It should be noted that a closely related approach can be found in [8], but the authors treat the variables $g_{ji}$ as implicit Lagrange multipliers for the continuity of the adjoint state over the interface $\Gamma$ (see [8, Remark 3.4.3]) and there is no convergence proof to the optimal solution.

**Remark 1** We have chosen to use an Augmented Lagrangian algorithm due to its property to limit the cases when the Hessian becomes ill-conditioned. Without the linear term $\int_{\Gamma_{ij}} (y_i(u_i, g) - y_j(u_j, g))\lambda_{ij}$, the algorithm becomes a quadratic penalization
<caption>**Algorithm 2.1:** New DD Algorithm based on Augmented Lagrangian</caption>

**Data:** Set initial tolerances $\{\eta_{ij}^0, \omega_{ij}^0\}_{ij} \subset (0,1)^2$, choose $\tau > 1$, and set $k = 1$. Choose an initial guess $u^0, g^0$, and initial parameters $\{\lambda_{ij}^0\}_{ij} \subset V_{\cap}^*$, $\{\rho_{ij}^0\}_{ij} \subset (1, \infty)$.

Compute $y^0(u^0, g^0)$ on each subdomain.

while $\sum_{\{ij, \Gamma_{ij} \neq \emptyset\}} \|y_i^{k-1} - y_j^{k-1}\|_{L^2(\Gamma)}^2 \geq \max_{ij} \eta_{ij}^k$ do
Solve approximately (8) for $u^k, g^k$ s.t. $\|\partial_{u,g} \hat{\mathcal{J}}^+(u^k, g^k)\|_{U^*} \leq \max_{ij} \omega_{ij}^k$;

for $i, j$ such that $\Gamma_{ij} \neq \emptyset$ do
if $\|y_i^k - y_j^k\|_{L^2(\Gamma)} \leq \eta_{ij}^k$ then
$\lambda_{ij}^{k+1} = \lambda_{ij}^k + \rho_{ij}^k(y_i^k - y_j^k)$; // update multiplier
$\rho_{ij}^{k+1} = \rho_{ij}^k$; // keep penalization
$\omega_{ij}^{k+1} = (\rho_{ij}^k)^{-1}\omega_{ij}^k$; // decrease tolerances
$\eta_{ij}^{k+1} = (\rho_{ij}^k)^{-1/2}\omega_{ij}^k$;

else
$\lambda_{ij}^{k+1} = \lambda_{ij}^k$; // keep multiplier
$\rho_{ij}^{k+1} = \tau \rho_{ij}^k$; // increase penalization
$\omega_{ij}^{k+1} = (\rho_{ij}^{k+1})^{-1}$; // update tolerances
$\eta_{ij}^{k+1} = (\rho_{ij}^{k+1})^{-1/2}$;

end

end

$k \leftarrow k + 1$

end

technique which forces the quadratic factors $\rho_{ij}$ to grow to infinity in order to ensure continuity. In this case, the Hessian will become ill-conditioned. However, with the addition of the linear term described above, we do not necessarily need to have $\rho$ growing to infinity, see the proof of Proposition 1 below.

## 3 Convergence analysis

In this section, we will analyze the convergence of Algorithm 2.1 toward the solution of (3). We will only treat the two subdomains case in detail, and comment on how these results can be extended to decomposition with more subdomains at the end of the current section.

Suppose we have two non-overlapping subdomains $\Omega_1, \Omega_2$ such that $\bar{\Omega} = \bar{\Omega}_1 \cup \bar{\Omega}_2$, $\Omega_1 \cap \Omega_2 = \emptyset$ and $\partial\Omega \cap \partial\Omega_j \neq \emptyset$. We also need the functional spaces

$$V_{\cap} := H_{00}^{\frac{1}{2}}(\Gamma), \quad U := L^2(\Omega_1) \times L^2(\Omega_2) \times H^{-\frac{1}{2}}(\Gamma), \qquad (9)$$
where $\Gamma_\cap = \partial\Omega_1 \cap \partial\Omega_2$ is the interface. The problem (5) with two subdomains then becomes

$$
\begin{equation} \tag{10}
\begin{aligned}
\min_{u,g} \hat{\mathcal{J}}(u,g) &:= \sum_{i=1}^{2} \frac{1}{2} \|y_i(u_i, g) - y_{\text{target}}\|_{L^2(\Omega_i)}^2 + \frac{\alpha}{2} \|u_i\|_{L^2(\Omega_i)}^2 \\
\text{s.t.} \quad & y_1|_{\Gamma_\cap} - y_2|_{\Gamma_\cap} = 0, \quad (u_1, u_2, g) \in \mathcal{U},
\end{aligned}
\end{equation}
$$

where the statically condensed solutions $y_i(u_i, g)$ satisfy

$$
\left\{
\begin{array}{ll}
-\Delta y_i = u_i + f & \text{in } \Omega_i, \\
y_i|_{\partial\Omega \cap \partial\Omega_i} = 0, & \text{(11)} \\
\partial_n y_i|_{\Gamma_\cap} = (-1)^{i+1} g, & i=1,2.
\end{array}
\right.
$$

The augmented Lagrangian associated to the reduced optimization problem (10) is
defined as

$$
\begin{equation}
\begin{split}
\min_{u,g} \hat{\mathcal{J}}^{+}(u,g; \lambda, \rho) ={}& \hat{\mathcal{J}}(u,g) \\
& + \int_{\Gamma} (y_1(u_1,g) - y_2(u_2,g))\lambda + \frac{\rho}{2} \|y_1(u_1,g) - y_2(u_2,g)\|_{L^2(\Gamma)}^2.
\end{split}
\tag{12}
\end{equation}
$$

The gradient of $\hat{\mathcal{J}}^+$ with respect to $u$ and $g$ can be computed by introducing the
Lagrangian associated to (10) and is given by

$$
\partial_{u,g} \hat{\mathcal{J}}^{+}(u, g; \lambda^{k}, \rho^{k}) = \begin{pmatrix} \alpha u_{1} - \tilde{p}_{1} \\ \alpha u_{2} - \tilde{p}_{2} \\ \tilde{p}_{2} |_{\Gamma} - \tilde{p}_{1} |_{\Gamma} \end{pmatrix},
$$

where $\tilde{p}_i$ are the adjoint variables that satisfy

$$
\left\{
\begin{array}{l}
-\Delta \tilde{p}_i + y_i = y_{\text{target}}, \\
\tilde{p}_i|_{\partial\Omega} = 0, \\
\partial_{n_i} \tilde{p}_i = (-1)^i (\lambda^k + \rho^k (y_1 - y_2)) \quad \text{on } \Gamma.
\end{array}
\right.
\tag{13}
$$

Note also that in this two subdomain case, $\eta_{ij}^k = \eta^k$ and $\omega_{ij}^k = \omega^k$.

**3.1 Existence of a unique solution**

Before studying the convergence of Algorithm 2.1, we first prove that (12) admits a minimizer and that, as $\rho \to +\infty$, the solution of the decoupled PDE constraint (11) at the optimum $(u, g)$ of (12) indeed satisfies the continuity at the interface. Note that theorem 1 remains true for a number $J > 2$ of subdomains.

**Theorem 1** For any $\lambda^k \in V^*_\cap$ and $\rho > 0$, the optimization problem (12) admits a unique solution $(\bar{u}, \bar{g})$. If the associated $(\overline{y_1}, \overline{y_2})$ satisfy (11), then, for any $\lambda^k \in L^2(\Gamma)$,
we have

$$ \|\overline{y_1} - \overline{y_2}\|_{L^2(\Gamma)}^2 \lesssim \frac{1}{\rho} + \frac{1}{\rho^2} \|\lambda^k\|_{L^2(\Gamma)}^2, \quad (14) $$

where the non explicitly mentioned constants do neither depend on $k$ nor on $\rho$.

*Proof* First note that the admissible set $\mathcal{U}$ is convex, and the map $(u_i, g) \mapsto y_i(u_i, g)$ is affine in $u_i$, linear in $g$ and injective. Therefore, we only need to prove the strict convexity of the associated function in the variables $(y, u)$,

$$ \mathcal{J}(y, u) := \frac{1}{2} \sum_{i=1}^{2} \left[ \|y_i - y_{\text{target}}\|_{L^2(\Omega_i)}^2 + \alpha \|u_i\|_{L^2(\Omega_i)}^2 \right] + \int_{\Gamma} (y_1 - y_2) \lambda^k + \frac{\rho^k}{2} \|y_1 - y_2\|_{L^2(\Gamma)}^2. $$

$\mathcal{J}$ is twice differentiable, and its Hessian matrix is clearly positive definite. This proves that $\mathcal{J}$ is strictly convex and [34, Theorem 1.46] then gives existence and uniqueness of the solution.

Let $(\bar{u}, \bar{g})$ be the unique solution to (12). We then have

$$ \mathcal{J}(\bar{y}(\bar{u}, \bar{g}), \bar{u}) \le \mathcal{J}(y(u, g), u), \text{ for all admissible } u, g. \quad (15) $$

For any fixed $u \in L^2(\Omega)$, let $y_\Omega \in H_0^1(\Omega)$ be the unique solution to $-\Delta y_\Omega = u$ and $g_\Omega$ such that the restrictions $y_{\Omega,i} := y_\Omega|_{\Omega_i}$ satisfy $\partial_{\mathbf{n}_i} y_i|_{\Gamma_\cap} = (-1)^{i+1}g_\Omega$. Since $y_\Omega \in H^1(\Omega)$, we have $y_{\Omega,1} = y_{\Omega,2}$ on $\Gamma_\cap$, and (15) then yields

$$ \begin{aligned} \int_{\Gamma} (\overline{y_1} - \overline{y_2}) \lambda^k &+ \frac{\rho}{2} \|\overline{y_1} - \overline{y_2}\|_{L^2(\Gamma)}^2 \le \mathcal{J}(y_\Omega(u, g_\Omega), u) \\ & = \frac{1}{2} \left[ \|y_\Omega - y_{\text{target}}\|_{L^2(\Omega)}^2 + \alpha \|f\|_{L^2(\Omega)}^2 \right] =: C. \end{aligned} $$

Young's inequality together with the Cauchy-Schwarz inequality then gives

$$ 2 \left| \int_{\Gamma} (\overline{y_1} - \overline{y_2}) \lambda^k \right| \le 2 \|\lambda^k\|_{L^2(\Gamma_\cap)} \|\overline{y_1} - \overline{y_2}\|_{L^2(\Gamma_\cap)} \le \frac{2}{\rho} \|\lambda^k\|_{L^2(\Gamma_\cap)}^2 + \frac{\rho}{2} \|\overline{y_1} - \overline{y_2}\|_{L^2(\Gamma_\cap)}^2. $$

Combining the previous inequalities, we get

$$ \rho \|\overline{y_1} - \overline{y_2}\|_{L^2(\Gamma)}^2 \le 2C - 2 \int_{\Gamma} (\overline{y_1} - \overline{y_2}) \lambda^k \le 2C + \frac{2}{\rho} \|\lambda^k\|_{L^2(\Gamma_\cap)}^2 + \frac{\rho}{2} \|\overline{y_1} - \overline{y_2}\|_{L^2(\Gamma_\cap)}^2, $$

from which we finally infer the estimate (14). $\square$

## 3.2 Convergence of the augmented Lagrangian

We now prove convergence of our new domain decomposition Algorithm 2.1. For ease of presentation, we denote by $x := (u_1, u_2, g)$. As in [40] (see also references therein), we rely on a so-called asymptotic KKT (AKKT) conditions (see e.g. [40, Definition 5.2]). We will prove that Algorithm 2.1 implies that the generated sequence complies with this AKKT condition and that this yields convergence to the optimal solution
of (10). As exposed in [40], this approach is common for analyzing the convergence of algorithms based on the augmented Lagrangian.

Before introducing AKKT conditions, we first give the (standard) KKT conditions for (10). We introduce the Lagrangian associated to (10) which reads

$$ \mathcal{L}(x, \lambda) = \hat{\mathcal{J}}(x) + \langle \lambda, c_o(x) \rangle_{V_{\cap}^*, V_{\cap}}, \quad (16) $$

where $c_o: x = (u_1, u_2, g) \in \mathcal{U} \mapsto y_1(u_1, g)|_{\Gamma} - y_2(u_2, g)|_{\Gamma} \in V_{\cap}$ is the linear map giving the continuity constraint and $y_i$ satisfies (11). We recall that some $x \in \mathcal{U}$ is a KKT point of (10) if

$$ \exists \lambda \in V_{\cap}^* \text{ such that } \partial_x \mathcal{L}(x, \lambda) = 0 \text{ and } c_o(x) = 0. \quad (17) $$

To compute the derivative of $\mathcal{L}$ with respect to $x$, we introduce the linear operator $M(x) : V_{\cap}^* \to \mathcal{U}^*$ such that $M(x)^* = \partial_x c_o(x)$. Since $y_i$ are affine functions of $x$, the linear map $\partial_x c_o(x)$ does not depend on $x$ and we use the notation $M := M(x)$. The differential of $\mathcal{L}$ with respect to $x$ can then be written as

$$ \partial_x \mathcal{L}(x, \lambda)[\delta x] = \partial_x \hat{\mathcal{J}}(x)[\delta x] + \langle \lambda, M^* \delta x \rangle_{V_{\cap}^*, V_{\cap}} = \langle \partial_x \hat{\mathcal{J}}(x) + M\lambda, \delta x \rangle_{\mathcal{U}^*, \mathcal{U}}. $$

We can then recast (17) so that $x \in \mathcal{U}$ is a KKT point for (10) if

$$ \exists \lambda \in V_{\cap}^* \text{ such that } \partial_x \hat{\mathcal{J}}(x) + M\lambda = 0 \text{ and } c_o(x) = 0. \quad (18) $$

It is usual to ask for $M^*$ to be surjective; this condition is called constraint qualification. As shown in the Appendix, for a two subdomain decomposition, one has **Theorem 2** The map $M^*: \mathcal{U} \to V_{\cap}$ is surjective and linear.

This constraint qualification will now serve to analyze the convergence of iterates $(x^k, \lambda^k)$ which we suppose to be AKKT points, defined as follows:

**Definition 1** We say that a feasible$^3$ point $x \in \mathcal{U}$ respects the AKKT condition, if there are sequences $x_k \to x$ in $\mathcal{U}$ and $(\lambda^k) \subset V_{\cap}^*$ such that

$$ \partial_x \mathcal{L}(x^k, \lambda^k) = \partial_x \hat{\mathcal{J}}(x^k) + M\lambda^k \to 0 \text{ as } k \to +\infty \text{ in } \mathcal{U}^*. \quad (19) $$

This convergence property implies boundedness of the sequence of multipliers:

**Lemma 3** Suppose (19) is verified at some point $x^*$ (not necessarily feasible), and let $(x^k, \lambda^k)$ be the associated sequence. Then $(\lambda^k)$ is bounded.

**Proof** First, note that $x \in \mathcal{U} \mapsto \partial_x \hat{\mathcal{J}}(x) \in \mathcal{U}^*$ and $\lambda \in V_{\cap}^* \mapsto M\lambda \in \mathcal{U}^*$ are strongly continuous. Due to theorem 2 and [10, Théorème II.2.20], we have the estimate

$$ \forall \lambda \in V_{\cap}^*, \quad \|\lambda\|_{V_{\cap}^*} \lesssim \|M\lambda\|_{\mathcal{U}^*}. \quad (20) $$

<footnote>$^3$A point $x$ is called feasible if $c_o(x) = 0$.</footnote>
Therefore, for any $n \in \mathbb{N}$, we get

$$ \begin{aligned} \|\lambda^n - \lambda^0\|_{V_{\cap}^*} &\lesssim \|M\lambda^n - M\lambda^0\|_{U^*} \\ &\lesssim \|\partial_x \hat{\mathcal{J}}(x^n) + M\lambda^n\|_{U^*} + \|\partial_x \hat{\mathcal{J}}(x^0) + M\lambda^0\|_{U^*} + \|\partial_x \hat{\mathcal{J}}(x^n) - \partial_x \hat{\mathcal{J}}(x^0)\|_{U^*} \\ &\le C. \end{aligned} $$

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

We now prove that Algorithm 2.1 produces a sequence $(x^k, \lambda^k)_k$ which is an AKKT point.

**Proposition 1** Let $(x^k, \lambda^k, \rho^k)$ be the sequence generated by Algorithm 2.1, and suppose that $x^k \to x^*$ in $U$. Then $x^*$ is an AKKT point.

*Proof* To prove that $x^*$ is an AKKT point, we start by noting that the cost function $\hat{\mathcal{J}}^+(x; \lambda, \rho)$ defined in (12) can be written as

$$ \hat{\mathcal{J}}^+(x; \lambda, \rho) = \hat{\mathcal{J}}(x) + \langle \lambda, c_o(x) \rangle_{V_{\cap}^*, V_{\cap}} + \frac{\rho}{2} \|c_o(x)\|_{L^2(\Gamma)}^2, $$

and thus $\partial_x \hat{\mathcal{J}}^+(x; \lambda, \rho) = \partial_x \hat{\mathcal{J}}(x) + M(\lambda + \rho c_o(x))$. Therefore, in Algorithm 2.1, $x^{k+1}$ complies with

$$ \|\partial_x \hat{\mathcal{J}}^+(x^{k+1}; \lambda^k, \rho^k)\|_{U^*} = \|\partial_x \hat{\mathcal{J}}(x^{k+1}) + M\lambda^{k+1}\|_{U^*} \le \omega^k. \quad (21) $$

As proved in [48, Lemma 3.1.1], $\lim \omega^k = \lim \eta^k = 0$. As a result, $x^*$ verifies (19) and it remains to prove that $x^*$ is a feasible point.

To prove that $c_o(x^k) \to 0$, we consider two cases: $(\rho^k)$ bounded and $\rho^k \to +\infty$.

Case 1: Suppose $(\rho^k)$ is bounded. Due to how $(\rho^k)$ is updated, this implies that there exists $k_0 \in \mathbb{N}$ such that, for all $k \ge k_0$, $\rho^k = \rho^{k_0}$. Therefore, for all $k \ge k_0$, we have $\|y_1(x^k) - y_2(x^k)\|_{L^2(\Gamma)} \le \eta^k$. Since $\eta^k \to 0$, this implies $y_1(x^k) - y_2(x^k) \to 0$ in $L^2(\Gamma)$.

Case 2: Let us now suppose $\rho^k \to +\infty$. We recall that $\lambda^{k+1} = \lambda^k + \rho^k(y_1^k - y_2^k)$. Condition (21) can then be recast as

$$ \left\| \left( \partial_x \hat{\mathcal{J}}(x^{k+1}) + M\lambda^k \right) + \rho^k M (y_1(x^k) - y_2(x^k)) \right\|_{U^*} \le \omega^k. $$

Therefore, multiplying by $(\rho^k)^{-1}$, we obtain

$$ \|M(y_1(x^k) - y_2(x^k))\|_{U^*} \le \omega^k (\rho^k)^{-1} + (\rho^k)^{-1} \left\| \partial_x \hat{\mathcal{J}}(x^{k+1}) + M\lambda^k \right\|_{U^*}. $$

Since $(\partial_x \hat{\mathcal{J}}(x^k))$ and $(\lambda^k)$ are bounded, we find $M(y_1(x^k) - y_2(x^k)) \to 0$. Using now (20), one finds $y_1(x^k) - y_2(x^k) \to 0$ in $V_{\cap}^*$.

Since the map $x \mapsto c_o(x)$ is continuous and we assumed that $x^k \to x^*$, we have $c_o(x^k) \to c_o(x^*)$ in $V_{\cap}$. For both cases considered above, we end up having $c_o(x^k) = y_1(x^k) - y_2(x^k) \to 0$ in $L^2(\Gamma)$ (Case 1) or in $V_{\cap}^*$ (Case 2). We recall that we have
$V_\cap \subset L^2(\Gamma) \subset V_\cap^*$ with continuous embedding and since $c_o(x^k)$ converges to 0 in either $L^2(\Gamma)$ or $V_\cap^*$, the uniqueness of the limit in these spaces ensure that $c_o(x^*) = 0$. □

We now prove that $x^*$ is the optimal solution of (10).

**Proposition 2** Let $(x^k, \lambda^k)$ be the sequence generated by Algorithm 2.1. Suppose that $x^k \to x^*$ in $U$, then $x^*$ is the optimal solution of (10) with associated multiplier $\lambda^*$ which is the weak-limit of $(\lambda^k)$.

*Proof* Proposition 1 gives that $x^*$ is an AKKT point. Since $V_\cap^*$ is reflexive, its unit ball is weakly compact. Therefore, theorem 3 implies that there exists $\lambda^*$ such that $\lambda^n \to \lambda^*$ in $V_\cap^*$. Therefore, taking the limit in (19), and since $x^*$ is feasible, $(x^*, \lambda^*)$ is such that

$$ \partial_x \mathcal{L}(x^*, \lambda^*) = \partial_x \hat{\mathcal{J}}(x^*) + M\lambda^* = 0, c_o(x^*) = 0, $$

which is exactly the KKT condition (18). Since the optimization problems (10) and the initial one (3) are equivalent, we obtain that $\lambda^*$ is given by (24) and that $x^*$ is the optimal solution of (10) which is unique since the latter is convex. □

**Remark 2** The multiplier $\lambda^*$, associated to the continuity constraint of the original optimization problem (10), is also directly linked to the continuity constraint of the adjoint variable.

The first order necessary conditions of optimality for (10) are given by

$$ \left\{ \begin{array}{l} \alpha u_i - p_i = 0 \quad \text{in } \Omega_i, \\ p_2|_{\Gamma} - p_1|_{\Gamma} = 0, \end{array} \right. , \quad i = 1, 2, \quad \left\{ \begin{array}{l} -\Delta p_i + y_i = y_{\text{target}}, \\ p_i|_{\partial \Omega} = 0, \\ \partial_{\mathbf{n}_i} p_i = (-1)^{i+1} \lambda^* \quad \text{on } \Gamma, \end{array} \right. \quad (22) $$

where $p_i$ are adjoint variables and $y_i = y_i(u_i, g)$ are solutions to (11) that also comply with $y_1(u_1, g)|_{\Gamma} - y_2(u_2, g)|_{\Gamma} = 0$.

Introducing now $y, p \in H_0^1(\Omega)$ and $u \in L^2(\Omega)$ such that

$$ y|_{\Omega_i} = y_i, \quad p|_{\Omega_i} = p_i, \quad u|_{\Omega_i} = u_i, $$

we get that $(y, p, u)$ satisfy

$$ \left\{ \begin{array}{l} -\Delta y = f + u \quad \text{in } \Omega \\ y \in H_0^1(\Omega), \end{array} \right. , \quad \left\{ \begin{array}{l} -\Delta p = -(y - y_{\text{target}}) \quad \text{in } \Omega, \\ p \in H_0^1(\Omega), \end{array} \right. , \quad \alpha u - p = 0 \quad \text{in } \Omega. \quad (23) $$

We emphasize that (23) are exactly the first order conditions of optimality of the undecomposed optimization problem (3). We thus obtain that

$$ \lambda^* = (-1)^{i+1} \partial_{\mathbf{n}_i} p|_{\Gamma} \quad (24) $$

where $p$ satisfies the coupled system (23).

**Remark 3 (Extension to more than 2 subdomains)** The extension of the convergence results to more subdomains mainly consists in proving that theorem 2 holds when more subdomains are involved. As shown in the appendix, the proof consists in proving that $M$ is injective (which can be easily adapted to a more complex decomposition) and with closed range. This second part seems more challenging for several
subdomains and depends on how the decomposition is made. For example, the proof of theorem 4 can be straightforwardly adapted to stripwise decomposition like those used in the numerical experiments from Section 4, mainly because we do not have subdomains with Neumann boundary conditions only. Furthermore, as long as one supposes $M^*$ surjective, all the results leading to proposition 2 can be easily adapted to a context with more subdomains. This extension to more subdomains can be seen more easily when one carries a similar analysis we did in section 3 once the problem (12) is discretized; see appendix B.

## 3.3 Fourier analysis

To get more insight into the convergence behavior of our new DD Algorithm 2.1, we study its rate of convergence in the case where the penalization parameter $\rho$ is kept constant during the iterations, and the domain is unbounded, $\Omega := \mathbb{R}^2$, with $(x_1, x_2)$ a system of coordinates on $\mathbb{R}^2$, in order to use Fourier analysis, which can give very detailed insight into the convergence of such algorithms, see e.g. [20, 24]. Let the subdomains $\Omega_1 := (-\infty, 0) \times \mathbb{R}$, $\Omega_2 := (0, +\infty) \times \mathbb{R}$, and the interface $\Gamma := \{0\} \times \mathbb{R}$. In this setting, we can study how the iterations $\{\lambda^k\}$ are converging. The necessary and sufficient optimality conditions for (12) are

$$ 
\left\{ 
\begin{array}{ll} 
-\Delta y_i^k = f + \alpha^{-1} p_i^k & \text{in } \Omega_i, \\ 
y_i^k|_{\partial\Omega} = 0, \\ 
\partial_{\mathbf{n}_i} y_i^k|_{\Gamma_\cap} = (-1)^{i+1} g^k, & i=1,2 
\end{array} 
\right. 
\quad 
\left\{ 
\begin{array}{l} 
-\Delta p_i^k + y_i^k = y_{\text{target}}, \\ 
p_i^k|_{\partial\Omega} = 0, \\ 
\partial_{\mathbf{n}_i} p_i^k = (-1)^i (\lambda^k + \rho(y_1^k - y_2^k)) \text{ on } \Gamma, \\ 
p_1^k = p_2^k \text{ on } \Gamma. 
\end{array} 
\right. 
\qquad (25) 
$$

Since we are interested in the error in the new DD Algorithm 2.1, without loss of generality, we will suppose that $f \equiv 0$ and $y_{\text{target}} \equiv 0$, so that we can study convergence to zero, i.e. the error equations.

We assume that $\lambda^0$ has enough regularity so that its Fourier transform along the vertical axis is well-defined. Let $\omega \in \mathbb{R}$ be the Fourier variable and $(\hat{y}, \hat{p})$ be the Fourier transform in the $x_2$ direction of $(y, p)$ along the interface. The Fourier transformed optimality conditions (25) then become

$$ 
\left\{ 
\begin{array}{ll} 
(-\partial_{x_1}^2 + \omega^2) \hat{y}_i^k = \alpha^{-1} \hat{p}_i^k & \text{in } \Omega_i, \\ 
\hat{y}_i^k|_{\partial\Omega} = 0, \\ 
\partial_{x_1} \hat{y}_i^k|_{\Gamma_\cap} = \hat{g}^k, & i=1,2 
\end{array} 
\right. 
\quad 
\left\{ 
\begin{array}{l} 
(-\partial_{x_1}^2 + \omega^2) \hat{p}_i^k + \hat{y}_i^k = 0, \\ 
\hat{p}_i^k|_{\partial\Omega} = 0, \\ 
\partial_{x_1} \hat{p}_i^k = - (\hat{\lambda}^k + \rho(\hat{y}_1^k - \hat{y}_2^k)) \text{ on } \Gamma, \\ 
\hat{p}_1^k = \hat{p}_2^k \text{ on } \Gamma. 
\end{array} 
\right. 
\qquad (26) 
$$
In order to solve this system of ordinary differential equations, we compute $\partial_{x_1}^4 \hat{y}_i^k$ and substitute,

$$ 
\begin{aligned}
\partial_{x_1}^4 \hat{y}_i^k &= \omega^2 \partial_{x_1}^2 \hat{y}_i^k - \alpha^{-1} \partial_{x_1}^2 \hat{p}_i^k \\
&= \omega^2 \partial_{x_1}^2 \hat{y}_i^k - \alpha^{-1} (\omega^2 \hat{p}_i^k + \hat{y}_i^k) \\
&= \omega^2 \partial_{x_1}^2 \hat{y}_i^k - \alpha^{-1} \hat{y}_i^k - \omega^2 (-\partial_{x_1}^2 \hat{y}_i^k + \omega^2 \hat{y}_i^k) \\
&= 2\omega^2 \partial_{x_1}^2 \hat{y}_i^k - (\omega^4 + \alpha^{-1}) \hat{y}_i^k.
\end{aligned}
$$

Using the boundary conditions at infinity, the solutions are thus of the form

$$ 
\begin{aligned}
\hat{y}_1^k &= C_1^k \exp \left( x_1 \sqrt{i\alpha^{-\frac{1}{2}} + \omega^2} \right) + C_3^k \exp \left( x_1 \sqrt{-i\alpha^{-\frac{1}{2}} + \omega^2} \right), \\
\hat{y}_2^k &= C_2^k \exp \left( -x_1 \sqrt{i\alpha^{-\frac{1}{2}} + \omega^2} \right) + C_4^k \exp \left( -x_1 \sqrt{-i\alpha^{-\frac{1}{2}} + \omega^2} \right),
\end{aligned}
$$

where $C_j^k$ are constants determined by the transmission conditions. This implies for the adjoint variables

$$ 
\begin{aligned}
\hat{p}_1^k &= i\alpha^{\frac{1}{2}} \left( -C_1^k \exp \left( x_1 \sqrt{i\alpha^{-\frac{1}{2}} + \omega^2} \right) + C_3^k \exp \left( x_1 \sqrt{-i\alpha^{-\frac{1}{2}} + \omega^2} \right) \right), \\
\hat{p}_2^k &= i\alpha^{\frac{1}{2}} \left( -C_2^k \exp \left( -x_1 \sqrt{i\alpha^{-\frac{1}{2}} + \omega^2} \right) + C_4^k \exp \left( -x_1 \sqrt{-i\alpha^{-\frac{1}{2}} + \omega^2} \right) \right).
\end{aligned}
$$

The transmission conditions at $x_1 = 0$ then yield

$$ 
\begin{aligned}
\partial_{x_1} \hat{y}_1^k |_{x_1=0} &= C_1^k \sqrt{i\alpha^{-\frac{1}{2}} + \omega^2} + C_3^k \sqrt{-i\alpha^{-\frac{1}{2}} + \omega^2} = \hat{g}^k, \\
\partial_{x_1} \hat{y}_2^k |_{x_1=0} &= -C_2^k \sqrt{i\alpha^{-\frac{1}{2}} + \omega^2} - C_4^k \sqrt{-i\alpha^{-\frac{1}{2}} + \omega^2} = \hat{g}^k, \\
\partial_{x_1} \hat{p}_1^k |_{x_1=0} &= -i\alpha^{\frac{1}{2}} \left( C_1^k \sqrt{i\alpha^{-\frac{1}{2}} + \omega^2} - C_3^k \sqrt{-i\alpha^{-\frac{1}{2}} + \omega^2} \right) = -\hat{\lambda}^{k+1}, \\
\partial_{x_1} \hat{p}_2^k |_{x_1=0} &= i\alpha^{\frac{1}{2}} \left( C_2^k \sqrt{i\alpha^{-\frac{1}{2}} + \omega^2} - C_4^k \sqrt{-i\alpha^{-\frac{1}{2}} + \omega^2} \right) = -\hat{\lambda}^{k+1}.
\end{aligned}
$$

Denoting by $A := \sqrt{i\alpha^{-\frac{1}{2}} + \omega^2}$ and $B := \sqrt{-i\alpha^{-\frac{1}{2}} + \omega^2}$, we thus need to solve the linear system

$$ 
\begin{pmatrix} 1 & 0 & 1 & 0 \\ 0 & -1 & 0 & -1 \\ 1 & 0 & -1 & 0 \\ 0 & 1 & 0 & -1 \end{pmatrix} \begin{pmatrix} AC_1^k \\ AC_2^k \\ BC_3^k \\ BC_4^k \end{pmatrix} = \begin{pmatrix} \hat{g}^k \\ \hat{g}^k \\ i\alpha^{-\frac{1}{2}}\hat{\lambda}^{k+1} \\ -i\alpha^{-\frac{1}{2}}\hat{\lambda}^{k+1} \end{pmatrix}.
$$

The inverse of the system matrix being simply half its transpose, we get

$$ 
\begin{aligned}
AC_1^k &= -AC_2^k = \frac{1}{2}(\hat{g}^k + i\alpha^{-\frac{1}{2}}\hat{\lambda}^{k+1}), \\
BC_3^k &= -BC_4^k = \frac{1}{2}(\hat{g}^k - i\alpha^{-\frac{1}{2}}\hat{\lambda}^{k+1}).
\end{aligned}
$$
<figure><img src="image_5.png" /><figcaption>Fig. 3: $|R|$ for $\rho = 3$, $\alpha = 1e-5$ (Left) and $\alpha = 1$ (Right).</figcaption></figure>

Note in particular that $C_3^k = -C_4^k$ and $C_1^k = -C_2^k$, and furthermore,

$$ 
\begin{aligned}
\hat{y}_1^k|_{x_1=0} - \hat{y}_2^k|_{x_1=0} &= C_1^k + C_3^k - (C_2^k + C_4^k) \\
&= 2(C_1^k + C_3^k) \\
&= (A^{-1} + B^{-1})\hat{g}^k + (A^{-1} - B^{-1})i\alpha^{-\frac{1}{2}}\hat{\lambda}^{k+1},
\end{aligned}
$$
$$ 
\begin{aligned}
\hat{p}_1^k|_{x_1=0} - \hat{p}_2^k|_{x_1=0} &= -i\alpha^{\frac{1}{2}}(C_1^k - C_3^k) + i\alpha^{\frac{1}{2}}(C_2^k - C_4^k) \\
&= -i\alpha^{\frac{1}{2}}(C_1^k - C_3^k - C_2^k + C_4^k) \\
&= 2i\alpha^{\frac{1}{2}}(C_3^k - C_1^k) \\
&= i\alpha^{\frac{1}{2}}(B^{-1} - A^{-1})\hat{g}^k + (A^{-1} + B^{-1})\hat{\lambda}^{k+1}.
\end{aligned}
$$

Since $\hat{p}_1^k|_{x_1=0} - \hat{p}_2^k|_{x_1=0} = 0$, this implies $\hat{g}^k = i\alpha^{-\frac{1}{2}}\frac{(A+B)}{A-B}\hat{\lambda}^{k+1}$ from which we infer

$$ \hat{y}_1^k|_{x_1=0} - \hat{y}_2^k|_{x_1=0} = 4i\alpha^{-\frac{1}{2}}(A-B)^{-1}\hat{\lambda}^{k+1}. $$

Denote by $D := 4i\alpha^{-\frac{1}{2}}(A-B)^{-1}$. Since in our new DD algorithm $\hat{\lambda}^k = \hat{\lambda}^{k+1} - \rho(\hat{y}_1^k|_{x_1=0} - \hat{y}_2^k|_{x_1=0})$, substituting the expression of $\hat{y}_1^k|_{x_1=0} - \hat{y}_2^k|_{x_1=0}$ gives $\hat{\lambda}^k = (1 - \rho D)\hat{\lambda}^{k+1}$, and we thus obtain for the convergence factor

$$ R := \frac{\hat{\lambda}^{k+1}}{\hat{\lambda}^k} = (1 - \rho D)^{-1} = \left( 1 - \frac{4i\rho}{\alpha^{1/2} \left( \sqrt{\omega^2 + i\alpha^{-1/2}} - \sqrt{\omega^2 - i\alpha^{-1/2}} \right)} \right)^{-1}, \quad (27) $$

where we substituted the expressions for $D$. Direct computations show that

$$ \lim_{\alpha \to 0} R = 0, \quad R = \frac{i\alpha^{1/2}}{4\rho} \left( \sqrt{\omega^2 + i\alpha^{-1/2}} - \sqrt{\omega^2 - i\alpha^{-1/2}} \right) + O\left(\frac{1}{\rho^2}\right). $$

Note that, for all $\alpha \ge 0$ and any frequency $\omega$, $|R|$ is strictly smaller than 1 as soon as $\rho$ is large enough which indicates geometric convergence of the new DD algorithm. We show a plot of $|R|$ for chosen parameters in fig. 3. We see that convergence is rather fast,
even for low frequencies $\omega$, and the method is a smoother: high frequencies converge extremely fast. Such behavior is also confirmed by our numerical results from Section 4 which show that the convergence of our decompose-then-optimize method is robust with respect to meshsize along with the number of subdomains. This Fourier analysis also indicates that these properties are going to hold regardless the discretization used.

# 4 Numerical experiments

We now test our algorithm on two examples. We will mainly pay attention to how the error evolves with the number of outer iterations and when the mesh is refined or the number of domains increases. An outer iteration here should be understood as the number of times one solves (8) and updates $\rho$ or $\lambda$; see the **while** loop in Algorithm 2.1. We discretized the Laplacian using Q1 Finite Elements on a structured cartesian mesh with step size $h$ in both $x_1$ and $x_2$ directions. The minimization problems in each iteration of our new DD algorithm are solved using the **quadprog** routine of MATLAB which uses an interior-point algorithm, and we used $\tau = 5$, $\rho_{ij}^0 = 3$, and $\lambda_{ij}^0$ is a random initial guess with values between 0 and 10. The options for **quadprog** are the default ones, namely 1e-8 for the **ConstraintTolerance**, except the **OptimalityTolerance** parameter which starts at $10^{-3}$ and is decreased with the iteration of the augmented Largangian algorithm until it reaches $10^{-10}$. This corresponds to a decrease of $\omega_{ij}^k$ appearing in Algorithm 2.1 but with a threshold placed on this decrease. All computations were performed in double precision floating-point arithmetic. Note that we use **quadprog** to solve the optimization problem (8) globally, including all subdomains in a single optimization procedure. The matrices defining the cost and the constraints are built by subdomain and combined in order to define the whole optimization problem. The code used for the computations is available online$^4$.

## 4.1 Badly insulated room

Our first example is a badly insulated room, for which a floor heating should be designed, and we study the convergence of our new DD method when only 5 iterations are performed. The room geometry is $\Omega := (-1, 1) \times (0, 1)$, with non-insulated walls that are at temperature zero, modeled by homogeneous Dirichlet conditions, and we want to see how a floor heating system would have to operate to heat it up to 20 degrees hence $y_{\text{target}}(x_1, x_2) = 20$. The room has in addition an open window at the bottom left where cold air acts as a heat sink, and similarly a skylight on the right in the middle, see Figure 4 (top left). We model this with $f(x_1, x_2) = -3000e^{-50((x_1-0.5)^2+(x_2-0.5)^2)} - 3000 \times \mathbf{1}_{[-0.75,-0.25]\times[0,0.25]}$, where $\mathbf{1}_O$ is the indicator function of the set $O$.

We decompose the room into two subdomains $\Omega_1 := (-1, 0) \times (0, 1)$ and $\Omega_2 := (0, 1) \times (0, 1)$, and apply our new domain decomposition algorithm in order to compute the optimal floor heating. We show in Figure 4 at the top the optimal control found for the floor heating by our new DD algorithm (middle), and the achieved room temperature (right), when using as penalization $\alpha = 1e-5$. We see that such a badly insulated room is not at all suitable for floor heating which needs to be evenly distributed; a more classical configuration with radiators along the walls, and especially

<footnote n="4">4 https://gitlab.osureunion.fr/avieira/liquofeti</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 4:** Badly insulated room example: heat sinks (top left), best floor heating found by our new DD algorithm (top middle), temperature distribution achieved (top right). Bottom left: convergence of the new DD algorithm for the control $u$ for different mesh sizes (log scale for the ordinate axis); bottom right: corresponding results for the solution $y$ (log scale for the ordinate axis).</figcaption></figure>

in front of the window would be much more suitable. Floor heating systems are only good for well insulated buildings! In the bottom row, we show the convergence of our new DD algorithm, on the left for the controls $u_j$ and on the right for the solutions $y_j$. The error here is the difference between the iterates of our new DD algorithm and the solution computed on the whole domain directly, measured in the infinity norm in volume. We see that convergence is very fast: the error is reduced by 6 orders of magnitude in five iterations, and convergence is robust in the mesh size $h$: when the mesh is refined, convergence gets actually a little faster before approaching a limit.

We next increase the number of subdomains, using a decomposition into equal vertical strips $\Omega_j = (x_j, x_{j+1}) \times (0,1)$, $x_j = -1 + 2j/J$, $j = 0, \dots, J$. Note that since we are taking vertical strips, there is always a part of the upper and lower boundary of the whole domain included in each subdomains, thus the optimal solution is never trivial in each of these (see Figure 4). We show in Figure 5 the errors for $J = 2, 4, 8, 16$ subdomains and meshsize $h = 1/64$. We see that convergence is still very fast and does not depend on the number of subdomains.
<figure><img src="image_7.png" /><figcaption>Fig. 5: Convergence of the new DD algorithm for the control u (left) and solution y (right) for different numbers of subdomains for the badly insulated room example. Both figures are with log scale for the ordinate axis.</figcaption></figure>

## 4.2 Convergence behavior with respect to $\alpha$

We now test the convergence properties of our new DD algorithm as $\alpha$ varies. To do this, we consider a model problem inspired by [27], namely

$$y_{\text{target}}(x_1, x_2) := C \sin(l\pi x_1) \sin(l\pi x_2), \text{ and } f(x_1, x_2) := 2Cl^2\pi^2 \sin(l\pi x_1) \sin(l\pi x_2),$$

for different parameters $C \in \mathbb{R}$, and $l \in \mathbb{N}$. The optimal solution is $u^* = 0$ and $y^* = y_{\text{target}}$, for all admissible choices of parameters $C, l, \alpha$. We solve this problem with our new DD Algorithm 2.1 using a decomposition into equal vertical strips as in the first physical example, and convergence parameters $\max_{ij} \omega_{ij}^k = \max_{ij} \eta_{ij}^k = 10^{-10}$, meaning that we stop the Augmented Lagrangian iterations when for all meaningful $i$, $\|y_i - y_{i+1}\|_{L^2(\{x_i\} \times (0,1))} \le 10^{-10}$ where $\{x_i\} \times (0,1)$ denotes each interface between two vertical strips. This example illustrates the convergence with respect to the outer iterations since the optimal solution is known, and also shows how the algorithm behaves when $\alpha$ is changing.

Results for different values of $\alpha, h$ and the number of subdomains are shown in table 1 and fig. 6. We see that our new DD algorithm remains robust and fast when changing the parameter $\alpha$, the stepsize $h$ and the number of subdomains. The number of outer iterations needed to achieve the prescribed precision is roughly constant in all our tests. Note nevertheless that, for large $\alpha$, the rate slightly changes with more subdomains.

## 5 Conclusions

We introduced a new concept for designing domain decomposition methods to solve optimal control problems: instead of following the optimize then decompose approach, we followed the decompose then optimize approach. Using then an augmented Lagrangian algorithm for the decomposed problem (or any other optimization technique) leads to new types of domain decomposition methods for such problems. We
<table><caption>**Table 1:** Errors on the optimal solution for the second example using different values of $\alpha$ and number of subdomains. The parameters were fixed as $l = 2$, $C = 1$.</caption><thead><tr><th>#subd.</th><th>h</th><th>α</th><th>Err. on y</th><th>Err. on u</th><th>Err. on g</th><th>Cost</th><th># it.</th></tr></thead><tbody><tr><td rowspan="10">2</td><td rowspan="5">1/64</td><td>1.0e+01</td><td>1.03e-09</td><td>1.81e-12</td><td>4.02e-08</td><td>4.18e-21</td><td>6</td></tr><tr><td>1.0e+00</td><td>1.07e-09</td><td>1.87e-11</td><td>4.02e-08</td><td>4.49e-21</td><td>6</td></tr><tr><td>1.0e-02</td><td>5.78e-10</td><td>8.58e-10</td><td>4.02e-08</td><td>1.14e-21</td><td>6</td></tr><tr><td>1.0e-04</td><td>1.50e-09</td><td>4.69e-08</td><td>4.03e-08</td><td>2.06e-21</td><td>5</td></tr><tr><td>1.0e-06</td><td>5.30e-09</td><td>1.62e-06</td><td>4.14e-08</td><td>2.67e-21</td><td>4</td></tr><tr><td rowspan="5">1/256</td><td>1.0e+01</td><td>1.11e-09</td><td>1.97e-12</td><td>4.92e-11</td><td>4.91e-21</td><td>6</td></tr><tr><td>1.0e+00</td><td>1.12e-09</td><td>1.98e-11</td><td>4.93e-11</td><td>4.98e-21</td><td>6</td></tr><tr><td>1.0e-02</td><td>7.35e-10</td><td>1.09e-09</td><td>4.65e-11</td><td>1.83e-21</td><td>6</td></tr><tr><td>1.0e-04</td><td>1.43e-09</td><td>4.45e-08</td><td>7.58e-11</td><td>1.84e-21</td><td>5</td></tr><tr><td>1.0e-06</td><td>4.68e-09</td><td>1.48e-06</td><td>6.34e-10</td><td>2.64e-21</td><td>4</td></tr><tr><td rowspan="10">4</td><td rowspan="5">1/64</td><td>1.0e+01</td><td>2.64e-10</td><td>5.06e-13</td><td>4.01e-08</td><td>2.06e-22</td><td>6</td></tr><tr><td>1.0e+00</td><td>2.11e-10</td><td>3.85e-12</td><td>4.01e-08</td><td>1.45e-22</td><td>6</td></tr><tr><td>1.0e-02</td><td>2.46e-10</td><td>2.88e-10</td><td>4.01e-08</td><td>2.00e-22</td><td>6</td></tr><tr><td>1.0e-04</td><td>1.42e-09</td><td>4.17e-08</td><td>4.03e-08</td><td>4.59e-21</td><td>5</td></tr><tr><td>1.0e-06</td><td>6.33e-09</td><td>1.72e-06</td><td>4.22e-08</td><td>9.48e-21</td><td>4</td></tr><tr><td rowspan="5">1/256</td><td>1.0e+01</td><td>3.12e-10</td><td>6.09e-13</td><td>4.20e-11</td><td>3.07e-22</td><td>6</td></tr><tr><td>1.0e+00</td><td>3.59e-10</td><td>6.68e-12</td><td>4.24e-11</td><td>4.29e-22</td><td>6</td></tr><tr><td>1.0e-02</td><td>2.14e-10</td><td>2.50e-10</td><td>4.13e-11</td><td>1.97e-22</td><td>6</td></tr><tr><td>1.0e-04</td><td>1.53e-09</td><td>4.47e-08</td><td>7.43e-11</td><td>4.68e-21</td><td>5</td></tr><tr><td>1.0e-06</td><td>4.71e-09</td><td>1.49e-06</td><td>6.58e-10</td><td>7.82e-21</td><td>4</td></tr><tr><td rowspan="10">8</td><td rowspan="5">1/64</td><td>1.0e+01</td><td>1.66e-09</td><td>6.25e-12</td><td>4.02e-08</td><td>3.00e-20</td><td>6</td></tr><tr><td>1.0e+00</td><td>1.83e-09</td><td>6.10e-11</td><td>4.02e-08</td><td>3.36e-20</td><td>6</td></tr><tr><td>1.0e-02</td><td>7.39e-10</td><td>2.23e-09</td><td>4.02e-08</td><td>4.28e-21</td><td>6</td></tr><tr><td>1.0e-04</td><td>7.98e-10</td><td>2.36e-08</td><td>4.02e-08</td><td>1.01e-21</td><td>5</td></tr><tr><td>1.0e-06</td><td>6.35e-09</td><td>2.01e-06</td><td>4.22e-08</td><td>2.23e-20</td><td>4</td></tr><tr><td rowspan="5">1/256</td><td>1.0e+01</td><td>8.14e-10</td><td>3.05e-12</td><td>4.59e-11</td><td>7.94e-21</td><td>6</td></tr><tr><td>1.0e+00</td><td>1.31e-09</td><td>4.43e-11</td><td>4.96e-11</td><td>1.78e-20</td><td>6</td></tr><tr><td>1.0e-02</td><td>3.64e-10</td><td>1.18e-09</td><td>4.26e-11</td><td>1.31e-21</td><td>6</td></tr><tr><td>1.0e-04</td><td>1.01e-09</td><td>3.25e-08</td><td>4.28e-11</td><td>1.61e-21</td><td>5</td></tr><tr><td>1.0e-06</td><td>5.09e-09</td><td>1.65e-06</td><td>6.69e-10</td><td>2.01e-20</td><td>4</td></tr><tr><td rowspan="10">16</td><td rowspan="5">1/64</td><td>1.0e+01</td><td>1.41e-09</td><td>5.91e-12</td><td>4.02e-08</td><td>4.37e-20</td><td>7</td></tr><tr><td>1.0e+00</td><td>1.05e-09</td><td>4.70e-11</td><td>4.02e-08</td><td>2.57e-20</td><td>7</td></tr><tr><td>1.0e-02</td><td>2.11e-10</td><td>6.16e-10</td><td>4.01e-08</td><td>4.69e-22</td><td>7</td></tr><tr><td>1.0e-04</td><td>4.50e-10</td><td>2.40e-08</td><td>4.02e-08</td><td>5.22e-22</td><td>5</td></tr><tr><td>1.0e-06</td><td>5.56e-09</td><td>1.55e-06</td><td>4.14e-08</td><td>1.43e-20</td><td>4</td></tr><tr><td rowspan="5">1/256</td><td>1.0e+01</td><td>1.25e-09</td><td>5.35e-12</td><td>5.00e-11</td><td>3.53e-20</td><td>7</td></tr><tr><td>1.0e+00</td><td>1.54e-09</td><td>6.75e-11</td><td>5.30e-11</td><td>5.62e-20</td><td>7</td></tr><tr><td>1.0e-02</td><td>2.78e-10</td><td>1.17e-09</td><td>4.24e-11</td><td>1.70e-21</td><td>7</td></tr><tr><td>1.0e-04</td><td>5.00e-10</td><td>2.56e-08</td><td>4.21e-11</td><td>6.82e-22</td><td>5</td></tr><tr><td>1.0e-06</td><td>4.59e-09</td><td>1.37e-06</td><td>5.92e-10</td><td>1.26e-20</td><td>4</td></tr></tbody></table>

studied a particular example of such a new domain decomposition method, namely when the Dirichlet trace jump is also added as a penalization term in the augmented Lagrangian formulation.

Many other new domain decomposition methods can be obtained this way, for example by also including in the penalization the Neumann jump, or the entire PDE constraint. Our choice led to a domain decomposition method which is robust under
<figure><img src="image_8.png" /><figcaption>Fig. 6: Convergence of the new DD algorithm with different values of α on the second example.</figcaption></figure>

mesh refinement, which lets us believe that this approach for designing new domain decomposition methods is a rather powerful one.

Based on the numerical results from Section 4, we give now an outlook on future research directions we foresee. First, a comment on the parallelizability of our new DD method is in order: in our prototype implementation, we solve (8) using a global solver encompassing all domains, hence our implementation was not parallel. Nevertheless, a well designed computation of the solution of (8) can become highly parallel. To see this, consider the simple case with two subdomains $\Omega_1$ and $\Omega_2$ having a common interface $\Gamma$. The necessary and sufficient conditions of optimality for (8) use the adjoint
functions $p_i^k$ defined by

$$ egin{cases} -\Delta p_i^k + y_i^k = y_{\text{target}}, \\ p_i^k|_{\partial\Omega} = 0, \\ \partial_{\mathbf{n}_i} p_i^k = (-1)^i (\lambda^k + \rho^k(y_1^k - y_2^k)) \text{ on } \Gamma. \end{cases} \quad (28) $$

The derivative of $\hat{\mathcal{J}}^+$ w.r.t. $g^k$ is $p_1^k|_{\Gamma} - p_2^k|_{\Gamma}$ and the derivative w.r.t. $u_i^k$ is $u_i^k - \alpha p_i^k$, $i = 1, 2$. Using a gradient descent method or a quasi-Newton method, problem (8) can then be solved in parallel as follows:

1. Starting with some initial guess $u_i^{k,l}, g^{k,l}$ at $l = 0$, we compute $y_i^{k,l}$ for $i = 1, 2$ in parallel solving (6).

2. The resulting $y_i^{k,l}$ are then used to compute $p_i^{k,l}, i = 1, 2$ in parallel by solving (28).

3. Using $y_i^{k,l}$ and $p_i^{k,l}, i = 1, 2$, we can compute $\hat{\mathcal{J}}^+$ and its derivatives.

4. We now set $l := l + 1$, update $u_i^{k,l}$ and $g^{k,l}$ using the gradient, and repeat the process by going back to step 1 until convergence.

In this gradient descent method the subdomains have to exchange interface data twice: first in order to compute the gap $y_1|_{\Gamma} - y_2|_{\Gamma}$ in Step 2, and then to compute the gap $p_1|_{\Gamma} - p_2|_{\Gamma}$ in Step 3. This thus corresponds to a typical domain decomposition solver by subdomain iteration, where both $y_i^{k,l}$ and $p_i^{k,l}$ can be computed independently on each subdomain. Note also that the gap $p_1^{k,l}|_{\Gamma} - p_2^{k,l}|_{\Gamma}$ only needs to be reduced to a threshold $\omega^k$ as prescribed in Algorithm 2.1, and thus the adjoint function is not necessarily continuous for the whole Algorithm 2.1 to converge. Thoroughly studying and testing this approach is however beyond the scope of the present manuscript and will be future work.

A naturally following second remark is that the new DD method is in fact an inner outer iteration method: the outer iteration is an update of $\rho$ or $\lambda$, and the inner iteration comes from solving (8) with a parallelizable iterative method as described above. Our numerical results from section 4 indicate that the number of outer iterations needed to converge does not depend on the number of subdomains. This is because (8) is an optimization problem defined on the entire domain, inducing global coupling through the $p_i$. Since we used quadprog in our prototype implementation to solve (8) globally, we did not measure the number of inner iterations needed to solve this optimization problem. We expect that these inner iterations increase with the number of subdomains as it is typically observed with one level DD algorithms [15]. An analysis of the convergence of these inner iterations and the addition of a coarse space is thus also a very important future research direction.

Finally, a major next step is to adapt these results to PDE-constrained optimization problems with non-linear PDEs as constraints like e.g. topology optimization [52, 54]. This approach is close to the direct multiple shooting method used for the numerical solution of optimal control problems of ODEs (see [31]), and the generalization of our new approach to non-linear equations seems straightforward. However, the analysis and the efficiency of this method applied to non-linear models is currently a wide open question.
A Surjectivity of operator M*

In this section, we will focus on proving that the derivative of the continuity con-
straint given by $c_o : (u_1, u_2, g) \in \mathcal{U} \mapsto y_1(u_1, g)|_{\Gamma} - y_2(u_2, g)|_{\Gamma} \in \mathcal{V} \cap$, where $y_i(u_i, g)$
solves (11), is surjective. This condition (called regularity condition or constraint qual-
ification) serves in order to prove that a minimizer of (10) respects the first order
conditions of optimality (or KKT conditions). We recall that we denote by $M^*$ the
derivative of $c_o$ with respect to $(u_1, u_2, g)$. We first derive an explicit expression for
$M$. To do this, we introduce the Lagrangian of (10),

$$
\begin{align*}
\mathcal{L}(y, u, g, p, \lambda) = & \sum_{i=1}^{2} \int_{\Omega_i} \nabla y_i \cdot \nabla p_i - (-1)^{i+1} \int_{\Gamma} g p_i - \int_{\Omega_i} (f_i + u_i) p_i \\
& + \frac{1}{2} \int_{\Omega_i} (y_i - y_{\text{target}})^2 + \frac{\alpha}{2} \int_{\Omega_i} u_i^2 + \int_{\Gamma} (y_1 - y_2)\lambda.
\end{align*}
$$

where we set $u := (u_1, u_2)$, $y := (y_1, y_2)$ and $p := (p_1, p_2)$. In order to derive an explicit expression of $M$, we compute the derivative of $\mathcal{L}$ and obtain

$$
\partial_{u,g} \mathcal{L}(y(u,g), u, g, p, \lambda) = \begin{pmatrix} \alpha u_1 - p_1 \\ \alpha u_2 - p_2 \\ p_1 - p_2 |_{\Gamma} \end{pmatrix} = \partial_{u,g} \hat{\mathcal{J}}(u,g) + M\lambda,
$$

where $M : \mathcal{V}_{\cap}^* \to \mathcal{U}^*$ is the adjoint of $M^*$. Using classical computations, we can also
prove that

$$
\partial_{u,g} \hat{\mathcal{J}}(u,g) = \begin{pmatrix} \alpha u_1 - \bar{p}_1 \\ \alpha u_2 - \bar{p}_2 \\ \bar{p}_2 |_{\Gamma} - \bar{p}_1 |_{\Gamma} \end{pmatrix} \quad \text{where} \quad \begin{cases} -\Delta \bar{p}_i + y_i = y_{\text{target}}, \\ \bar{p}_i |_{\partial \Omega} = 0, \\ \partial_{\mathbf{n}} \bar{p}_i = 0 \text{ on } \Gamma. \end{cases}
$$

Therefore, we get

$$
M\lambda = \partial_{u,g}\mathcal{L}(y(u,g), u, g, p, \lambda) - \partial_{u,g}\hat{\mathcal{J}}(u,g) = \begin{pmatrix} q_1 \\ q_2 \\ q_1|_{\Gamma} - q_2|_{\Gamma} \end{pmatrix},
$$

where $q_i = \bar{p}_i - p_i \in H^1(\Omega_i)$ is the weak solution of

$$
\left\{
\begin{array}{l}
-\Delta q_i = 0, \\
q_i|_{\partial \Omega} = 0, \\
\partial_{\mathbf{n}} q_i = (-1)^{i+1} \lambda \text{ on } \Gamma.
\end{array}
\right.
\qquad (29)
$$

Note that $M$ does not depend on the controls.

**Lemma 4** The operator $M : \lambda \in \mathcal{V}_{\cap}^* \mapsto M\lambda \in \mathcal{U}^*$ is injective with closed range.
*Proof Injectivity:* Let $\lambda$ be such that $M\lambda = 0$. The associated $q_i$ are then harmonic and verify the transmission condition (on $\Gamma$) $q_1|_{\Gamma} - q_2|_{\Gamma} = 0$ and $\partial_{n_1}q_1 + \partial_{n_2}q_2 = 0$. As a result, $q_i = q|_{\Omega_i}$ where $q \in H_0^1(\Omega)$ is harmonic. Therefore $q = 0$ and thus $\lambda = \partial_{n_1}q_1 = 0$.

*Closed range:* Let $M\lambda_n = (q_{1,n}, q_{2,n}, \varphi_n)^{\intercal}$ be a sequence of images such that $q_{i,n}$ converges toward some $q_i$ in $L^2(\Omega_i)$ and $\varphi_n$ converges toward some $\varphi$ in $H_{00}^{1/2}(\Gamma)$. To show that $M$ has closed range, we have to prove that there exists $\lambda \in H_{00}^{1/2}(\Gamma)^*$ such that $M\lambda = (q_1, q_2, \varphi)^{\intercal}$.

**Step 1:** Since $q_{i,n}$ is harmonic on $\Omega_i$, we have

$$ \forall \psi_i \in C_c^\infty(\Omega_i) : \int_{\Omega_i} q_{i,n} \Delta \psi_i dx = 0. $$

Using the $L^2$ convergence of $q_{i,n}$ toward $q_i$, we obtain that $q_i$ satisfies

$$ \forall \psi_i \in C_c^\infty(\Omega_i) : \int_{\Omega_i} q_i \Delta \psi_i dx = 0, $$

and Weyl's Lemma (see e.g. [55, p. 78, Theorem 18.G]) ensures that $q_i \in C^\infty(\Omega_i)$ and satisfies $\Delta q_i = 0$ pointwise in $\Omega_i$.

**Step 2:** Let us consider the spaces

$$ H_{i,\Delta} := \{ \Phi_i \in H^1(\Omega_i) \mid \Delta\Phi_i \in L^2(\Omega_i), \Phi_i|_{\partial\Omega_i \setminus \Gamma} = 0 \}. $$

For any $\varphi \in H_{00}^{1/2}(\Gamma)$, we denote by $\tilde{\varphi}$ its extension by zero to $\partial\Omega_i$. The latter is in $H^{1/2}(\partial\Omega_i)$ and the surjectivity of the trace operator gives the existence of some $\mathcal{E}\tilde{\varphi} \in H^1(\Omega_i)$ such that $\mathcal{E}\tilde{\varphi}|_{\partial\Omega_i} = \tilde{\varphi}$ and $\|\mathcal{E}\tilde{\varphi}\|_{H^1(\Omega_i)} \lesssim \|\tilde{\varphi}\|_{H_{00}^{1/2}(\partial\Omega_i)}$. From the Green's formula, for all $\Phi_i \in H_{i,\Delta}$,

$$ \langle \partial_{n_i} \Phi_i, \varphi \rangle_{H_{00}^{1/2}(\Gamma)^* \times H_{00}^{1/2}(\Gamma)} = \int_{\Omega_i} \nabla \Phi_i \cdot \nabla \mathcal{E}\tilde{\varphi} dx + \int_{\Omega_i} \Delta \Phi_i \mathcal{E}\tilde{\varphi} dx, $$

and thus, we can prove, using the Cauchy-Schwarz inequality, that

$$ \left| \langle \partial_{n_i} \Phi_i, \varphi \rangle_{H_{00}^{1/2}(\Gamma)^* \times H_{00}^{1/2}(\Gamma)} \right| \lesssim \|\varphi\|_{H_{00}^{1/2}(\Gamma)} \left( \|\Phi_i\|_{H^1(\Omega_i)} + \|\Delta\Phi_i\|_{L^2(\Omega_i)} \right). $$

As a result, taking the supremum over all $\varphi$ such that $\|\varphi\|_{H_{00}^{1/2}(\Gamma)} = 1$, we prove that the linear map

$$ \partial_{n_i} : \Phi_i \in H_{i,\Delta} \mapsto \partial_{n_i} \Phi_i \in H_{00}^{1/2}(\Gamma)^* $$

is continuous. In addition, for any harmonic $\Phi_i \in H_{i,\Delta}$, we have the bound

$$ \| \partial_{n_i} \Phi_i \|_{H_{00}^{1/2}(\Gamma)^*} \lesssim \| \Phi_i \|_{H^1(\Omega_i)} \lesssim \| \nabla \Phi_i \|_{L^2(\Omega_i)}, $$

where we used the Poincaré inequality to get the last upper bound.
**Step 3:** Both $q_{i,n}$ satisfy the weak formulation

$$ \forall \psi_i \in V_i : \int_{\Omega_i} \nabla q_{i,n} \cdot \nabla \psi_i dx = \langle \partial_{\mathbf{n}_i} q_{i,n}, \psi_i \rangle_{H_{00}^{1/2}(\Gamma)^* \times H_{00}^{1/2}(\Gamma)}, $$

where $V_i := \{ \psi_i \in H^1(\Omega_i) \mid \psi_i|_{\partial\Omega_i \setminus \Gamma} = 0 \}$. Now taking $\psi_i = q_{i,n}$, one gets

$$ \begin{aligned} \sum_i \|\nabla q_{i,n}\|_{L^2(\Omega_i)}^2 &= \langle \partial_{\mathbf{n}_i} q_{1,n}, q_{1,n} \rangle_{H_{00}^{1/2}(\Gamma)^* \times H_{00}^{1/2}(\Gamma)} + \langle \partial_{\mathbf{n}_2} q_{2,n}, q_{2,n} \rangle_{H_{00}^{1/2}(\Gamma)^* \times H_{00}^{1/2}(\Gamma)} \\ &= \langle \partial_{\mathbf{n}_1} q_{1,n}, q_{1,n} - q_{2,n} \rangle_{H_{00}^{1/2}(\Gamma)^* \times H_{00}^{1/2}(\Gamma)} \\ &= \langle \partial_{\mathbf{n}_1} q_{1,n}, \varphi_n \rangle_{H_{00}^{1/2}(\Gamma)^* \times H_{00}^{1/2}(\Gamma)}, \end{aligned} $$

where we used that $\partial_{\mathbf{n}_1} q_{1,n}|_{\Gamma} + \partial_{\mathbf{n}_2} q_{2,n}|_{\Gamma} = 0$ and that, due to the definition of $M$, we have $\varphi_n = q_{1,n} - q_{2,n}$. Using now **Step 2**, we obtain that

$$ \sum_i \|\nabla q_{i,n}\|_{L^2(\Omega_i)}^2 \le \|\partial_{\mathbf{n}_i} q_{1,n}\|_{H_{00}^{1/2}(\Gamma)^*} \|\varphi_n\|_{H_{00}^{1/2}(\Gamma)} \lesssim \|\nabla q_{1,n}\|_{L^2(\Omega_1)} \|\varphi_n\|_{H_{00}^{1/2}(\Gamma)}. $$

Using finally Young's inequality, we get

$$ \sum_i \|\nabla q_{i,n}\|_{L^2(\Omega_i)}^2 \lesssim \|\varphi_n\|_{H_{00}^{1/2}(\Gamma)}^2. $$

Since $\varphi_n$ is assumed to converge toward $\varphi$, it is bounded and thus the sequences $(q_{i,n})_n$ are also bounded (uniformly with respect to $n$) in $H^1(\Omega_i)$. We can then extract subsequences that converge weakly in $H^1(\Omega_i)$ toward $q_i$. The trace operator being compact, we obtain that $q_i|_{\partial\Omega_i \setminus \Gamma} = 0$ and $\varphi = q_1|_{\Gamma} - q_2|_{\Gamma}$. Since $q_i \in V_i$ is harmonic, it also belongs to $H_{i,\Delta}$ and its normal derivative can be defined as an element of $H_{00}^{1/2}(\Gamma)^*$.

**Step 4:** We now identify the limit of the sequence $\lambda_n$. Let $\eta \in H_{00}^{1/2}(\Gamma)$, since its extension by zero (denoted by $\tilde{\eta}$) over $\partial\Omega_i$ is in $H^{1/2}(\partial\Omega_i)$, we have some $\mathcal{E}_i \tilde{\eta} \in V_i$ such that $\mathcal{E}_i \tilde{\eta}|_{\partial\Omega_i} = \tilde{\eta}$ and $\|\mathcal{E}_i \tilde{\eta}\|_{H^1(\Omega_i)} \lesssim \|\tilde{\eta}\|_{H_{00}^{1/2}(\Gamma)}$. Using that $\lambda_n = (-1)^{i+1} \partial_{\mathbf{n}_i} q_{i,n}$ and the weak formulation satisfied by $q_{i,n}$, we have

$$ \begin{aligned} \langle \lambda_n, \eta \rangle_{H_{00}^{1/2}(\Gamma)^* \times H_{00}^{1/2}(\Gamma)} &= \frac{1}{2} \langle \partial_{\mathbf{n}_1} q_{1,n}, \mathcal{E}_1 \tilde{\eta} |_{\Gamma} \rangle_{H_{00}^{1/2}(\Gamma)^* \times H_{00}^{1/2}(\Gamma)} \\ &\quad - \frac{1}{2} \langle \partial_{\mathbf{n}_2} q_{2,n}, \mathcal{E}_2 \tilde{\eta} |_{\Gamma} \rangle_{H_{00}^{1/2}(\Gamma)^* \times H_{00}^{1/2}(\Gamma)} \\ &= \frac{1}{2} \int_{\Omega_1} \nabla q_{1,n} \cdot \nabla \mathcal{E}_1 \tilde{\eta} \,dx - \frac{1}{2} \int_{\Omega_2} \nabla q_{2,n} \cdot \nabla \mathcal{E}_2 \tilde{\eta} \,dx. \end{aligned} $$

Now passing to the limit (after extracting a subsequence), we obtain that

$$ \lim_{n \to +\infty} \langle \lambda_n, \eta \rangle_{H_{00}^{1/2}(\Gamma)^* \times H_{00}^{1/2}(\Gamma)} = \frac{1}{2} \int_{\Omega_1} \nabla q_1 \cdot \nabla \mathcal{E}_1 \tilde{\eta} \,dx - \frac{1}{2} \int_{\Omega_2} \nabla q_2 \cdot \nabla \mathcal{E}_2 \tilde{\eta} \,dx $$
$$ = \left\langle \frac{1}{2} (\partial_{\mathbf{n}_1} q_1 - \partial_{\mathbf{n}_2} q_2), \eta \right\rangle_{H_{00}^{1/2}(\Gamma)^* \times H_{00}^{1/2}(\Gamma)}. $$

We have then proved that $(\lambda_n)_n \subset H_{00}^{1/2}(\Gamma)^*$ has a subsequence converging toward $\lambda = \frac{1}{2}(\partial_{\mathbf{n}_1} q_1 - \partial_{\mathbf{n}_2} q_2)$. We emphasize that each $q_i$ is unique since it is defined as the $L^2$-limit of $q_{i,n}$ and then the limit of the subsequence of $(\lambda_n)$ is also unique. Urysohn's subsequence principle finally proves that the whole sequence $\lambda_n$ converges toward $\lambda$. Since $M\lambda = (q_1, q_2, q_1|\Gamma - q_2|\Gamma)^t = (q_1, q_2, \varphi)^t$ this proves that the range of $M$ is closed.

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

Applying now [10, Théorème II.20], we obtain theorem 2.

## B Surjectivity of operator $M^*$ when discretized with finite element

Using a standard Galerkin method to discretize the continuous optimal control problem (10) following the notations in [19], we can go left the black arrow in the middle of Figure 1 and obtain the discrete optimal control problem

$$ \begin{array}{l} \min \sum_{i=1}^{2} \frac{1}{2}(\mathbf{y}_i - \mathbf{y}_{\text{target}_i})^{\top} N_i (\mathbf{y}_i - \mathbf{y}_{\text{target}_i}) + \frac{\alpha}{2} \mathbf{u}_i^{\top} N_i \mathbf{u}_i \\ \text{s.t.} \quad \begin{cases} K_i \mathbf{y}_i = \mathbf{F}_i + N_i \mathbf{u}_i + (-1)^{i+1} B_i^{\top} \mathbf{g}, & i = 1, 2, \\ B_1 \mathbf{y}_1 - B_2 \mathbf{y}_2 = 0, \end{cases} \end{array} \quad (30) $$

where $\mathbf{y}_i \in \mathbb{R}^{n_i^s+n^I}$ is the vector of degrees of freedom of the finite element approximation of $y_i$, $n_i^s$ is the number of interior nodes in $\Omega_i$ and $n^I$ the number of interface nodes on $\Gamma$, and $B_i$ is of the form $B_i = [0_i, I_i]$, $i = 1, 2$, where $0_i$ is an $n^I \times n_i^s$ zero matrix, and $I_i$ is the $n^I \times n_i^I$ identity matrix. $K_i$ is the stiffness matrix and $N_i$ the mass matrix, that are both invertible.

Solving this discrete optimal control problem (30) with the same new DD algorithm based on the augmented Lagrangian approach leads to the decoupled augmented discrete Lagrangian

$$ \begin{aligned} \min \sum_{i=1}^{2} & \frac{1}{2}(\mathbf{y}_i(\mathbf{u}_i, \mathbf{g}) - \mathbf{y}_{\text{target}_i})^{\top} N_i(\mathbf{y}_i(\mathbf{u}_i, \mathbf{g}) - \mathbf{y}_{\text{target}_i}) + \frac{\alpha}{2} \mathbf{u}_i^{\top} N_i \mathbf{u}_i \\ & + (\lambda^k)^{\top}(B_1 \mathbf{y}_1(\mathbf{u}_1, \mathbf{g}) - B_2 \mathbf{y}_2(\mathbf{u}_2, \mathbf{g})) + \frac{\rho^k}{2} \|B_1 \mathbf{y}_1(\mathbf{u}_1, \mathbf{g}) - B_2 \mathbf{y}_2(\mathbf{u}_2, \mathbf{g})\|_2^2, \end{aligned} \quad (31) $$

where $\mathbf{y}_i(\mathbf{u}_i, \mathbf{g})$ satisfy

$$ K_i \mathbf{y}_i = \mathbf{F}_i + N_i \mathbf{u}_i + (-1)^{i+1} B_i^{\top} \mathbf{g}, \quad i = 1, 2. $$

As in the continuous framework, the convergence of our new DD method given by Algorithm 2.1 depends on the surjectivity of the derivative of $M_h^* : (\mathbf{u}_1, \mathbf{u}_2, \mathbf{g}) \mapsto$
$B_1y_1(u_1,g) - B_2y_2(u_2,g)$ where $y_i(u_i,g)$ verify $K_iy_i = F_i + N_iu_i + (-1)^{i+1}B_i^Tg$, $i = 1,2$ (see (20) and A for this property at the continuous level).

**Proposition 3** The derivative of $M_h^*$ w.r.t. $(u_1, u_2, g)$ is onto.

*Proof* Note first that

$$ B_1y_1(u_1,g) - B_2y_2(u_2,g) = (B_1K_1^{-1}N_1 - B_2K_2^{-1}N_2 B_1K_1^{-1}B_1^T - B_2K_2^{-1}B_2^T) \begin{pmatrix} u_1 \\ u_2 \\ g \end{pmatrix} \\ + B_1K_1^{-1}F_1 - B_2K_2^{-1}F_2. $$

Therefore, the derivative of $M_h^*$ is simply given by

$$ \partial M_h^* = (B_1K_1^{-1}N_1 - B_2K_2^{-1}N_2 B_1K_1^{-1}B_1^T - B_2K_2^{-1}B_2^T). $$

In order to prove that $\partial M_h^*$ is surjective, we will prove that $(\partial M_h)$ is injective. Let $\lambda \in \mathbb{R}^{n^I}$. Suppose $(\partial M_h)\lambda = 0$. This implies that

$$ \begin{gathered} K_1^{-T} B_1^T N_1 \lambda = 0 \\ -K_2^{-T} B_2^T N_2 \lambda = 0 \iff B_1^T N_1 \lambda = 0 \\ B_2^T N_2 \lambda = 0 \iff \lambda = 0. \end{gathered} \\ (B_1 K_1^{-T} B_1^T - B_2 K_2^{-T} B_2^T) \lambda = 0 \end{gathered} $$

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

The proof of the convergence of our new discrete DD algorithm to the solutions of (30) can then be done as in the Section 3, or can be found in [6]. Note also that in this discrete context, the extension to a striped decomposition with more subdomains becomes easier, since only the injectivity of $M$ needs to be proved.

# Declarations

## Funding

A. Vieira and P.H. Cocquet have been supported by the "Agence Nationale de la Recherche" (ANR), Project O-TO-TT-FU number ANR-19-CE40-0011 (2019-2023). P.H. Cocquet is supported by the ANR project DGMod (ANR-23-CE56-0007, 2024-2027). M.J. Gander is supported by the Swiss National Science Foundation. The authors also thank the CIRM for their hospitality and giving us the opportunity to work together during a "Recherche en résidence".

## Ethics declaration

Not applicable.
References

[1] Attouch, H., Soueycatt, M.: Augmented Lagrangian and proximal alternating direction methods of multipliers in Hilbert spaces. Applications to games, PDE's and control. Pacific Journal of Optimization 5(1), 17–37 (2008)

[2] Avery, P., Farhat, C.: The FETI family of domain decomposition methods for inequality-constrained quadratic programming: Application to contact problems with conforming and nonconforming interfaces. Computer Methods in Applied Mechanics and Engineering 198(21-26), 1673–1683 (2009)

[3] Bergounioux, M.: Augmented Lagrangian method for distributed optimal con- trol problems with state constraints. Journal of Optimization Theory and Applications 78(3), 493–521 (1993)

[4] Bergounioux, M., Haddou, M.: A SQP-augmented Lagrangian method for optimal control of semilinear elliptic variational inequalities. In: Control and Estimation of Distributed Parameter Systems, pp. 57–72. Birkhäuser Basel, Basel (2003)

[5] Bergounioux, M., Kunisch, K.: Augmented Lagrangian techniques for ellip- tic state constrained optimal control problems. SIAM Journal on Control and Optimization 35(5), 1524–1543 (1997)

[6] Bertsekas, D.P.: Constrained Optimization and Lagrange Multiplier Methods. Academic press, New York (2014)

[7] Bonnans, J.F., Laurent-Varin, J.: Computation of order conditions for sym- plectic partitioned Runge-Kutta schemes with application to optimal control. Numerische Mathematik 103, 1–10 (2006)

[8] Bounaim, A.: Méthodes de décomposition de domaine: Application à la résolution de problèmes de contrôle optimal. PhD thesis, Université Joseph-Fourier-Grenoble I (1999)

[9] Bourgat, J.-F., Glowinski, R., Le Tallec, P., Vidrascu, M.: Variational formula- tion and algorithm for trace operator in domain decomposition calculations. In: Chan, T., Glowinski, R., Périaux, J., Widlund, O. (eds.) Domain Decomposition Methods, pp. 3–16. SIAM, Philadelphia, PA (1989)

[10] Brézis, H.: Analyse Fonctionnelle: Théorie et Applications. Collection Mathématiques appliquées pour la maîtrise. Masson, Paris (1983)

[11] Casas, E.: Control of an elliptic problem with pointwise state constraints. SIAM Journal on Control and Optimization 24(6), 1309–1318 (1986)

[12] Chaudet-Dumas, B., Deteix, J.: Shape derivatives for an augmented Lagrangian formulation of elastic contact problems. ESAIM: Control, Optimisation and
[13] Chen, Z., Zou, J.: An augmented Lagrangian method for identifying discontinuous parameters in elliptic systems. SIAM Journal on Control and Optimization **37**(3), 892–910 (1999)

[14] Chyba, M., Hairer, E., Vilmart, G.: The role of symplectic integrators in optimal control. Optimal control applications and methods **30**(4), 367–382 (2009)

[15] Ciaramella, G., Gander, M.J.: Iterative Methods and Preconditioners for Systems of Linear Equations. Society for Industrial and Applied Mathematics, Philadelphia, PA (2022)

[16] Ciaramella, G., Halpern, L., Mechelli, L.: Convergence analysis and optimization of a Robin Schwarz waveform relaxation method for periodic parabolic optimal control problems. Technical report, MOX Lab, Politecnico di Milano (2022)

[17] Ciaramella, G., Kwok, F., Müller, G.: Nonlinear optimized Schwarz preconditioner for elliptic optimal control problems. In: Domain Decomposition Methods in Science and Engineering XXVI. Springer, Hong Kong, China (2020)

[18] Discacciati, M., Gervasio, P., Quarteroni, A.: The interface control domain decomposition (ICDD) method for elliptic problems. SIAM Journal on Control and Optimization **51**(5), 3434–3458 (2013)

[19] Farhat, C., Roux, F.-X.: A method of finite element tearing and interconnecting and its parallel solution algorithm. International Journal for Numerical Methods in Engineering **32**(6), 1205–1227 (1991)

[20] Gander, M.J.: Optimized Schwarz methods. SIAM Journal on Numerical Analysis **44**(2), 699–731 (2006)

[21] Gander, M.J., Kwok, F.: Schwarz methods for the time-parallel solution of parabolic control problems. In: Domain Decomposition Methods in Science and Engineering XXII, pp. 207–216 (2016). Springer

[22] Gander, M.J., Lu, L.-D.: Dirichlet-Neumann and Neumann-Neumann methods for elliptic control problems. In: Domain Decomposition Methods in Science and Engineering XXVII, pp. 207–216 (2023). Springer

[23] Gander, M.J., Lu, L.-D.: Dirichlet-Neumann and Neumann-Neumann methods for parabolic control problems. in preparation (2023)

[24] Gander, M.J., Zhang, H.: Schwarz methods by domain truncation. Acta Numerica **31**, 1–134 (2022)

[25] Gander, M.J., Kwok, F., Mandal, B.C.: Convergence of substructuring methods for elliptic optimal control problems. In: Domain Decomposition Methods in
Science and Engineering XXIV 24, pp. 291–300 (2018). Springer

[26] Gander, M.J., Kwok, F., Wanner, G.: Constrained optimization: From Lagrangian mechanics to optimal control and PDE constraints. Optimization with PDE Constraints: ESF Networking Program’OPTPDE’, 151–202 (2014)

[27] Gong, W., Kwok, F., Tan, Z.: Convergence analysis of the Schwarz alternating method for unconstrained elliptic optimal control problems (2022)

[28] Gunzburger, M., Lee, J.: A domain decomposition method for optimization problems for partial differential equations. Computers & Mathematics with Applications 40(2-3), 177–192 (2000)

[29] Hager, W.W.: Runge-Kutta methods in optimal control and the transformed adjoint system. Numerische Mathematik 87, 247–282 (2000)

[30] Hairer, E., Lubich, C., Wanner, G.: Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations. Springer Series in Computational Mathematics. Springer, Berlin (2006)

[31] Hesse, H.K.: Multiple shooting and mesh adaptation for pde constrained optimization problems. PhD thesis, Ruprecht-Karls-Universität - Heidelberg (2008)

[32] Hestenes, M.R.: Multiplier and gradient methods. Journal of Optimization Theory and Applications 4(5), 303–320 (1969)

[33] Hintermüller, M., Kunisch, K.: Feasible and noninterior path-following in constrained minimization with low multiplier regularity. SIAM Journal on Control and Optimization 45(4), 1198–1221 (2006)

[34] Hinze, M., Pinnau, R., Ulbrich, M., Ulbrich, S.: Optimization with PDE Constraints vol. 23. Springer, Dordrecht (2008)

[35] Hou, L.S., Lee, J.: A Robin-Robin non-overlapping domain decomposition method for an elliptic boundary control problem. International Journal of Numerical Analysis & Modeling 8(3) (2011)

[36] Ito, K., Kunisch, K.: The augmented Lagrangian method for equality and inequality constraints in Hilbert spaces. Mathematical programming 46(1), 341–360 (1990)

[37] Ito, K., Kunisch, K.: The augmented Lagrangian method for parameter estimation in elliptic systems. SIAM Journal on Control and Optimization 28(1), 113–136 (1990)

[38] Ito, K., Kunisch, K.: Augmented Lagrangian–SQP methods for nonlinear optimal control problems of tracking type. SIAM Journal on Control and Optimization 34(3), 874–891 (1996)
[39] Iusem, A., Gárciga Otero, R.: Inexact version of proximal point and augmented Lagrangian algorithms in Banach spaces. Numerical Functional Analysis and Optimization 22(5-6), 609-640 (2001)

[40] Kanzow, C., Steck, D., Wachsmuth, D.: An augmented Lagrangian method for optimization problems in Banach spaces. SIAM Journal on Control and Optimization 56(1), 272-291 (2018)

[41] Khanh, P.Q., Nuong, T.H., Théra, M.: On duality in nonconvex vector optimization in Banach spaces using augmented Lagrangians. Positivity 3(1), 49-64 (1999)

[42] Kunisch, K., Tai, X.-C.: Nonoverlapping domain decomposition methods for inverse problems. In: Proceedings of 9th International Conference on Domain Decomposition Methods, Editors P. Bjorstard, M. Espedal and D. Keyes, John Wiley and Sons (1997)

[43] Lagnese, J.E., Leugering, G.: Domain in Decomposition Methods in Optimal Control of Partial Differential Equations vol. 148. Birkhäuser, Basel (2004)

[44] Leugering, G.: Domain decomposition of optimal control problems for dynamic networks of elastic strings. Computational Optimization and Applications 16, 5-27 (2000)

[45] Leugering, G.: Nonoverlapping domain decomposition for instantaneous optimal control of friction dominated flow in a gas-network. Technical report (2020). preprint

[46] Lions, J.L., Magenes, E.: Non-homogeneous Boundary Value Problems and Applications vol. 1. Springer, Heidelberg (1972)

[47] Magoulès, F., Roux, F.-X.: Lagrangian formulation of domain decomposition methods: A unified theory. Applied Mathematical Modelling 30(7), 593-615 (2006). Parallel and Vector Processing in Science and Engineering

[48] Maruhn, J.H.: An augmented Lagrangian algorithm for optimization with equality constraints in Hilbert spaces. Master's thesis, Virginia Tech (2001)

[49] Picard, É.: Sur l'application des méthodes d'approximations successives à l'étude de certaines équations différentielles ordinaires. Journal de mathématiques pures et appliquées 9, 217-271 (1893)

[50] Powell, M.J.: A method for nonlinear constraints in minimization problems. Optimization, 283-298 (1969)

[51] Rao, A.V.: A survey of numerical methods for optimal control. Advances in the Astronautical Sciences 135(1), 497-528 (2009)
[52] Vieira, A., Cocquet, P.-H.: The boussinesq system with nonsmooth boundary conditions: Existence, relaxation, and topology optimization. SIAM Journal on Control and Optimization 60(6), 3457–3484 (2022)

[53] Vieira, A., Cocquet, P.-H.: Optimized schwarz method for coupled direct-adjoint problems applied to parameter identification in advection-diffusion equation. In: Domain Decomposition Methods in Science and Engineering XXVII, pp. 417–425. Springer, Cham (2024)

[54] Vieira, A., Bastide, A., Cocquet, P.-H.: Topology optimization for steady-state anisothermal flow targeting solids with piecewise constant thermal diffusivity. Applied Mathematics & Optimization 85(3), 1–32 (2022)

[55] Zeidler, E.: Nonlinear Functional Analysis and Its Applications: II/B: Nonlinear Monotone Operators. Springer, New York (2013)