N ov 2 01 5 FUNCTIONAL A POSTERIORI ERROR ESTIMATES FOR TIME-PERIODIC PARABOLIC OPTIMAL CONTROL PROBLEMS

This paper is devoted to the a posteriori error analysis of multiharmonic finite element approximations to distributed optimal control problems with time-periodic state equations of parabolic type. We derive a posteriori estimates of functional type, which are easily computable and provide guaranteed upper bounds for the state and co-state errors as well as for the cost functional. These theoretical results are confirmed by several numerical tests that show high efficiency of the a posteriori error bounds.


Introduction
We consider optimal control problems with time-periodic parabolic state equations.Problems of this type often arise in different practical applications, e.g., in electromagnetics, chemistry, biology, or heat transfer, see also [1,2,15].Moreover, optimal control problems are the subject matter of lots of different works, see, e.g., [30,14,36,6] and the references therein.The multiharmonic finite element method (MhFEM) is well adapted to the class of parabolic time-periodic problems.Within the framework of this method, the state functions are expanded into Fourier series in time with coefficients depending on the spatial variables.In numerical computations, these series are truncated and the Fourier coefficients are approximated by the finite element method (FEM).This scheme leads to the MhFEM (also called harmonic-balanced FEM), which was successfully used for the simulation of electromagnetic devices described by nonlinear eddy current problems with harmonic excitations, see, e.g., [38,3,4,5,10] and the references therein.Later, the MhFEM has been applied to linear time-periodic parabolic boundary value and optimal control problems [17,18,24,28,37] and to linear time-periodic eddy current problems and the corresponding optimal control problems [19,20,21].In the MhFEM setting, we are able to establish inf-sup and sup-sup conditions, which provide existence and uniqueness of the solution to parabolic timeperiodic problems.For linear time-periodic parabolic problems, MhFEM is a natural and very efficient numerical technology based on the decoupling of computations related to different modes.
This paper is aimed to make a step towards the creation of fully reliable error estimation methods for distributed time-periodic parabolic optimal control problems.We consider the multiharmonic finite element (MhFE) approximations of the reduced optimality system, and derive guaranteed and fully computable bounds for the discretization errors.For this purpose, we use the functional a posteriori error estimation techniques earlier introduced by S. Repin, see, e.g., the papers on parabolic problems [32,13] as well as on optimal control problems [11,12], the books [33,29] and the references therein.In particular, our functional a posteriori error analysis uses the techniques close to those suggested in [32], but the analysis contains essential changes due to the MhFEM setting.In [27], the authors already derived functional a posteriori error estimates for MhFE approximations to parabolic time-periodic boundary value problems.Similar results are now obtained for the MhFE approximations to the state and co-state, which are the unique solutions of the reduced optimality system.It is worth mentioning that these a posteriori error estimates for the state and co-state immediately yield the corresponding a posteriori error estimates for the control.In addition to these results, we deduce fully computable estimates of the cost functional.In fact, we generate a new formulation of the optimal control problem, in which (unlike the original statement) the state equations are accounted in terms of penalties.It is proved that the modified cost functional attains its infimum at the optimal control and the respective state function.Therefore, in principle, we can use it as an object of direct minimization, which value on each step provides a guaranteed upper bound of the cost functional.
The paper is organized as follows: In Section 2, we discuss the time-periodic parabolic optimal control problem and the corresponding optimality system.The multiharmonic finite element discretization of this space-time weak formulation is considered in Section 3. Section 4 is devoted to the derivation of functional a posteriori error estimates for the optimality system associated with the multiharmonic setting.In Section 5, we present new results related to guaranteed and computable bounds of the cost functional.In the final Section 6, we discuss some implementation issues and present the first numerical results.

