A posteriori error estimates of spectral method for the fractional optimal control problems with non-homogeneous initial

Abstract: In this paper we consider an optimal control problem governed by a space-time fractional diffusion equation with non-homogeneous initial conditions. A spectral method is proposed to discretize the problem in both time and space directions. The contribution of the paper is threefold: (1) A discussion and better understanding of the initial conditions for fractional differential equations with Riemann-Liouville and Caputo derivatives are presented. (2) A posteriori error estimates are obtained for both the state and the control approximations. (3) Numerical experiments are performed to verify that the obtained a posteriori error estimates are reliable.


Introduction
The subject of fractional calculus has gained considerable popularity and importance during the past three decades, due to, on one side, its well-recognized applicability in diverse and widespread fields of science and engineering, and on the other side, its attractive complementary of the integer order calculus in mathematical sciences. Interested reader in fractional calculus and fractional partial differential equations (FPDEs) can refer to [1][2][3].
In recent years the research of optimal control problem governed by FPDEs has formed a hot topic not only in theoretical investigation, but also in numerical methods. In [4,5] the authors discussed the wellposedness of optimal control problem governed by time fractional diffusion equation with control constraint and state constraint, respectively. The authors in [6][7][8][9] studied the finite element method for optimal control problem governed by the fractional Laplace equation, the prior error estimate and the posterior error estimate were derived. A spectral Galerkin approximation of optimal control problem governed by Riesz fractional differential equation was investigated in [10]. The integral definition of the fractional Laplacian and a linear-quadratic optimal control problem for the fractional heat equation with control constraints were analyzed in [11]. Finite element method for optimal control problem for a fractional semilinear PDE with control constraints was studied in [12]. Numerical method for solving a class of fractional optimal control problems with canonical equality and inequality constraints was considered in [13,14]. Li et al. [15] presented a numerical method for solving nonlinear fractional optimal control problems with multiple states and controls. A spectral optimization methods for the time fractional diffusion inverse problem was discussed in [16]. Spectral Galerkin approximation of time fractional optimal control problems with integral constraint on the state was studied in [17]. Du et al. [18] studied the optimal control problem of fractional diffusion equations based on the finite difference fast projection gradient algorithm. Wu and Huang [19] introduced a time parallel algorithm to solve the optimal control problem of fractional diffusion equation in order to reduce the computation time. Also we refer the interested reader in optimal control problem for FPDEs to see [20][21][22][23][24][25] and its references for some recent works in the subject.
To the best of our knowledge, only few attempts have been made to develop the posteriori error estimation for problems involving fractional derivative, especially the one for optimal control problem governed by FPDEs. An a posteriori error analysis for an optimal control problem involving the fractional Laplacian was introduced in [7]. An a posteriori error estimator for PDE constrained optimization problem involving a nondifferentiable cost functional, fractional diffusion, and control-constraints was proposed and analyzed in [26]. In an earlier work of Ye and Xu [27], a posteriori error estimates for spectral method was carried out for the fractional optimal control problem governed by the time fractional diffusion equation with Riemann-Liouville fractional derivative. It is known within the community of fractional calculus that with the Riemann-Liouville definition, only homogeneous conditions (under an integral form) is possible. As an extension of the above mentioned paper [27], in this work we shall mainly be interested in a posteriori error estimates for the fractional optimal control problem with non-homogeneous conditions. As we are going to see, the non-homogeneous case presents a considerable difference on the solution property, and introduces additional non-trivial difficulties into the analysis of the a posteriori error of the optimal control problem.
The literature reports several definitions for the fractional derivative, including the expressions derived by Riemann-Liouville, Caputo, Riesz, Riesz-Caputo, Grunwald-Letnikov, Hadamard, and so on. Although the way to impose suitable initial conditions for different definitions has been discussed in the literature, there are still some issues need to be clarified. We will focus on the initial value problem for Riemann-Liouville and Caputo fractional differential equations in Section 2, where some new initial conditions are derived. Besides, as compared to the posteriori error analysis in [27], the presence of a singular term t −α in the right hand side of the state equation associated to the non-homogeneous initial condition leads to additional difficulties. Hence special techniques such as Embedding Theorem should be adopt to achieve reliable a posteriori error estimator.
The rest of this paper is organized as follows. In Section 2, we discuss the issue about the initial condition for Riemann-Liouville and Caputo fractional ordinary differential equation, and propose a consistent way to impose the "well-defined" initial condition. The formulation of the optimal control problem with non-homogeneous initial conditions and the optimality conditions are presented in Section 3. Section 4 describes the spectral discretization of the optimal control problem. In Section 5, we establish the main result on the posteriori error analysis for the considered problem. Numerical experiment is presented in Section 6, showing that the posteriori error estimates are reliable. Some concluding remarks are given at the end of the paper.

