Non-Overlapping Domain Decomposition for 1D Optimal Control Problems Governed by Time-Fractional Diffusion Equations on Coupled Domains: Optimality System and Virtual Controls

: We consider a non-overlapping domain decomposition method for optimal control problems of the tracking type governed by time-fractional diffusion equations in one space dimension, where the fractional time derivative is considered in the Caputo sense. We concentrate on a transmission problem defined on two adjacent intervals, where at the interface we introduce an iterative non-overlapping domain decomposition in the spirit of P.L. Lions for the corresponding first-order optimality system, such that the optimality system corresponding to the optimal control problem on the entire domain is iteratively decomposed into two systems on the respective sub-domains; this approach can be framed as first optimize, then decompose . We show that the iteration involving the states and adjoint states converges in the appropriate spaces. Moreover, we show that the decomposed systems on the sub-domain can in turn be interpreted as optimality systems of so-called virtual control problems on the sub-domains. Using this property, we are able to solve the original optimal control problem by an iterative solution of optimal control problems on the sub-domains. This approach can be framed as first decompose, then optimize . We provide a mathematical analysis of the problems as well as a numerical finite difference discretization using the L1-method with respect to the Caputo derivative, along with two examples in order to verify the method.


Introduction
Fractional diffusion equations have been discussed in a many articles in the context of anomalous diffusion (see, e.g., [1,2]), viscoelastic fluid models for weakly singular kernels (see, e.g., [3,4]), the flow of nanofluids with nanoparticles in a vertical channel [5], subdiffusion models in spiny neuronal dendrites [6,7], and many other important applications.In essence, hereditary processes with weakly singular memory are prone to being modelled in this way.Not surprisingly, there are many results concerning the well-posedness of such systems (see, e.g., [8,9] and references therein).Clearly, fractional optimal control problems genuinely arise in applications of time-fractional diffusion equations.Such fractional optimal control problems on standard domains were initiated by O.P. Agrawal in [10,11], where he formulated the necessary optimality condition based on the Riemann-Liouville fractional derivative.Fractional optimal control problems for time-fractional diffusion equations in the Riemann-Liouville sense were discussed by Mophou et al. in [12,13].There, the authors considered the well-posedness of the underlying time-fractional diffu-sion equaiton and derived the optimality system using the Lagrange method.Space-time fractional optimal control problems have been discussed in [14], among others.
A more recent development concerns applications on networked systems, where the domain is a union of simple domains, although one than can potentially have multiple connections.Such domains can be framed using metric graphs.To date, the literature has almost exclusively considered locally one-dimensional graphs.Applications relevant for the material of this article concern the flow of gas in pipe networks (see the website (accessed on 10 January 2024) https://www.trr154.fau.de/trr-154-en/(accessed on 10 January 2024) and the the large number of publications therein) or water networks (see, e.g., [15]).We refer readers to [16,17] and the works by Mophou et al. [18,19] for more such problems involving fractional derivatives.
Clearly, the complexity of such graphs demands domain decomposition methods.The very nature of metric graphs naturally suggests non-overlapping domain decomposition techniques at the multiple nodes.In the last decades, there has been a significant amount of research on domain decomposition of optimal control problems for classical (integer-order) partial differential equations on Euclidean domain.Such decomposition has been discussed for both the time and spatial domains as well as on the continuous level and for discretized problems; see the textbooks [20,21] for general descriptions and the website of the domain decomposition organization http://www.ddm.org/(accessed on 10 January 2024).Spatial decompositions of parabolic problems have been discussed in [22][23][24], while in [25][26][27][28][29][30][31] the authors focused on time-domain decompositions.
However, the number of articles shrinks drastically when it comes to such methods for the underlying problems on metric graphs; see [32].Domain decomposition for a fractional-time parabolic transmission problem based on the overlapping Schwarz method was considered in [33].To the best of our knowledge, this article is the first work on non-overlapping domain decompositions associated with transmission problems for timefractional diffusion equations, let alone corresponding optimal control problems for such systems.The results obtained in the following will be extended to metric graphs in a forthcoming publication.
To fix the ideas, we focus on two intervals, namely, I i := (0, ℓ i ), i = 1, 2, connected at x = 0 such that they stretch from the common point x = 0 to the ends at x = ℓ i , i = 1, 2; hence, Ω = ∪ 2 i=1 I i .By admitting the above reduction and considering x ∈ (0, ℓ i ) as the coordinate, we can formulate the following distributed optimal control problem for the time-fractional diffusion equation.To this end, we introduce the tracking-type cost function Then, the fractional optimal control problem (FOCP) is to minimize J(y, u) with respect to (y, u) while satisfying the following system of equations.
In the above problem, C D α 0,t denotes the Caputo fractional derivative of order α, 0 < α < 1 with respect to t, (y 0 i ) 2 i=1 = y 0 ∈ D(−L) (see Section 2 for definitions) is the initial state, and ( denotes the distributed controls.The formulation is deliberately written in this way in order to compare it with the problem formulation on the star graph, as in [17], and more importantly, to prepare for the domain decomposition procedure.We note that it would in fact be possible to handle different or even space-dependent coefficients in this framework as well.Equation (2) 2,3 represent the transmission conditions at the interface.

Remark 1.
In order to explain the connection between the transmission problem (2) and problems on, say, a star graph consisting of k edges I i := (0, ℓ i ), i = 1, . . ., k coupled at the point x = 0, we remark that the transmission conditions (2) 2,3 are extended to See, [17].Alluding to electrical networks, the second condition is often called the Kirchhoff condition in the context of problems on metric graphs.
The goal of this article is to introduce a non-overlapping domain decomposition of the optimal control problem (1) in the spirit of [32], which goes back to the ideas of P.L. Lions [34].We introduce and analyse two methods, which we will then show to be equivalent.The first one can be framed as first optimize, then decompose, while the second can be regarded as first decompose, then optimize.In the first case it is understood that the optimization problem is replaced by the necessary (and here, in fact, also sufficient) optimality conditions.These optimality conditions have been derived in [17] (Theorem 4.2) and are given by where we denote by C D α t,T the backward (right-sided) Caputo fractional derivative and postpone the mathematical definitions to the next section.Solving the coupled forward and backward equation is sometimes called the all-at-once-approach.The non-overlapping domain decomposition we are going to analyse in this article is described as follows.
Assume at an iteration level n that we are given data λ n ij (t), ρ n ij (t), and y n i (0, t), p n i (0, t), t ∈ (0, T).On the interval I i , we solve the problem where σ, µ ≥ 0, σ + µ > 0 are process parameters and the index j refers to the mutual other interval, i.e., i = 1, j = 2 and i = 2, j = 1.The history of the coupling is encoded in λ n ij , ρ n ij and the corresponding update at each iteration step is according to The iteration is initialized by choosing λ 0 ij , ρ 0 ij and y 0 i , p 0 i , e.g., as zero.This procedure aims to decompose the optimality system, and as such can be seen as an example of the first optimize, then decompose principle.A crucial property of the proposed procedure is that the system (4) and ( 5) is in fact equivalent to an iteration of optimal control problems on the individual edges.In these local optimal control problems, a crucial so-called virtual control h plays a major role.We formulate this method as follows.Suppose that we are given the iterates y n i , p n i , i = 1, 2, and consequently the updates in (5); then, we are looking for y n+1 i and the actual control u n+1 i as well as the virtual control(s) h ij minimizing the augmented cost function Thus, the optimal control problem on the individual edge subject to The motivation for this approach originates in the idea of using an extra control, the virtual control, in order to achieve the matching of states in the limit.The sequence of optimal control problems (7), (8) can be seen as a first decompose, then optimize approach.The interesting feature of this interpretation is that, after appropriate discretization, standard software (e.g., CASADI-OPT (accessed on 10 January 2024) https://web.casadi.org/using IPOPT https://coin-or.github.io/Ipopt/or SNOPT https://ccom.ucsd.edu/~optimizers/solvers/snopt) can be used to obtain the solution.We will demonstrate in the sequelae that as n → ∞, the iterates (y n i , u n i ) obtained from the above algorithm converge to the exact solution (y, u) of the considered optimal control problem (1) and (2).Finally, we provide some numerical evidence.We note, however, that the emphasis of this article is on the continuous level, i.e., the partial differential equation (PDE), rather than on the discrete level.In the language of PDE-constrained optimization, this follows the first optimize, then discretize approach.For α = 1, the proposed domain decomposition algorithms have been surveyed in [35].
Remark 2. We remark that the iteration update for λ n+1 ij , µ n+1 ij can be modified by a convex combination of the current update (with factor (1 − ϵ)) and the previous one (with factor ϵ).Moreover, we can use the current update for domain #1 in the update for domain #2.This results in a Gauß-Seidel-type procedure, while the original method is reminiscent of a Jacob-type iteration.

Preliminaries
In order to prepare for the main body of this article, we first recall some basic definitions for fractional derivatives, relevant properties, and space settings.For easier reference, we refer the reader to [17], even though the notations are classic.Definition 1.For a continuous function f on [a, b], the left and right fractional integrals of order α > 0 are respectively defined by where Γ(.) is the Euler gamma function.
Definition 2. For a function f with absolutely continuous derivatives up to order n − 1, n − 1 ≤ α < n, n ∈ N in (a, b), the left and right Riemann-Liouville fractional derivatives of order α > 0 are respectively defined by Definition 3.For functions as in Definition 2, the left and right Caputo fractional derivatives of order α > 0 are respectively defined by then the left and right Caputo fractional derivatives of order α > 0 are respectively provided by where n is as defined in Definition 2.
Next, we recall the following spaces Ω: which are Hilbert spaces when equipped with the inner products [36] ⟨y, where L 2 (0, ℓ i ) and H m (0, ℓ i ) are the standard spaces.We introduce the following operator L on the Hilbert space L 2 (Ω): .
The following lemma relates to an integration-by-parts formula for the fractional diffusion equation with Caputo derivative, and is by now well known, cf.[17], among others.
Using the relation between the Riemann-Liouville and Caputo derivatives and integration by parts (Lemma 1), we obtain the following lemma.

Existence and Uniqueness of Solution to TFDEs
For the well-posedness, we refer to [17] (Section 3) for detailed proofs.Here, we first consider the following TFDE on Ω: with y 0 ∈ D(−L) and g ∈ L 2 (Ω × (0, T)).
The existence results for TFDE (9a)-(9e) can be related to the Mittag-Leffler function.
We can convert (10) into a forward-running system by applying a time shift.In view of Remark 4, we then have the following result.Proposition 1.Let 0 < α < 1; then, TFDE (10) admits a unique solution y ∈ L 2 ((0, T); D(−L)).Furthermore, there exists a positive constant C 1 such that

Existence of Optimal Solution and the Optimality System
In this section, we briefly summarize the results on well-posedness for the optimal control problem (1).We first recall the cost function where z d ∈ L 2 (Ω × (0, T)) and κ, ν > 0. Now, defining u ∈ U := L 2 (Ω × (0, T)), the FOCP ( 1) and ( 2) can be written as follows: A solution to (14) is denoted as the optimal solution ( ȳ, ū) and the corresponding control is denoted as the optimal control.For the existence, we refer to [17].
For the corresponding optimality system of the considered optimal control problem (1)-( 2), we refer to the following theorem [17].
Theorem 2. Let ( ȳ, ū) be an optimal solution for (1) and (2), i.e., ( ȳ, ū) satisfies (14); then, there exists a unique p ∈ L 2 ((0, T); D(−L)) such that ( ȳ, ū, p) satisfies the following optimality system: and Remark 5. We notice that control constraints can be handled in this context with little extra work.Therefore, if we introduce the set of constraints U d := {u : , where both the functions u l = (u l i ) 1≤i≤k and u r = (u r i ) 1≤i≤k belong to L 2 (Ω × (0, T)), then instead of equality in (17), we obtain the following variational inequality using standard variational theory:

Non-Overlapping Domain Decomposition
We now consider the proposed domain decomposition procedure (4), (5).First of all, we go a step back from (5) and set Then, the update (5) easily follows and vice versa.We use (5) in the algorithm, as it avoids explicit computation of the normal derivatives of y n i and p n i .We use (19) to show consistency with the exact solution y, p of (3).To this end, let us assume convergence.Then, we may delete the iteration index n in (4), (5).We explicitly write Adding the first two and second two equations yields a system of two linear equations in the differences y 1 (0, t) − y 2 (0, t) and p 1 (0, t) − p 2 (0, t), which has the unique zero solution iff σ 2 + µ 2 ̸ = 0.When the continuity of traces is established, the Kirchhoff condition immediately follows.
From the linearity of the system, the errors ỹn i (x, t) := y n i (x, t) − y i (x, t) and pn i (x, t) := p n i (x, t) − p i (x, t) satisfy the same system, except with zero initial conditions and with z d deleted.Thus, we are looking at We notice that the updates now read We are going to show that ỹn i , pn i tend to zero as n → ∞.In the subsequent analysis, we omit the tilde for the sake of simplicity.We introduce the space X := (L 2 ((0, T)) 4 and We define the fixed point mapping T : X → X such that A simple calculation shows that where F n is provided below by (31).In order to show that T is non-expansive, we need to show that F n is positive.In fact, we will show more.
Lemma 3. We have the following identities: Proof.For the first two identities, we multiply the forward and the adjoint equation by y i , p i , respectively, and integrate by parts.The second two identities (3) follow from multiplying the forward and adjoint equation by p i , y i , respectively.
With the identities in Lemma 3 at hand, we are now able to evaluate F n .
Lemma 4. We have Proof.Indeed, according to Lemma 3, we have Going back to the definition of F in (26), we arrive at the desired result.
In order to estimate the fractional derivatives in F , we prove the following results.

Lemma 5.
T 0 T 0 Proof.We first note (see [37] Lemma 1) that for any absolutely continuous v we have Therefore, in view of (34) and using Fubini's theorem (two times) and the Leibniz integral rule, we obtain where the last equality in the above expression follows from the integration by parts.For (33), we can use a similar inequality for the backward operator and follow the proof as above to obtain the desired inequality.
We are now in position to prove our convergence result.
Theorem 3. Let the optimization penalties κ ≥ 0, µ > 0 and the initial condition y 0 ∈ L 2 (Ω) be given, and assume that we have initial transmission data λ 0 ij , ρ 0 ij ∈ L 2 (0, T), i, j = 1, 2. We consider the errors ỹi = y n i − y i and pi = p n i − p i , where y i , p i solve the original optimality system (3) and the iterates y n+1 i , p n+1 i solve (4).Then, the iteration (4) converges under the following conditions.
Proof.Recalling ( 26) and (31), we set X n+1 = T X n .We then iterate the equality Now, whether F n is positive (definite) for each n depends on the constellation of the parameters σ, µ, κ, ν, in which case we can conclude that F n → 0 as n → ∞.

Case (i.)
In this case, Therefore, we have the stated convergence.
Case (ii.)According to Lemma 4, in particular Formula (31), along with Lemma 5, we know that the terms involving the fractional derivatives are at least non-negative.Therefore, the first space-time integrals in (31) provide the convergence in L 2 (0, T, H 1 (0, ℓ i )).
The problem then reduces to estimating the mixed zero-order terms T 0 y i p i dxdt.We use a generalized Cauchy-Schwartz inequality ab ≤ ϵa 2 + 1 ϵ b 2 and absorb the corresponding squares into the quadratic terms in (31) by adjusting the process parameters σ, µ appropriately, which provides the desired result in this case.
We now show the equivalence of the iteration process in ( 4) and ( 6).Theorem 4. Let the optimization penalties κ ≥ 0, µ > 0 and the initial condition y 0 ∈ L 2 (Ω) be given.Assume that we have the initial transmission data λ 0 ij , ρ 0 ij ∈ L 2 (0, T), i, j = 1, 2. Then, the domain decomposition procedure for the optimality system in (4) is equivalent to the sequence of virtual optimal control problems in (6).
Proof.Recalling the cost function in (6) and the corresponding Lagrange function We integrate by parts in order to reveal the Robin boundary conditions for y i .Taking variations, we arrive at the same adjoint equation as in (4) (with the actual iteration index deleted for the sake of simplicity).The corresponding resulting Robin boundary condition for p i reads as follows: and the action of the virtual control h ij is revealed by the variation of L i with respect to ĥij , i.e., ∂ h ij L i (y i , u i , h ij , p i )( ĥij ) = 0, which leads to Putting these equations together provides the Robin conditions in (4).After each optimization step, the coupling parameters λ n+1 ij , ρ n+1 ij need to be evaluated by the update rule (5).In the given case of purely distributed controls, the adjoint variable needed for the update is simply provided by νu i .

Finite Difference Approximation
We now provide a numerical scheme for solving the resulting optimality system and the corresponding decomposed system (4) using an implicit finite difference method.We change the sign of p, which, as the system is linear, does not change the solutions, only the signs of the controls and the coupling.Hence, in this case, the optimality system ( 15)-( 17) reads as follows.
Remark 6.We note that the above estimate (41) holds only for smooth functions, i.e., when y ∈ C 2 [0, t m ].However, due to the weakly singular kernel of the Caputo fractional derivative, the solution of fractional parabolic equations generally does not possess a smooth solution throughout the closed domain even in the case of smooth input data.In particular, we recall from Stynes et al. [39] (Theorem 2.1) that the typical solution of time-fractional parabolic equations (TFPEs) involving Caputo fractional derivatives possesses a singularity near the initial time t = 0. Thus, by taking into account this initial singularity, the authors of [39] proposed the L1 method on a graded mesh and obtained the following error estimate for y ∈ C[0, t m ] ∩ C 2 (0, t m ]: where t m ∈ Ω τ := {t m : t m = T m M r , 0 ≤ m ≤ M}, with r ≥ 1 as the mesh grading constant. Therefore, the order of accuracy of the L1 method can be obtained according to (41) in the case of a smooth solution and (42) in the case of a non-smooth solution.However, as the emphasis of this article is to introduce (for the first time) a non-overlapping domain decomposition procedure for such problems and not to dwell on the numerical analysis, for the sake of simplicity and numerical comparison we have only considered the discrete approximation of the Caputo fractional derivative according to (41), i.e., for those solutions which are sufficiently smooth in time.However, in the case of a non-smooth solution (initial singularity) of TFPEs on metric graphs, a graded mesh could be employed to obtain the desired order of accuracy in time while following the spatial discretization according to the one carried out here, which we will explicitly describe in a forthcoming article.
As for the spatial operator, we use the standard the second-order central finite difference for approximation of the spatial derivative and denote Y m+1 i,r and P m+1 i,r as the approximations for the state variable y i (x i,r , t m+1 ) and adjoint variable p i (x i,r , t m+1 ), respectively.Therefore, using (40), we obtain the following finite difference scheme for the state equation in the optimality system (39): with the following additional conditions: • initial conditions Y i,r (0) = y 0 i (x r ), r = 0, 1, . . ., R, i = 1, 2; (44) • boundary conditions at (x = ℓ i ) • continuity condition at (x = 0) with an unknown function Ȳ; and • discretization of the Kirchhoff condition at x = 0, as follows: Now, considering Equation (2) 1 at x = 0 while approximating the fractional derivative by (40) and y ′′ i (0) using the above formula, we obtain Finally, summing up over all edges in the above expression and using Equation (2) 4 , we obtain where P m 1,0 = P m 2,0 = Pm .• Robin-type conditions: − ∂ ∂x y i (0, t) + σy i (0, t) = g i (t), i = 1, 2, in a similar manner as above; hence, we obtain Regarding the accuracy of the proposed difference scheme (43)-(47) for approximating the state Equation ( 2), we have the following result (see, [40] (Theorem 3.2)).Theorem 5. Let y i (x i,r , t m+1 ) be the exact solution of (2) at grid point (x i,r , t m+1 ) on the Interval I i , i = 1, 2, and let y m+1 i,r be the solution of the difference scheme (43)-(47); then, for sufficiently small ∆t and ∆x i , we have where , e m+1 i,r Next, in order to obtain the difference scheme for the adjoint equation, we introduce where L α T denotes the discrete right-sided fractional differential operator.Hence, in view of (49) and following the same algorithm (43)-(47) as above, we obtain the following discrete version of the optimality system.
The discrete optimality system (50) is solved for the value of Y i and P i at all the grid points for i = 1, 2. Finally, based on the above discretization, we now introduce the discrete decoupling scheme.For the sake of simplicity, we ignore the iteration index at the state and adjoint variable, and instead keep only the index at the iteration histories λ ij , ρ ij , obtaining (51)

Examples
We consider two sets of examples, one based on the direct decomposition of the optimality system according to (4) and one on the iteration of virtual optimal control problems (6).
Example 1.In the first example, we consider an exact solution for comparison.We rearrange the intervals I i = (0, 0.5), i = 1, 2 such that their union, after a proper change of variables, can be represented by Ω = (0, 1).We take y exact (t, x) = t sin(πx); correspondingly, we take the right-hand side f (t, x) = ( Γ(2) Γ(2−α) t 1−α + tπ 2 ) sin(πx).Now, we take as the target the solution z d = y exact , meaning that the adjoint is effectively zero everywhere and the optimal state is the exact solution.We take the parameters κ = 1e.4,ν = 1, µ = 0.5, σ = 0.2, ϵ = 0.5 (see Remark 2) and start with zero initial data.We start the iterations with zero values.For the discretization, we choose Nx = 32 spatial points and Nt = 100 time points.In the first part of the experiment, we take the domain decomposition algorithm based on the optimality system (DD for the OS) (4), (5), while in the second part of the experiment we take the domain decompostion for the virtual optimal control problem (DD for the virtual OCP) (6), (5).In both cases, we take the same process parameters in order to ensure a good comparison.In the first part we solve the local optimality systems in each iteration, while in the second part we invoke an optimization routine for the corresponding virtual control problem.For the latter, we take IPOPT as the optimizer, which uses the MUMPS linear solver.The package we use is CASADI-OPT https://web.casadi.org/(accessed on 10 January 2024).Clearly, the optimal choice of the parameters σ, µ depends on the penalty parameters κ, ν, and the choice may differ depending on the method; here, we keep both sets in order to make for a better comparison.Of course, when using the optimization approach for the virtual control problem, the IPOPT solver again employs a number of parameters, which makes a very precise numerical comparison of this matter difficult.In Figure 1, we show on the left (Figure 1a) the states on the two sub-domains and on the right (Figure 1b) the states at the interface x = 0.5.As a result, there is no visible difference between the results obtained from the iterations (4) and (6).The corresponding errors of the states and the derivatives are shown in Figure 2.
Example 2. In the second example, we take the target z d = 1.The sub-domains I i = (0, 1), i = 1, 2 are such that the overall domain can be represented by Ω = (−1, 1) after an appropriate change of variables.Here, using the same parameter as in the previous example, we compare the solutions with the computed solutions of the global optimality system.See Figure 3 for the states and adjoints and Figure 4 for the states and errors at the interface.In Figure 5, we show the corresponding errors for the virtual control problem.It is interesting to see that for high penalties the fully distributed control concentrates on the action where the initial and boundary data are prescribed.In addition, it is apparent that the errors of the iterations at the interface decline more slowly after a few iterations due to the accuracy of the finite difference solution.

Conclusions and Future Work
We have introduced a non-overlapping domain decomposition method for optimal control problems associated with time-fractional diffusion equations in one space dimension on a two-link domain.By decomposing the resulting optimality system of the considered optimal control problem by an iterative procedure, we have obtained the iteration of optimal control problems on the individual edges (intervals) and solved the alternating sequence of Robin-type virtual optimal control problems.In addition, we have shown the consistency and convergence of the proposed algorithm and demonstrated the equivalence between the domain decomposition procedure for the optimality system and the sequence of virtual optimal control problems.Finally, we have discussed the discretization of the resulting optimality system and the corresponding decomposed system of the considered problem, and shown some numerical evidence.In future, we will consider the non-overlapping domain decomposition method for optimal control problems with boundary controls on general metric graphs.

Figure 2 .
Figure 2. Errors at the interface: example 1.(a) Errors at the interface for the DD of the OS.(b) Errors at the interface for the DD of the virtual OCP.

Figure 5 .
Figure 5. Errors at the interface for the DD of the virtual OCP, different discretizations: example 2. (a) Errors at interface Nx = 10.(b) Errors at interface Nx = 40.