A Time-Periodic Parabolic Optimal Control Problem
Let Q T := Ω × (0, T ) denote the space-time cylinder and Σ T := Γ × (0, T ) its mantle boundary, where the spatial domain Ω ⊂ R d , d = {1, 2, 3}, is assumed to be a bounded Lipschitz domain with boundary Γ := ∂Ω, and (0, T ) is a given time interval.Let λ > 0 be the regularization or cost parameter.We consider the following parabolic time-periodic optimal control problem: subject to the parabolic time-periodic boundary value problem Here, y and u are the state and control functions, respectively.The coefficients σ(•) and ν(•) are uniformly bounded and satisfy the conditions where σ, σ, ν and, ν are constants.As usual, the cost functional (1) contains a penalty term weighted with a positive factor λ.This term restricts (in the integral sense) values of the control function u.
We can reformulate the problem (1)-( 2) in an equivalent form.For this purpose, we introduce the Lagrangian which has a saddle point, see, e.g., [14,36] and the references therein.The proper sets for which this saddle point problem is correctly stated are defined later.Since the saddle point exists, the corresponding solutions satisfy the system of necessary conditions Using the second condition, we can eliminate the control u from the optimality system (5), i.e., From ( 6) it appears very natural to choose u and p from the same space.Moreover, we arrive at a reduced optimality system, written in its classical formulation as In order to determine generalized solutions of (7), we define the following spaces (here and later on, the notation is similar to that was used in [25,26]): which are endowed with the norms , respectively.Here, ∇ = ∇ x is the spatial gradient and ∂ t denotes the generalized derivative with respect to time.Also we define subspaces of the above introduced spaces by putting subindex zero if the functions satisfy the homogeneous Dirichlet condition on Σ T and subindex per if they satisfy the periodicity condition v(x, 0) = v(x, T ).All inner products and norms in L 2 related to the whole space-time domain Q T are denoted by Here, are the Fourier coefficients, T denotes the periodicity and ω = 2π/T is the frequency.Since the problem has time-periodical conditions, these representations of the exact solution and respective approximations are quite natural.In what follows, we also use the spaces where t v is defined in the Fourier space by the relation Here, v k = (v c k , v s k ) T for all k ∈ N, see also [28].These spaces can be considered as Hilbert spaces if we introduce the following (equivalent) products: The above introduced spaces allow us to operate with a "symmetrized" formulation of the problem (7) presented by (15).The seminorm and the norm of the space H 1, 1 2 per (Q T ) are defined by the relations respectively.Using Fourier type series, it is easy to define the function "orthogonal" to v: Obviously, u ⊥ k Ω = u k Ω and we find that Henceforth, we use the following subsidiary result (which proof can be found in [27]): are valid for all y ∈ H 0,1 per (Q T ) and v ∈ H Also, we recall the orthogonality relations (see [28,37]) σ∂ t y, y = 0 and σy ⊥ , y = 0 t y ⊥ = 0 and ν∇y, ∇y and the identity where ξ, ∂ We note that for functions presented in terms of Fourier series the standard Friedrichs inequality holds in Q T .Indeed, In order to derive the weak formulations of the equations (7), we multiply them by the test functions z, q ∈ H 1, 1 2 0,per (Q T ), integrate over Q T , and apply integration by parts with respect to the spatial variables and time.We arrive at the following "symmetric" space-time weak formulations of the reduced optimality system (7): Given the desired state y d ∈ L 2 (Q T ), find y and p from for all test functions z, q ∈ H 1, 1 2 0,per (Q T ).We can represent (15) in a somewhat different form.For this purpose, it is convenient to introduce the bilinear form Then (15) reads as B((y, p), (z, q)) = (y d , 0), (z, q) ∀ (z, q) ∈ (H

Multiharmonic Finite Element Approximation
In order to solve the optimal control problem (1)-( 2) approximately, we discretize the optimality system (15) by the MhFEM (see [28]).Using the Fourier series ansatz (8) in (15) and exploiting the orthogonality of cos(kωt) and sin(kωt), we arrive at the following problem: Find for all test functions z k , q k ∈ V.The system (17) must be solved for every mode k ∈ N.For k = 0, we obtain a reduced problem: Find for all test functions z c 0 , q c 0 ∈ V .The problems ( 17) and ( 18) have unique solutions.In order to solve these problems numerically, the Fourier series are truncated at a finite index N and the unknown Fourier coefficients In this work, we use continuous, piecewise linear finite elements on a regular triangulation T h to construct V h and its basis (see, e.g., [7,9,16,35]).This leads to the following saddle point system for every single mode k = 1, 2, . . ., N : which has to be solved with respect to the nodal parameter vectors . The matrices M h , M h,σ , and K h,ν denote the mass matrix, the weighted mass matrix and the stiffness matrix, respectively.Their entries are defined by the integrals and the right hand side vectors have the form and y s dk = Ω y s dk ϕ j dx j=1,...,n .
For k = 0, the problem ( 18) generates a reduced system of linear equations, i.e., Fast and robust solvers for the systems ( 19) and ( 20) can be found in [18,23,28,37], which we use in order to obtain the MhFE approximations

Functional A Posteriori Error Estimates for the Optimality System
Now we are concerned with a posteriori estimates of the difference between the exact solution (y, p) and the respective finite element solution (y N h , p N h ).First, we present the inf-sup and sup-sup conditions for the bilinear form (16).
Proof.The right hand-side inequality in (23) results from the triangle and Cauchy-Schwarz inequalities and, the Friedrichs inequality (14).Indeed, , where μ2 = max{1, 1 λ , ν, σ} max{1, C 2 F + 1}.The left hand-side inequality in (23) is proved quite similarly to the previous case.We select the test functions and use the σ-and ν-weighted orthogonality relations (11).Then, we find that In view of the estimate . This justifies the left hand-side inequality in (23).
Let (η, ζ) be an approximation of (y, p), which is a bit more regular with respect to the time variable than the exact solution (y, p).Namely, we assume that η, ζ ∈ H 1,1 0,per (Q T ) (this assumption is of course true for the MhFE approximations y N h and p N h defined in ( 21)).Our goal is to deduce a computable upper bound of the error e := (y, p) (15), it follows that the integral identity holds for all z, q ∈ H 1, 1 2 0,per (Q T ).The linear functional It can be viewed as a quantity measuring the accuracy of (24) for any pair of test functions (z, q).Therefore, getting an upper bound of the error is reduced to the estimation of We reconstruct F (η,ζ) using the identity Also, we use the identities which hold for any z, q ∈ H 1 0 (Ω) and any For ease of notation, the index x in div x will be henceforth omitted, i.e., div = div x denotes the generalized spatial divergence.Using the Cauchy-Schwarz inequality yields the estimate Applying (14), we find that Hence, and by the inf-sup condition in (23), we obtain Thus, we arrive at the following result: , and the bilinear form B(•, •) defined by ( 16) meets the inf-sup condition (23).Then, where e = (y, p) .
Remark 1.For computational reasons, it is useful to reformulate the majorants in such a way that they are given by quadratic functionals, see, e.g., [11].This is done by introducing parameters α, β, γ > 0 and applying Young's inequality.For the error majorant Finally, we note that the inf-sup condition (22) implies an estimate similar to (28) for the error in terms of and the bilinear form B(•, •) defined by ( 16) satisfies (22).Then, where e is defined in Theorem 1 and The multiharmonic approximations.Since the desired state y d belongs to L 2 it can be represented as a Fourier series.Henceforth, we assume that the approximations η and ζ of the exact state y and the adjoint state p, respectively, are also represented in terms of truncated Fourier series as well as the vector-valued functions τ and ρ, i.e., where all the Fourier coefficients belong to the space L 2 (Ω).In this case, and the L 2 (Q T )-norms of the functions R 1 , R 2 , R 3 and R 4 defined in (27) can be represented in the form, which exposes each mode separately.More precisely, we have where div Remark 2. Since y d is known, we can always compute the remainder term of truncation with any desired accuracy.
It is important to outline that all the L 2 -norms of R 1 , R 2 , R 3 and R 4 corresponding to every single mode k = 0, . . ., N are decoupled.It is useful to introduce quantities related to each mode.For k ≥ 1, we denote them by  28) and ( 29), respectively, can be represented in somewhat new forms that contain quantities associated with the modes, namely, Indeed, let the error majorants vanish.Then, , for all k = 1, . . ., N d , so that collecting the N d +1 harmonics leads to multiharmonic representations for η, ζ, τ and ρ of the form (33) satisfying the equations Since η and ζ also meet the boundary conditions, we conclude that η = y and ζ = p.