Issue on the initial condition
In this section we discuss the issue how to correctly impose initial conditions for the time fractional differential equations involving Riemann-Liouville and Caputo fractional derivatives, respectively.
Let I = (0, T ), T > 0. We begin with the following fractional differential equation: where R D α t is the Riemann-Liouville fractional differential operator of order α given by The following initial condition of integral form for a given u 0 should be provided to guarantee the uniqueness of the solution, where 0 < α < 1, I α t is the Riemann-Liouville fractional integral operator of order α defined by In fact, a direct calculation using [3,Lemma 5.2] gives the expression of the solution to (2.1): Clearly, the solution uniquely exists subject to the condition (2.2). This is in agreement with the common sense in two extreme cases: α → 1 and α → 0. In the former case, it was shown [1]: Thus the Eq (2.1) becomes and the initial condition (2.2) becomes In the latter case, i.e., α → 0, it can be easily verified from (2.3) that u(t) → f (t) for integrable f . One of the main inconveniences with the Eq (2.1) is that there will be not much choice in imposing a well-defined initial condition u 0 for 0 < α < 1. In fact, it has been shown in [28] the following lemma.
In order to be able to use non-homogeneous initial conditions, a common practice is to consider the fractional equation with Caputo derivative as follows: where C D α t is the Caputo differential operator defined by It has been widely believed that the Eq (2.6) subject to the standard form of the initial condition is well-posed. However, this may not be true if we are interested in weak solutions or even singular solutions. For example, if we consider the right hand side function f (t) = t δ−1/2 ∈ L 2 (I), δ > 0, it is not difficult to check that the Eq (2.6) admits a unique solution which belongs to H α+δ−ε (I) for any ε > 0, but has non finite values at t = 0 for 0 < α < 1/2, 0 < δ < 1/2 − α. In fact the initial condition (2.7) can only be made sense for f ∈ H δ (I) with δ > 1/2 − α.
For the purpose to impose the initial condition in a consistent way and allow for weaker solutions, we propose to consider the fractional differential equation under Riemann-Liouville definition: It is seen in the above equation that an "initial condition" is imposed in a weak sense in the right hand side, similar to the Neumann condition for a second order elliptic problem. This is motivated by the relationship which is true for all functions v such that all terms in the above identity are well-defined. A direct computation shows that the Eq (2.8) subject to the initial condition

Formulation of the problem and optimization
Inspired by the discussion in the previous section, we are interested in an optimal control problem governed by the space-time fractional diffusion equation (STFDE) with non-homogeneous initial condition under the form: We first introduce some notations that will be used throughout the paper. In all that follows, we use the expression A B to denote A ≤ cB with c being a generic positive constant. For a domain O, which may be Λ, I or Ω = Λ × I, we recall that L 2 (O), H s (O), and H s 0 (O) stand for the usual Sobolev spaces, equipped with the norms · 0,O and · s,O respectively. Let X be the Sobolev space with norm · X , we introduce the space Particularly, when X stands for H µ (Λ) or H µ 0 (Λ), the norm of the space H s (I; X) will be denoted by · µ,s,Ω . Hereafter, in cases where no confusion would arise, the domain symbols I, Λ, Ω may be dropped from the notations.
We now consider the following optimal control problem: where g and h are given convex functionals, the state variable u is the solution of STFDE (3.1), and the control variable q satisfies In order to formulate the weak form of the state Eq (3.1), we recall the state space which can be found in Li and Xu [28]. Then we have the following weak formulation of the state Eq (3.1): with the bilinear form A(·, ·) being defined by T respectively denote the left and right Riemann-Liouville fractional derivative of order α 2 , R a ∂ β 2 x and R x ∂ β 2 b respectively denote the left and right Riemann-Liouville fractional derivative of order β 2 . It has been shown in [28] that the following continuity and coercivity hold and the problem (3.3) is well-posed. Let K denote the admissible set associated to the constraints (3.2b) as then the optimal control problem (3.2) reads: find (q * , u(q * )) ∈ K × B α 2 , β 2 (Ω), such that Assume that the functional g(·) is bounded below, and h(q) → +∞ as q 0,Ω → +∞, it is known (see, e.g., [29,30]) that the optimal control problem (3.5) admits a unique solution (q * , u(q * )) ∈ K × B α 2 , β 2 (Ω). The well-posedness of the state equation ensures the existence of a control-to-state mapping q → u = u(q) defined through (3.3). Then the optimal control problem (3.5) equivalently can be rewritten in the following form: find q * ∈ K, such that for the reduced functional J(q) := J(q, u(q)) over L 2 (Ω). The first order necessary optimality conditions for (3.6) take the form which is proved in the following lemma.
where z(q) = z is the solution of the following adjoint state equation with t ∂ α T (0 < α < 1) denoting the right Caputo fractional derivative. Proof. The proof of this lemma is very close to that given for Lemma 2.1 in [21], but for the reader's convenience, we therefore give the details as follows. It is apparent from the chain rule that We now in a position to compute u (q)(ζ). For simplicity, let δu denote the derivative of u = u(q) in the direction ζ, that is Then it is readily seen that δu is the solution of the following problem: (3.13) To prove (3.8), we multiply each side of (3.9) by δu, then integrate the resulted equation on the domain Ω to yield Taking into account the boundary conditions in (3.13) and (3.11), by means of the fractional integration by parts demonstrated in [28], it holds Once again we use the fractional integration by parts to find Noticing the terminal condition in (3.10) and the initial condition in (3.13), we get Finally, combining (3.14), (3.15), and (3.17), we obtain This, together with (3.12), leads to (3.8).
The weak form of (3.9)-(3.11) reads: find z ∈ B . In what follows we will need the mapping q → u(q) → z(q), where for any given q, u(q) is defined by (3.3), and once u(q) is known, z(q) is defined by (3.19).