Functional A Posteriori Estimates for Cost Functionals of Parabolic Time-Periodic Optimal Control Problems
This section is aimed at deriving guaranteed and computable upper bounds for the cost functional and establishing their sharpness.This is important because, in optimal control, we cannot in general compute the (exact) cost functional since the (exact) state function is unknown.By using a posteriori estimates for the state equation we overcome this difficulty.Similar results for elliptic optimal control problems can be found, e.g., in [11,33].Let y = y(v) be the corresponding state to a control v.The cost functional J (y(v), v) defined in (1) has the form where We wish to deduce majorants for the cost functional J (y(u), u) of the exact control u and corresponding state y(u) by using some of the results presented in [27], which are obtained for the time-periodic boundary value problem (2).In [27], the following functional a posteriori error estimate for problem (2) has been proved: where and v is a given function in L 2 (Q T ).Now, adding and subtracting η in the cost functional J (y(v), v) as well as applying the triangle and Friedrichs inequalities yields the estimate Since we conclude that Together with (34) this leads to the estimate By introducing parameters α, β > 0 and applying Young's inequality, we can reformulate the estimate such that the right-hand side is given by a quadratic functional (the latter functional is more convenient from the computational point of view).We have where The majorant (35) provides a guaranteed upper bound of the cost functional, which can be computed for any approximate control and state functions.Moreover, minimization of this functional with respect to η, τ , v and α, β > 0 yields the exact value of the cost functional.This important result is summarized in the following theorem: Theorem 3. The exact lower bound of the majorant J ⊕ defined in (35) coincides with the optimal value of the cost functional of problem ( 1)-( 2), i.e., inf η∈H 1,1 0,per (QT ),τ ∈H(div,QT ) v∈L 2 (QT ),α,β>0 Proof.The infimum of J ⊕ is attained for the optimal control u, the corresponding state function y(u) and the exact flux (ν∇y(u)).In this case, R 1 and R 2 vanish, and for α tending to zero, values of J ⊕ tend to the exact value of the cost functional.
Corollary 2. From Theorem 3, we obtain the following estimate: Now, it is easy to derive a posteriori estimates for the cost functional in the setting of multiharmonic approximations.Let η be the MhFE approximation y N h to the state y.Since the control v can be chosen arbitrarily in (35), we choose a MhFE approximation u N h for it as well.More precisely, we can compute the MhFE approximation u N h for the control from the MhFE approximation p N h of the adjoint state, since u N h = −λ −1 p N h , by solving the optimality system, from which we obtain y N h as well.We now apply (37) and select η = y N h and v = u N h .Next, we need to make a suitable reconstruction of τ , which can be done by different techniques, see, e.g., [33,29] and the references therein.In the treatment of our approach it is natural to represent τ in the form of a multiharmonic function τ N h .Then the majorant where T .Note that the remainder term (31) remains the same in (38).Since all the terms corresponding to every single mode k in the majorant J ⊕ are decoupled, we arrive at some majorants J ⊕ k , for which we can, of course, introduce positive parameters α k and β k for every single mode k as well.Then the majorant (38) can be written as where α N +1 = (α 0 , . . ., α N +1 ) T , β N = (β 0 , . . ., β N ) T , and Next, we have to reconstruct the fluxes τ c 0h and τ kh for all k = 1, . . ., N , which we denote by τ kh = R flux h (ν∇y kh ).This can be done by various techniques.In [27], we have used Raviart-Thomas elements of the lowest order (see also [31,8,34]), in order to regularize the fluxes by a post-processing operator, which maps the L 2 -functions into H(div, Ω).Collecting all the fluxes corresponding to the modes together yields the reconstructed flux ).After performing a simple minimization of the majorant J ⊕ with respect to α N +1 and β N , we finally arrive at the a posteriori estimate where ᾱN+1 and βN denote the optimized positive parameters.
It is worth outlining that the majorant J ⊕ provides a guaranteed upper bound for the cost functional, and, due to Theorem 3, the infimum of the majorant coincides with the optimal value of the cost functional.
Remark 4. In this work, problems with constraints on the control or the state are not considered, but inequality constraints imposed on the Fourier coefficients of the control can easily be included into the MhFE approach, see [17], and, hence, may be considered in the a posteriori error analysis of parabolic time-periodic optimal control problems as well.

Numerical Results
We compute and analyze the efficiency of the above derived a posteriori estimates for different cases, namely, 1. the desired state is periodic and analytic in time, but not time-harmonic, 2. the desired state is analytic in time, but not time-periodic, and 3. the desired state is a non-smooth function in space and time.Note that convergence and other properties of numerical approximations generated by the MhFEM have been studied in [18,28] for the same three cases.The optimal control problem (1)-( 2) is solved on the computational domain Ω = (0, 1) × (0, 1) with the Friedrichs constant C F = 1/( √ 2π) using a uniform simplicial mesh and standard continuous, piecewise linear finite elements.The material coefficients are supposed to be σ = ν = 1.In the first two examples, we choose the cost parameter λ = 0.1, and in the third one, we choose λ = 0.01.Both choices are common.
Getting sharp error bounds requires an efficient construction of η, ζ, τ and ρ in order to compute sharp guaranteed bounds from the majorants.We choose MhFE approximations (30) for η and τ as well as for ζ and ρ.In order to obtain suitable fluxes τ , ρ ∈ H(div, Q T ), we reconstruct them by the standard lowest-order Raviart-Thomas (RT 0 -) extension of normal fluxes.We refer the reader to [27], where the authors have discussed this issue thoroughly.
In order to solve the saddle point systems (19) for k = 1, . . ., N and ( 20) for k = 0, we use the robust algebraic multilevel preconditioner of Kraus (see [22]) for an inexact realization of the block-diagonal preconditioners where D = √ λK h,ν + kω √ λM h,σ + M h , and in the minimal residual method, respectively.The preconditioners (43) and (44) were presented and discussed in [18,28,37].The numerical experiments were computed on grids of mesh sizes 16×16 to 256×256.The algorithms were implemented in C++, and all computations were performed on an average class computer with Intel(R) Xeon(R) CPU W3680 @ 3.33 GHz.Note that the presented CPU times t sec in seconds include the computational times for computing the majorants, which are very small in comparison to the computational times of the solver.
In where T = 2π/ω and ω = 1.This function is time-periodic and analytic, but not time-harmonic.Hence, the truncated Fourier series approximation of y d has to be computed for applying the MhFEM as presented in Section 3.For that, the Fourier coefficients of y d can be computed analytically, and, then, they are approximated by the FEM.Next the systems (19) and (20) are solved for all k ∈ {0, . . ., N }.We choose the truncation index N = 8.Finally, we reconstruct the fluxes by a RT 0 -extension and compute the corresponding majorants.The exact state is given by y(x, t) = e t sin(t) 3 sin(x 1 π) sin(x 2 π).In Table 1, we present the CPU times t sec , the majorants and J ⊕ 0 as defined in (40) as well as the corresponding efficiency indices and J ⊕ k as defined in (41) for k = 1 and, finally, the corresponding efficiency indices obtained on grids of different mesh sizes.Similar results are obtained for larger k as well, which is illustrated in Table 3.This table compares the results for the modes k ∈ {0, . . ., 8} computed on the 256 × 256-mesh and presents the overall functional error estimates.For that, the remainder term E N is precomputed exactly, see Remark 2. It can be observed that the values of the majorants M ⊕ k |•| and J ⊕ k decrease for increasing k, but that the values of the efficiency indices are all about the same, which is a demonstration for the robustness of the method with respect to the modes.Note that the overall efficiency index for N = 6 is large (I M eff = 3.15) compared to the efficiency indices corresponding to the single modes.The reason for that is that the remainder term E 6 = 640.25 is still quite large, and hence, more modes are needed.For N = 8, the remainder term is E 8 = 106.07,which leads to a much better overall efficiency index (I M eff = 1.69).The value for the cost functional is however sufficiently accurate with N = 6.This example demonstrates that the a posteriori error estimates clearly show what amount of modes would be sufficient for representing the solution with a desired accuracy.20) for all k ∈ {0, . . ., N } with first N = 6 and second N = 8 being the truncation index.Finally, we compute also the solutions for N = 10.The exact state is given by y(x, t) = e t sin(t) sin(x 1 π) sin(x 2 π).The results related to computational expenditures and efficiency indices are quite similar to those for Example 1.Therefore, we present only numerical results in the form similar to Table 3 (see Table 4).In this numerical experiment, we again observe good and satisfying efficiency indices for M ⊕ k |•| .The remainder terms for N = 6, N = 8 and N = 10 are E 6 = 44094.84,E 8 = 19869.30and E 10 = 10597.20,respectively.The efficiency index for M ⊕ |•| with N = 10 improves a lot compared to the index with N = 6.Note that -as in the first example -the efficiency indices for the cost functional are approximately one.This again demonstrates the accurateness of the majorants for the cost functional.
In Example 3, we set , where χ denotes the characteristic function in space and time.Let T = 1, then ω = 2π.Again the coefficients of the Fourier expansion associated with y d can be found analytically.They are ,1] 2 (x), and y s dk (x) = 0 for all k ∈ N.For k = 0, y c d0 (x) = χ [ 1 2 ,1] 2 (x)/2.Since the exact solution cannot be computed analytically, we compute its MhFE approximation on a finer mesh (512 × 512-mesh).Since the modes y c dk (x) = 0 for all even k ∈ N, it suffices to show the results for odd modes as well as for k = 0. Table 5 presents the results with truncation index N = 23, since the results regarding the efficiency indices are similar for higher modes.The computational times presented include the times for computing the approximations of the exact modes on the finer mesh.The numerical results again show the efficiency of the majorants for both, the discretization error and the cost functional.This is especially observed for the majorant related to the cost functional, which is very close to the exact value (in spite of a really complicated y d ).The majorant M ⊕ |•| exposes an overestimation but anyway provides realistic estimates of the errors in the state and control functions measured in terms of the combined error norm.