Space-time spectral discretization
We define the polynomial space where P M and P N respectively denote the space of polynomials with degree no more than M and N, L stands for the parameter pair (M, N).
Then we arrive at the spectral approximation to the state Eq (3.3) as follows: find u L (q) ∈ S L such that According to [28], the above discrete state equation admits a unique solution. The following estimate, also derived in [28], will be used in the analysis later on.
Suppose u ∈ H α 2 (I; H µ (Λ)) ∩ H ν (I; H β 2 0 (Λ)), 0 < α < 1, 1 < β < 2, ν > 1, µ ≥ 1, then we have We now define the semidiscrete optimal problem: find q * ∈ K, such that where for q ∈ L 2 (Ω) we set J L (q) := J(q, u L (q)) with u L (q) being defined in (4.1). The unique solution q * of problem (4.3) is characterized by the variational inequality where similar as Eq (3.8) (4.5) The above semidiscrete adjoint state z L (q) ∈ S L is the solution of the following problem: A(ϕ L , z L (q)) = (g (u L (q)), ϕ L ) Ω , ∀ϕ L ∈ S L . (4.6) By defining the following finite dimensional subspace for the control variable: we arrive at the full discrete optimal control problem: find q * L ∈ K L , such that where the unique solution q * L satisfies the first order optimality condition (4.8)

A posteriori error estimates
In order to obtain the posteriori error estimates for the spectral approximation (4.7), we need some approximation operators and their approximation properties. Let r and s be two real numbers such that r n + 1 2 , 0 ≤ s ≤ r, there exits an operator Π s,0 r,N : Let Π N be the L 2 -orthogonal projection operator from L 2 (I) onto P N (I). Equivalently, it means that, for any function v ∈ L 2 (I), Π N v belongs to P N (I) and satisfies (Π N v − v, w N ) I = 0, ∀w N ∈ P N (I).
Then for any nonnegative real number m, the following estimate holds The L 2 -orthogonal projector Π M in Λ is defined similarly. We now proceed to analyze the approximation error of the proposed space-time spectral method. The proof of the main result will be accomplished with a series of lemmas which we present below.

Lemma 5.1. If g(·) is convex, and h(·) is strictly convex such that
Then for all p, q ∈ L 2 (Ω), we have Proof. The proof of this result is quite similar to that given for Lemma 4.2 in [21] and so is omitted.
Proof. Using (3.7), (5.3) and (4.8), we can show that for arbitrary p L ∈ K L , By substituting (5.6) into (5.5), and using the Young inequality, we get where ε > 0 is a sufficient small constant. Then, we choose ε = c 4 to derive Finally, by using the approximation estimate (5.2) we obtain (5.4).
Remark 5.1. The above results can be simplified in the case where the cost functional is quadratic, i.e., whereū is a given observation data. In this case, we have Lemma 5.3. Let q * L denote the solution of (4.7) with the associated discrete state u L (q * L ) and adjoint state z L (q * L ). Assume g (·) is Lipschitz continuous. Then it holds A(e z , e z ) = A(e z − e L z , e z ) + A(e L z , e z ) = A(e z − e L z , z L (q * L )) − A(e z − e L z , z(q * L )) + A(e L z , e z ) ,Ω e z 0,Ω .