Conclusions and Outlook
In [27], the authors derived functional-type a posteriori error estimates for MhFE approximations to linear parabolic time-periodic boundary value problems.In this work, we extend this technique to the derivation of a posteriori error estimates for MhFE solutions of the corresponding distributed optimal control problem which leads to additional challenges in the analysis.The reduced optimality system is nothing but a coupled parabolic time-periodic PDE system for the state and the adjoint state.We are not only interested in computable a posteriori error bounds for the state, the adjoint state and the control, but also for the cost functional.In case of linear time-periodic parabolic constraints, the approximation via MhFE functions leads to the decoupling of computations related to different modes.Due to this feature of the MhFEM, we can in principle use different meshes for different modes and independently generate them by adaptive finite element approximations of the respective Fourier coefficients.To assure the quality of approximations constructed in this way, we need fully computable a posteriori estimates, which provide guaranteed bounds of global errors and reliable indicators of errors associated with the modes.Then, by prescribing certain bounds, we can finally filter out the Fourier coefficients, which are important for the numerical solution of the problem.This technology will lead to an adaptive multiharmonic finite element method (AMhFEM) that will provide complete adaptivity in space and time.The development and the analysis of such an AMhFEM goes beyond the scope of this paper, but will heavily be based on the results of this paper as described above.It is clear that the functional a posteriori estimates derived here for time-harmonic parabolic optimal control problems can also be obtained for distributed time-harmonic eddy current optimal control problems as studied in [19,20,21].
Example 1, we set the desired state

Table 1 .
on grids of different mesh sizes.Moreover, Table2presents the CPU The
(19)in(t) sin(x 1 π) sin(x 2 π), where T = 2π/ω with ω = 1.It is easy to see that this function is time-analytic, but not timeperiodic.As in the first example, we compute the MhFE approximation of the desired state and solve the systems(19)and (