(5.13)
On the other hand, Hence, if 0 < α < 1 2 , we take and if 1 2 ≤ α < 1, we set As a result, it holds for the above chosen parameters that . (5.20) Taking (5.14) into consideration, we get With the help of the preceding three lemmas we can now prove the following main result.
Theorem 5.1. Let q * be the solution of (3.6) with the associated state u(q * ) and adjoint state z(q * ). q * L is the solution of (4.7) with the associated discrete state u L (q * L ) and adjoint state z L (q * L ). Under the assumptions in Lemmas 5.2, 5.3, and 5.4, the following estimate holds:
Proof. We begin with estimating q * − q * L 0,Ω . According to (5.4), it suffices to estimate z L (q * L ) − z(q * L ) 0,Ω . Combining (5.8) and (5.11) we obtain Using the above estimate, inequality · 0,Ω · B α 2 , β 2 (Ω) , and Lemma 5.2, we get It follows from (3.3) and (3.19) that u(q * ) − u(q * L ) satisfies and z(q * ) − z(q * L ) satisfies A(ϕ, z(q * ) − z(q * L )) = (g (u(q * )) − g (u(q * L )), ϕ), ∀ϕ ∈ B Taking v = u(q * ) − u(q * L ) in (5.25) gives Furthermore, using the triangle inequalities and (5.11), we get Similarly, taking ϕ = z(q * ) − z(q * L ) in (5.26) we have In view of the Lipschitz continuity of g (·), we obtain It follows from (5.23), (5.29) and the triangle inequality: Remark 5.2. In Theorem 5.1, we have obtained an a posteriori upper bound, expressed by the reliable indicator η, which will be checked in our numerical experiments. For the time being we are unable to derive an optimal estimate. As for the classical elliptic and parabolic equations, we guess such an optimal estimate will have to make use of some Jacobi-weighted Sobolev spaces and polynomial inverse inequalities. We leave this issue for future work.

Numerical results
We carry out in this section a series of numerical experiments and present results to validate the obtained error estimates. We use the projection gradient optimization algorithm proposed in [27] to solve the optimization problem. In all the calculations, we take Λ = (−1, 1), T = 1, and We consider the problem (3.2) with the exact analytical solutions as follows u(q * ) = sin πx cos πt, z(q * ) = sin πx sin π(1 − t), q * = max{0, z(q * )} − z(q * ).
The right-hand side function f and the observation stateū are numerically calculated through (3.3) and (3.19) using u(q * ), z(q * ), and q * .  We first investigate the impact of the initial guess on the convergence of the projection gradient optimization algorithm. We start by considering q (0) = 25q * . In Figure 1, we present the convergence history of the gradient of the objective function versus the iteration number with M = 18, N = 18, α = 0.6, β = 1.6. We see that the iterative method converges within a dozen iterations. We then take q (0) to be constant c with c = 0 or 15, and repeat the same computation as the previous test. The result is given in Figure 2. These results seem to tell that the initial guess has no significant effects on the convergence of the projection gradient iterative algorithm. In the following, the initial guess is set to 0.
Next we check the convergence behavior of numerical solutions with respect to the polynomial degrees M. In Figures 3 to 5, we plot the errors as functions of the polynomial degrees M with α = 0.6, N = 20 for β = 1.2, 1.6, 1.9. The errors show an exponential decay, since in this semilog representation one observes that the error variations are essentially linear versus the degrees of polynomial.   We then investigate the temporal errors, that is, for fix M, we check the convergence behavior of numerical solutions with respect to the polynomial degrees N in time direction. In Figures 6 to 8, we plot the errors as functions of N with β = 1.6, M = 20 for three values α = 0.2, 0.6, 0.9. The straight line of the error curves indicates that the convergence in time is also exponential.   We observe that the indicator η has almost the same exponential decay rate as the real error e. Also observed in the figures is that the a posteriori error is an upper bound. This is consistent with the theoretical results.

Conclusions
In the present work, we have first discussed the commonly used initial conditions for differential equations involving Riemann-Liouville and Caputo derivatives, and proposed a consistent way to impose the non-homogeneous initial condition. Then an optimal control problem governed by the space-time fractional diffusion equation with non-homogeneous initial conditions was investigated. A posteriori upper bound was derived for the space-time spectral approximation. Numerical experiment has been carried out to confirm the theoretical results. In a future work, we will focus on the efficiency of the error estimates, i.e. a posteriori error lower bound, to obtain equivalent a posteriori error estimates for the considered optimal control problems. It will be also interesting to look for a strategy that allows to approximate the optimal control problem with adaptive spectral method.