Non-commutative Discretize-then-Optimize Algorithms for Elliptic PDE-Constrained Optimal Control Problems 1footnote 11footnote 1This project was supported in part by the grant NSF-DMS 1522672 of the United States..

# Non-commutative Discretize-then-Optimize Algorithms for Elliptic PDE-Constrained Optimal Control Problems 111This project was supported in part by the grant NSF-DMS 1522672 of the United States..

Jun Liu Zhu Wang Department of Mathematics and Statistical Sciences, Jackson State University, Jackson, MS 39217, USA. Department of Mathematics, University of South Carolina, Columbia, SC 29208, USA.
###### Abstract

In this paper, we analyze the convergence of several optimize-then-discretize and discretize-then-optimize algorithms, based on either a second-order or a fourth-order finite difference discretization, for solving elliptic PDE-constrained optimization or optimal control problems. To ensure the convergence of a discretize-then-optimize algorithm, one well-accepted criterion is to choose or redesign the discretization scheme such that the resultant discretize-then-optimize algorithm commutes with the corresponding optimize-then-discretize algorithm. In other words, both types of algorithms would give rise to exactly the same discrete optimality system. However, such an approach is not trivial. In this work, by investigating a distributed control problem governed by an elliptic equation, we first show that enforcing such a stringent condition of commutative property is only sufficient but not necessary for achieving the desired convergence. We then introduce some suitable semi-norm penalty/regularization terms to recover the lost convergence due to the inconsistency caused by the loss of commutativity. Numerical experiments are carried out to verify our theoretical analysis and also validate the effectiveness of our proposed regularization techniques.

###### keywords:
PDE-constrained optimization, elliptic optimal control, discretize-then-optimize, optimize-then-discretize, finite difference method, trapezoidal rule, Simpson rule, regularization
journal: Computational Optimization and Applications

## 1 Introduction

In the past few decades, PDE-Constrained optimization and optimal control problems Lions1971 (); Hinze2009 (); Troltzsch2010 (); De_los_Reyes_2015 () have gained many efforts from the scientific computing community, due to the increasingly broad applications and tremendous computational challenges. Generally speaking, there are two different pathways of constructing feasible numerical algorithms: (i) optimize-then-discretize (OD) approach and (ii) discretize-then-optimize (DO) approach gunzburger2003perspectives (); Herzog_2010 (); Hinze_2011 (). In OD approach, one first derives the first-order necessary continuous optimality conditions analytically, and then discretizes the continuous optimality system with appropriate discretization schemes, such as finite difference method strikwerda2007finite (); LeVeque2007 () and finite element method hackbusch1992elliptic (); Johnson2009 (). After that, the resultant discretized linear/nonlinear systems can be solved by various well-designed efficient iterative solvers Saad2003 (); Briggs2000 (); Trottenberg2001 (); Borzi2009 (); hackbusch2013multi (). While in DO approach, one first discretizes the original continuous problem directly to obtain a fully discretized large-scale finite-dimensional optimization problem, which is then solved by any existing numerical optimization algorithms Betts2001 (); Nocedal2006 () that can handle a large-scale size of decision variables. Intuitively, the DO approach is more straightforward in practice since the optimization processes can be automatically operated by using black-box optimization solvers Betts2001 (); Betts2005 (); GilMS05 (); snopt76 (). Moreover, the DO approach provides more flexibility in handling additional complex constraints, bounds, and regularization terms Kameswaran2008 ().

Although it has been widely used in practice with very satisfactory outcomes, the above described OD approach indeed has at least the following two widely known drawbacks gunzburger2003perspectives (). First, the continuous KKT PDE system (3) may be difficult to derive for more general constraints, such as complicated nonlinear state gradient constraints. This difficulty may limit its broader applications. Second, the discretized KKT linear system (4) may become non-symmetric for more general discretization schemes, such as a standard 9-point high-order finite difference scheme (as shown later). Without symmetry, the system (4) will not correspond to the optimality system of any discrete optimization problem, which hence yields inconsistent gradients of the objective functional that lead to troubles in using gradient-based optimization algorithms for solving nonlinear problems. To overcome the above disadvantages of the OD approach, one may resort to the DO approach, which ensures the generation of a symmetric discretized KKT linear system that is easier to be solved computationally. However, one major disadvantage of the DO approach is that its control approximations may be spuriously oscillating. Therefore, in this paper, we attempt to systematically address this problem through developing new regularization techniques. The proposed regularization methods guarantee the convergence of non-commutative DO discretization schemes, which are especially attractive to those applications that can not be treated easily by the OD approach and often demands advanced functional analysis skills from the practitioners.

It is worthwhile to mention that such spurious numerical instabilities have been observed, but not thoroughly resolved, across various disciplines and application scenarios. In 1979, Huntley Huntley1979 () observed “spurious oscillations in the optimal control given by the Riccati algorithm are shown to stem from the use of Simpson’s rule to approximate the integrals” from their numerical results without giving further mathematical justification. Similar checkerboards numerical instabilities were thoroughly discussed and investigated in topology optimization Jog1996 (); Sigmund1998 (). In optimal control of ODEs, similar numerical oscillations were noticed due to inappropriate discretization within the context of the pseudo-spectral (PS) methods Fahroo2008 (); Garg_2010 () as well as Runge-Kutta methods Hager_2000 (); DontchevHagerVeliov00 (). As for optimal control of PDEs, similar numerical oscillations has been carefully investigated in Heinkenschloss2010 (); Leykekhman2012 () when applying discontinuous Galerkin methods to advection-diffusion PDE constrained optimal control problems with mild interior and boundary layers. However, there is no very satisfactory solution to this challenging problem in general. A currently well-accepted strategy is to redesign the discretization scheme (such as adding stabilization terms Becker2007 ()) such that discretization and optimization commute in the sense that OD and DO approaches are essentially the same Apel2012 (). However, such a strategy has a very limited applicability in dealing with high-order discretizations, adaptive meshes as well as more complicated constraints, since it becomes impractical to adjust the corresponding discretization scheme in such a dedicated manner.

On the other hand, for better efficiency and more accuracy in solving PDEs, significant progress has been made in developing high-order (compact) finite difference discretization schemes Lele_1992 (); Shu_2003 (); Shu_2009 (). One would agree that PDE-constrained optimization problems can benefit much more from a successful high-order discretization scheme, since it dramatically reduces the number of decision variables to optimize and hence the overall computational cost. Nevertheless, high-order discretization schemes have been less frequently discussed or utilized in the emerging field of PDE-constrained optimization, possibly due to low regularity of the solution with active control/state constraints. Bearing this regularity limitation in mind, we believe high-order discretization schemes are still very desirable because there are many applications turn out to be of high regularity or having certain singularities but can be treated with suitable high-order schemes and numerical techniques. For example, Wachsmuth and Wurst Wachsmuth2016 () recently applied finite element method to elliptic boundary control problems with point-wise control constraints and established its exponential order of accuracy with special mesh refinement strategies. Rösch and Wachsmuth R_sch_2017 () developed a higher-order finite element discretization based on a new mass lumping strategy for a control constrained elliptic optimal control problem, which achieves up to fourth-order accuracy on locally refined meshes. In the framework of OD approach, the first high-order finite difference discretization scheme for solving elliptic optimal control problem was proposed by Borzì Borzi2007 (), who applied the extended nine-point approximation to the Laplacian that achieves a formal fourth-order accuracy in the case of no control constraints. As mentioned in Borzi2007 () without proof, that finite difference scheme, indeed, numerically attains a fourth-order accuracy even in the presence of active control constraints.

In this work, we investigate the convergence properties of the DO algorithms with either a second-order or a fourth-order finite difference discretization scheme in alignment with the corresponding OD algorithms, and further propose an innovative regularization technique to ensure the anticipated convergence of the discussed DO algorithms (see Tables 1 and 2 for a quick overview). To the best of our knowledge, fourth-order finite difference discretization schemes within the framework of DO approach have not been discussed in the literature, which makes the convergence analysis provided in this paper a major contribution of this work. Considering the convergence properties of finite difference discretization schemes have been rarely discussed Borz2005 () in the field of numerical optimal control of PDEs, our convergence analysis is also of independent interests to the large PDE constrained control community.

The rest of this paper is organized as follows. In Section 2, we present two optimize-them-discretize algorithms based second-order and fourth-order finite difference schemes, respectively, for solving a 2D prototype elliptic optimal control problem. In Section 3, we describe and analyze several different discretize-then-optimize algorithms for solving the prototype elliptic optimal control problem, with second-order or fourth-order finite difference discretization for the state equation, and trapezoidal and Simpson’s rule for the objective functional, respectively. A simple 2D numerical example demonstrates the possible convergence failure when using Simpson’s rule. In Section 4, we propose to resolve the convergence issue of Simpson’s rule by adding semi-norm regularization terms to the discrete objective functional, which is inspired by the carefully comparing the discrete KKT system arising from OD and DO approaches. In Section 5, a few numerical experiments are conducted to validate our theoretical conclusions and illustrate the effectiveness of our proposed algorithms. Finally, some concluding remarks are drawn in Section 6.

## 2 Optimize-then-Discretize algorithms with finite difference discretizations

As a motivating example, we consider the following 2D elliptic optimal control problem Lions1971 (); Hinze2009 (); Troltzsch2010 () of

 minu∈UJ(z,u) =12∫Ω(z−g)2dx+α2∫Ωu2dx, (1)

subject to

 (2)

where the source term , desired state , is a regularization parameter, and is the Laplacian operator. If , the above optimal control problem admits a unique optimal solution Lions1971 (), which satisfies the following first-order necessary optimality KKT system

 (KKT)⎧⎪⎨⎪⎩−Δz−u=finΩ,z=0on∂Ω,−Δp+z=ginΩ,p=0on∂Ω,αu−p=0inΩ, (3)

where is called adjoint state. Due to the strict convexity of (1-2), the above optimality KKT system (3) is also sufficient. In other words, the solution to the KKT system (3) gives the unique optimal solution to (1-2). Notice that the derivation of the KKT system (3) using either calculus of variation or Lagrange functional techniques corresponds to the ‘optimize’ step in the OD approach, while the ‘discretize’ step is to discretize the KKT system (3) with an appropriate discretization. Such an OD approach will lead to a large-scale sparse linear system that often requires efficient iterative solvers, and its resulting discrete approximate solutions are expected to converge to the optimal solutions under certain regularity assumptions.

Given a step size , the 2D space domain is partitioned uniformly by grid points with , where and . Denote the set of all interior grid points of the partitioned space domain by . In OD approach, upon discretizing the Laplacian with a second-order 5-point central finite difference scheme, we achieve the following fully discretized KKT linear system

 (OD-2)⎧⎨⎩−Δhzh−uh=fh,−Δhph+zh=gh,αuh−ph=0, (4)

where denotes the discrete Laplacian matrix corresponding to the stencil

 h−2⎡⎢⎣0101−41010⎤⎥⎦

and , , , , are the vectorized discrete approximation of the associated variables on . Here we assume the homogeneous Dirichlet boundary conditions are already enforced through the discrete Laplacian matrix LeVeque2007 (). As another major contribution of this paper, we are also particularly interested in studying a fourth-order discretization scheme and the corresponding discretized KKT system. In OD approach, by discretizing the KKT system (3) with a fourth-order compact 9-point central finite difference scheme hackbusch1992elliptic (), we arrived at the following discretized KKT linear system

 (OD-4)⎧⎨⎩Fhzh−Rhuh=Rhfh,Fhph+Rhzh=Rhgh,αuh−ph=0, (5)

where the stencil

denotes the discrete negative Laplacian matrix and the weights of averaging the right-hand-side, respectively.

For any two grid functions and defined on , we define the discrete weighted Euclidean inner product

 (v,w)=h2N−1∑i=1N−1∑j=1vijwij.

Based on this inner product, we can define the corresponding induced discrete vector norm

We also define the standard infinity vector norm

 ∥v∥∞=max1≤i,j≤N−1|vij|.

We first present the following lemma that will be used in proving our convergence results. Both the following inequalities are the special cases of the discrete Sobolev embedding inequalities Bramble_1966 (); brenner2007mathematical (), but their proofs within the context of finite difference methods are not widely available in the literature, to our best knowledge. For completeness, we provide our elementary proofs adapted from the classic book (Samarskii2001, , page 281).

###### Lemma 2.1.

For any grid function defined on that vanishes on the boundary nodes ( or ), there exists a generic positive constant , independent of , such that there hold

1.  ∥v∥∞≤C∥Δhv∥,
2.  ∥v∥∞≤C∥Fhv∥.
###### Proof.

See Appendix A for the detailed proof. ∎

Based on (3) and (4), one can prove the above scheme (4) has a formal second-order accuracy under suitable regularity assumptions, as stated in the following theorem.

###### Theorem 2.1.

Suppose the exact solution triple of (3) satisfies , then the finite difference scheme (OD-2) or (4) has a second-order accuracy with respect to the infinity norm .

###### Proof.

See Appendix B for the detailed proof. ∎

Similarly, the formal fourth-order accuracy of the scheme (5) is proved in the following theorem.

###### Theorem 2.2.

Suppose the exact solution triple of (3) satisfies , then the finite difference scheme (OD-4) or (5) has a fourth-order accuracy with respect to the infinity norm .

###### Proof.

See Appendix C for the detailed proof. ∎

###### Remark 2.1.

In the above both Theorem 2.1 and Theorem 2.2, the assumptions on the regularity of solution triple are rather restrictive, considering the fact that and in most realistic applications. In other words, the conclusions hold only when the data and and the domain are sufficiently smooth. Nevertheless, in this paper we focus on investigating the subtle convergence properties of discretize-then-optimize algorithms with currently available finite difference discretizations, and our chosen test problems indeed have very smooth solutions. We will give a numerical example when the assumptions are not met.

It is possible to slightly weaken the current regularity assumptions while retaining the same order of accuracy (possibly in different norms), but with technically more involved cell-average of the right-hand-side and integral representations (instead of Taylor series expansions); we refer to (hackbusch1992elliptic, , p. 68) and (Jovanovic2014, , p. 127) for further discussion. However, we will not go into further detail on this issue, since it is beyond the scope of this paper. To minimize the needed regularity assumptions (also depends on the smoothness of domain), one may choose to use finite element discretizations, which will be left as our future work.

###### Remark 2.2.

We also point out that the above convergence proofs are mainly based on summation by parts, the discrete analogue of integration by parts. In some cases, it is more convenient to show the discretization scheme is both consistent and stable under certain norm, which will automatically imply convergence via the Lax-Richtmyer equivalence theorem. The consistency is relatively easier to verify through analyzing the truncation errors, and the stability can often be proved by estimating the eigenvalues/singular values of the underlying coefficient matrix, as shown in the following.

We first rewrite the approximation errors system (32) into

 Lh[epez]:=[Ih/αΔh−ΔhIh][epez]=[−FG]. (6)

Let denotes the Hermitian part of , i.e.,

 H(Lh)=[Ih/α00Ih].

Let and denote the smallest singular value and eigenvalue, respectively. Then it follows from a well-known singular value inequality (Horn1994, , p. 151)

 σmin(Lh)≥λmin(H(Lh))=min(1/α,1)

that (notice that )

 ∥L−1h∥2=1σmin(Lh)≤1λmin(H(Lh))=max(α,1), (7)

which proves the stability of the scheme (4) under the spectral norm (induced by the Euclidean norm).

## 3 Discretize-then-Optimize algorithms with finite difference discretizations

We now briefly describe how the OD approach works for the same optimal control problem (1-2). With the same uniform mesh and the given homogeneous boundary conditions, the objective functional (1) can be approximated by the trapezoidal rule as

 Jh(zh,uh) =12(zh−gh)T(zh−gh)+α2uThuh, (8)

which has a formal second-order accuracy. Note that, to better match the DO approach, we omit the scaling factor of appearing in the numerical quadrature for simplicity of notations, which will not affect our analysis in the following discussion. Similar to the above OD approach, we discretize the state equation (2) by the second-order accurate 5-point finite difference scheme and get

 −Δhzh−uh=fh. (9)

This is the ‘discretize’ step of the DO approach. Obviously, the obtained discretization (8-9) gives rise to a large-scale finite-dimensional linearly constrained optimization problem, which can be then solved by any existing optimization algorithms that work efficiently with a large-scale size of decision variables. By the strict convexity of (8-9), the obtained numerical approximations in fact approximately solve the corresponding first-order necessary optimality KKT system regardless of the used optimization algorithms. The discrete optimality KKT system can be obtained by forming a Lagrange functional (through introducing the discrete Lagrange multiplier or adjoint state )

 L(zh,uh,ph) =12(zh−gh)T(zh−gh)+α2uThuh+pTh(−Δhzh−uh−fh) (10)

and then setting its gradient to be zero, which gives

 (DO-2-Trap)⎧⎪⎨⎪⎩Lph=−Δhzh−uh−fh=0,Lzh=−Δhph+zh−gh=0,Luh=αuh−ph=0. (11)

Coincidentally, this KKT system (11) from the DO approach is identical to the above KKT system (4) obtained from the OD approach. In this case, the DO approach and OD approach are said to be commutative and the discretization scheme has the optimal convergence properties. However, we will demonstrate in the following that such a ‘coincidence’ does not happen in general. Indeed, we only simply replace the trapezoidal rule used in approximating the objective functional (1) with the more accurate, fourth-order Simpson’s rule, which gives

 Jh(zh,uh) =12(zh−gh)TQh(zh−gh)+α2uThQhuh, (12)

where denotes a positive diagonal matrix corresponding to the quadrature weights of the 2D compound Simpson rule davis2014methods (). With the common Kronecker product notation , can be explicitly expressed to be

 Qh=19diag([4,2,4,2,⋯]⊗[4,2,4,2,⋯]).

Using the new discretized functional (12) and the same second-order discretization (9) of the state equation, we can derive the corresponding discrete optimality KKT system

 (DO-2-Simp)⎧⎨⎩−Δhzh−uh=fh,−Δhph+Qhzh=Qhgh,αQhuh−ph=0. (13)

Clearly, (13) is identical to (11) if equals to an identity matrix. But, the matrix is not an identity matrix when the Simpson’s rule is used. The question then arises: will the approximate solution of (13) converge to the exact solution of (3) in this case?

In the following, we will illustrate what would happen in numerical simulations by considering a simple 2D example. Let and choose such that the exact optimal solution reads

 z(x,y)=sin(πx)sin(πy),u(x,y)=sin(2πx)sin(2πy)/α.

In Figure 1, we show both surface and color-map of the computed optimal control with the trapezoidal and Simpson’s rules, respectively. Very different from the trapezoidal rule case (11), we observe that the Simpson’s rule case (13) gives a badly behaved approximation (with spurious checkerboard oscillations). We performed a mesh refinement test on both approaches and list the convergence results in Table 1.

Recall that (12) is derived from the Simpson’s rule, which yields more accurate approximation of the objective functional than (8) that uses the trapezoidal rule. Hence, after combining with the state equation discretized by the central finite difference scheme, we would reasonably expect the computed solution approximation of (13) has at least the same level of accuracy as that of (11). Unfortunately, such a reasonable expectation fails to hold and the computed solution of (13) does not converge to the exact solution of the original continuous problem as shown in Table 1. From another perspective, this failure is not completely surprising since the discretize KKT system (13) based on the Simpson’s rule dramatically deviates from the discretize KKT system (4) given by the OD approach. Indeed, we can reformulate (13) into the following form (aligned with the system (4))

 ⎧⎪⎨⎪⎩−Δhzh−uh=fh,−Δhph+zh=gh+(Qh−Ih)gh+(Ih−Qh)zh,αuh−ph=α(Ih−Qh)uh, (14)

where the extra residual terms (highlighted in bold) on the right hand side can be shown to be of since . In other words, the discrete system (13) or (14) is inconsistent to the continuous KKT system (3) and hence the obtained numerical approximation will not converge to the exact solution of (3). However, in the following we will consider a high-order discretization scheme for the state equation, in which we numerically show that such a consistency between the discrete and continuous KKT systems is only a sufficient condition but not a necessary one.

Next, we use a fourth-order 9-point finite difference scheme for the discretization of the state equation, and consider both the trapezoidal rule and Simpson’s rule for the objective functional approximation, respectively. Note that the former gives the following discretized constrained optimization

 (15)

 (16)

Then, by forming the Lagrange functional and setting its gradient to be zero, we obtain the corresponding discretized KKT system

 (DO-4-Trap)⎧⎨⎩Fhzh−Rhuh=Rhfh,Fhph+zh=gh,αuh−Rhph=0, (17)

for the trapezoidal rule and

 (DO-4-Simp)⎧⎨⎩Fhzh−Rhuh=Rhfh,Fhph+Qhzh=Qhgh,αQhuh−Rhph=0, (18)

for the Simpson’s rule. Clearly, both the KKT systems (17) and (18) based on the DO approach are obviously different from the KKT system (5) obtained from the OD approach with the same 9-point discretization scheme for the state equation. However, according to our numerical tests on the 2D elliptic control problem listed in the Table 1, the system (17) based on the trapezoidal rule gives convergent approximation with a fourth-order accuracy, while the system (18) based on the Simpson’s rule fails to converge. It illustrates that the commutativity between the DO and OD approaches is not a necessary condition. To intuitively understand the observed fourth-order accuracy of the system (17), we can reformulate it into the following form (in view of matching (5))

 (19)

It is straightforward to check that and under mild regularity assumptions, which hence indicates that (17) differs from (5) only with some perturbation. Therefore, we expect the system (17) to deliver at least overall second-order accuracy, given the system (5) has a fourth-order accuracy. In fact, our numerical results and convergence analysis show that the obtained approximation and in (17) actually has a fourth-order accuracy, as given in the following theorem.

###### Theorem 3.1.

Suppose the exact solution triple of (3) satisfies , then the finite difference scheme (DO-4-Trap) or (17) has a fourth-order accuracy in and and a second-order accuracy in with respect to the infinity norm .

###### Proof.

See Appendix D for the detailed proof. ∎

In Table 1, we summarized the above discussed four different schemes within the framework of DO approach and their commutative properties with the OD approach, the observed convergence or divergence based on our numerical tests. A key observation is that the commutativity between the DO and OD approach seems to be only a sufficient but not necessary condition for an employed discretization scheme to become convergent within the DO approach, e.g., (17) is convergent.

The next section is devoted to better understanding and efficiently resolving the observed possible convergence failure of the DO approach with the Simpson’s rule. The central question we would like to address is how to achieve the expected convergence of the DO approach, even when the underlying discretization scheme does not guarantee the commutative property of the ‘optimize’ and ‘discretize’ processes. We propose to add well-chosen regularization penalty terms to the objective functional, which turns out to be work quite well in improving the convergence of the DO approach, which is shown to be very effective in removing the spurious oscillations based on our numerical experiments.

###### Remark 3.1.

With the similar arguments as in Remark 2.2, we can also prove the stability of the scheme (13), which, however, is not consistent to the continuous KKT optimality system (3). In other words, the scheme (13) is not convergent. In fact, the reduced coefficient matrix of the system (13) has the form

 ^Lh=[Q−1h/αΔh−ΔhQh], (20)

which leads to the stability result

 ∥^L−1h∥2≤1λmin(H(Lh))=max(16α9,94).

## 4 Regularized Discretize-then-Optimize algorithms

As shown in Figure 1, the direct use of Simpson’s rule leads to spurious checkerboard oscillations in the computed control approximations . To suppress such spurious oscillations, one of the natural approaches is to penalize the undesirable non-smoothness of the computed control approximations , which can be straightforwardly implemented through adding to the original discrete objective functional a semi-norm regularization term such as

 γ∥∇huh∥22=γ(∇huh,∇huh)=γ(uh,−Δhuh)=γuTh(−Δh)uh, (21)

where denotes the discrete gradient defined by forward finite difference. However, based on our numerical experiments, adding such a regularization term has very limited effects on suppressing the observed spurious oscillations, and its overall performance is highly sensitive to the manually chosen regularization parameter . In particular, the optimal order of accuracy cannot be obtained even by tuning carefully, although we acknowledge that such a regularization term indeed promotes smoother approximations . In Figure 2, we plot the computed control approximations from adding the regularization term (21) with different values of . Among them, the case gives smooth approximation of the control variable but is of an incorrect magnitude, while the case gives the most satisfactory control approximation but it consists of spurious oscillations. In all our tested cases, we did not observe any uniform convergence with a fixed order of accuracy as the mesh size is refined. It seems that adding such a regularization term only mildly alleviates but not completely eliminates the spurious oscillations, which, of course, does not lead to satisfactory approximations, as illustrated in Figure 2.

Based on the above discussion, better regularization techniques are needed in order to eliminate the oscillations and achieve the optimal order of accuracy. To this end, we turn our attention to scrutinize the difference between the discretized KKT systems (5) and (18). By comparing (5) and (18), we conclude that they become identical if is replaced by . Recall that comes from the used Simpson’s quadrature rule, which is a diagonal matrix. What is the mathematical implication behind replacing by ? One possible interpretation is to decompose into the sum of two different stencils

 Rh=112⎡⎢⎣010181010⎤⎥⎦=⎡⎢⎣000010000⎤⎥⎦+112⎡⎢⎣0101−41010⎤⎥⎦,

where the first stencil corresponds to the trapezoidal rule and the second one is a scaled discrete Laplacian associated to a negative semi-norm term. Thus, by subtracting a scaled semi-norm term from , we would achieve a matrix close to . Inspired by this key observation, we propose a regularized optimal control strategy to improve the non-communicative DO-2-Simp and DO-4-Simp schemes by introducing the following objective function:

 ^Jh=12(zh−gh)T(Qh−γΔh)(zh−gh)+α2uTh(Qh−γΔh)uh. (22)

This is equivalent to adding two semi-norm regularization terms and , i.e.,

 ^Jh=12(zh−gh)TQh(zh−gh)+α2uThQhuh+12(zh−gh)T(−γΔh)(zh−gh)+α2uTh(−γΔh)uh, (23)

which will give rise to the following discrete KKT system (corresponding to the fourth-order scheme)

 (DO-4-Simp-Reg)⎧⎪⎨⎪⎩Fhzh−Rhuh=Rhfh,Fhph+(Qh−γΔh)zh=(Qh−γΔh)gh,α(Qh−γΔh)uh−Rhph=0. (24)

Obviously, this new regularized KKT system is very different from (5), since for any . Based on our numerical investigation (also see the following Example 2), the optimal choice of the free parameter is . In this case, the regularized KKT system (24) achieves an optimal fourth-order accuracy in control approximations.

Applying the same regularization technique to the second-order scheme as in (13), we can similarly get

 (25)

By taking , our numerical experiments show that the regularized KKT system (25) produces an expected second-order accuracy in control approximations.

Though originally motivated by the observed numerical failure of Simpson’s rule, we also tested the proposed regularization terms with the case of trapezoidal quadrature rule, e.g., (11) and (17). More specifically, we consider the following two discretized KKT systems

 (26)

and

 (DO-4-Trap-Reg)⎧⎪⎨⎪⎩Fhzh−Rhuh=Rhfh,Fhph+(Ih−γΔh)zh=(Ih−γΔh)gh,α(Ih−γΔh)uh−Rhph=0, (27)

associated to the second-order scheme () and fourth-order scheme (), respectively. Our numerical experiments show that the proposed regularized schemes are also able to output the optimal control approximations. In this sense, our proposed -regularized DO approach is robust for using in general settings.

In summary, we report in Table 2 the convergence results of the above discussed schemes. Compared to Table 1, we have successfully achieved the expected convergence and optimal order of accuracy for all schemes, without enforcing the stringent condition of commutativity between the OD and DO approaches.

It is also worthwhile to point out that the optimal value of parameter in our introduced regularization terms have been found: and for the fourth-order scheme and second-order scheme, respectively, which makes the proposed regularization approach to be essentially parameter-free. Although we have tried to justify our choice of such regularization terms in the above discussion, a complete understanding of the observed convergence in numerical results requires more analysis, which is left as an open problem at the moment. Next, we will illustrate the above algorithms in several numerical experiments.

## 5 Numerical experiments

In this section, we provide numerical examples to demonstrate the effectiveness of our proposed methods. All simulations are implemented using MATLAB R2016a on a laptop PC and the optimality KKT systems are solved by MATLAB’s built-in backslash sparse direct solver. All approximation errors are measured in the infinity norm in comparison with the known exact solution (if available).

Example 1. We first consider a simple 1D example. Let and choose such that the exact optimal solution reads

 z(x)=sin(πx),u(x)=sin(2πx)/α.

In Figure 3, we plot the computed control by using the second-order scheme for state equation with trapezoidal and Simpson’s rule for objective functional approximation, respectively. In comparison to the trapezoidal rule (left plot), we observe strong spurious oscillations in the obtained approximation by Simpson’s rule (right plot). The computed control by the second-order scheme with regularized trapezoidal and Simpson’s rule are plotted in Figure 4, respectively. With our proposed regularization term, the spurious oscillations of Simpson’s rule have been eliminated and the obtained approximations by both trapezoidal and Simpson’s rule are indistinguishable. Furthermore, the added regularization does not degrade the approximation accuracy of the trapezoidal rule. The case with the fourth-order scheme gives very similar results and we hence choose not to duplicate the plots for simplicity.

Example 2. Let and choose such that the exact optimal solution reads

 z(x,y)=sin(πx)sin(πy),u(x,y)=sin(2πx)sin(2πy)/α.

In Figure 5, we plot the computed control by the second-order scheme with regularized trapezoidal and Simpson’s rule. Compared to Figure 1 that is associated to the control without regularizations, the spurious oscillations caused by Simpson’s rule have been clearly eliminated. To validate our theoretical analysis, we report in Table 3 the approximation errors and expected second-order of accuracy of computed optimal control with both trapezoidal and Simpson’s rules. Notice the approximation errors by Simpson’s rule without regularization seem to be of and do not converge at all. Similar results from the fourth-order scheme are shown in Figure 6 and reported in Table 4, where the expected fourth-order accuracy of Simpson’s rule is observed after adding the regularization term. In all the cases, introducing the proposed regularization term eliminates the spurious oscillations, while leading to the expected accuracy of numerical schemes. We also point out that the observed fourth-order accuracy obtained by the fourth-order scheme (17) with the trapezoidal rule relies on the coincidental commutative property between and , which may not be valid with other discretization schemes.

Example 3. In this example, we illustrate the performance of our proposed algorithm with the case of a non-attainable discontinuous target function. Let , and choose

 g(x,y)={sin(πx)sin(πy),if |x−12|≥14 and |y−12|≥14,0,otherwise.

Notice the lower regularity of the target function will bring some difficulties to any standard finite difference schemes, since they usually require higher regularity of the solutions for obtaining the formal order of accuracy in infinity norm. To avoid confusing or misleading readers, we choose not to report the approximation errors and order of accuracy as in Example 2 because the exact optimal solution is not available. Instead, we compare our second-order and fourth-order algorithms by visualizing the obtained approximations in Figure 7 and Figure 8, respectively. Clearly, the obtained approximations of optimal control with regularization are free of the numerical oscillations observed in the case without regularization.

In Figure 9, we plot the target state and computed optimal state with and , respectively. As we expected, a smaller leads to better tracking property in approaching the optimal state to the target state , but at the cost of solving a discrete KKT system with a larger condition number.

## 6 Conclusions

Discretize-then-optimize algorithms are widely used in numerically solving PDE-constrained optimization and optimal control due to their relatively straightforward implementation by maximizing the reuse of existing numerical optimization algorithms. However, the subtle interaction between discretization and optimization complicates the rigorous convergence analysis of those non-commutative discretize-then-optimize algorithms, especially for the cases with high-order discretizations. By working with a prototype linear-quadratic elliptic PDE optimal control problem, we compared and analyzed the convergence of both optimize-then-discretize and discretize-then-optimize algorithms with several second-order and fourth-order finite different discretizations. Several new theoretical conclusions on the convergence of both OD and DO algorithms are obtained with some elementary error analysis techniques.

The DO algorithms may lead to different discrete KKT systems from the OD algorithms, but they may still be convergent provided the obtained discretization remains stable and consistent. To resolve the observed problematic numerical oscillations when approximating the objective functional with Simpson’s rule, we proposed to penalize the non-smoothness of both terms in the objective functional by adding a well-chosen semi-norm regularization term. Numerical experiments show the effectiveness of our proposed regularization techniques in eliminating spurious numerical oscillations. Our future work includes the application of our proposed regularization techniques to cases with control and/or state constraints Casas_1986 (); Bergounioux1999 (); Borz2005 () and with other popular discretization schemes, such as finite element method and discontinuous Galerkin methods.

## Acknowledgments

The authors would like to thank Dr William W. Hager for his revision suggestions that lead to better readability and Dr Buyang Li for his constructive comments on our proofs. The first author also would like to thank Drs Qiang Du and Xiaochuan Tian for their insightful discussion and helpful suggestions, which was made possible by the support from the NSF/CBMS Conference“Nonlocal Dynamics: Theory, Computation and Applications”, held at Illinois Institute of Technology, Chicago, during June 4-9, 2017.

## Appendix A Proof of Lemma 2.1: two discrete Sobolev embedding inequalities.

###### Proof.

(1) In matrix Kronecker product notation, we explicitly have

 −Δh=(I⊗A)+(A⊗I)

where is a identity matrix and is the 1D discrete Laplacian tridiagonal matrix

 A =1h2⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣2−10⋯00−12−10⋯000⋱⋱⋱000⋯−12−100⋯0−12⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦∈R(N−1)×(N−1).

It is well-known that the eigenvalues of are given by

 λk=4h2sin2(kπh2),k=1,2,⋯,(N−1),

which satisfy the following obvious inequalities (provided )

 4k2=k2π2(2π)2≤k2π2sin2(kπh2)(kπh2)2=λk≤4h2,

where we have used the fact that the function is monotonically decreasing on with

 1≥sin(θ)θ≥sin(π/2)π/2=2π.

By the spectral properties of the Kronecker sum, the eigenvalues of are given by

 λk,l=4h2(sin2(kπh2)+sin2(lπh2)),k,l=1,2,⋯,(N−1)

with the following similar inequalities (provided )

 4(k2+l2)≤λk,l≤8h2.

Let be the normalized unit eigenvector of corresponding to each eigenvalue , and the set of all eigenvectors forms a complete orthonormal basis of . In fact, the -th component of the eigenvector has the following explicit expression

 wk,li,j=2sin(ikπh)sin(jlπh).

Given any grid function defined on that vanishes on the boundary nodes, there is a unique representation

 v=N−1∑k=1N−1∑l=1bk,lwk,l

and there holds

 (−Δh)v=N−1∑k=1N−1∑l=1bk,l(−Δh)wk,l=N−1∑k=1N−1∑l=1bk,lλk,lwk,l,

 ∥Δhv∥2=∥(−Δh)v∥2=N−1∑k=1N−1∑l=1b2k,lλ2k,l.

From , for each -th component, we have and hence

 |vi,j|=∣∣ ∣∣N−1∑k=1N−1∑l=1bk,lwk,li,j∣∣ ∣∣≤N−1∑k=1N−1∑l=1|bk,lwk,li,j|≤2N−1∑k=1N−1∑l=1|bk,l|. (28)

Applying the Cauchy-Schwarz inequality for the above sum, it gives

 |vi,j|2 ≤4(N−1∑k=1N−1∑l=1|bk,l|)2=4(N−1∑k=1N−1∑l=1|bk,lλk,l|1|λk,l|)2 (29) ≤4(N−1∑k=1N−1∑l=1b2k,lλ2k,l)(N−1∑k=1N−1∑l=11λ2k,l)=4∥Δhv∥2(N−1∑k=1N−1∑l=11λ2k,l). (30)

From the above estimates of eigenvalues , we can obtain

 1λ2k,l≤116(k2+l2)2,

 N−1∑k=1N−1∑l=11λ2k,l≤N−1∑k=1N−1∑l=1116(k2+l2)2≤116∞∑k=1∞∑l=11(k2+l2)2=:^C<∞,

where the last series can be easily shown to be convergent by using either comparison test or integral test. Using this estimate in (30) it follows that

 |vi,j|2≤4^C∥Δhv∥2,

which further implies

 ∥v∥∞=max1≤i,j≤N−1|vi,j|≤C∥Δhv∥.

This proves the part (1) of the conclusion .

(2) In matrix Kronecker product notation, we can explicitly write the discrete Laplacian matrix as

 Fh=((I⊗A)+(A⊗I))−h26(A⊗A).

Again, by the spectral properties of the Kronecker product, the eigenvalues of are given by

 μk,l=43h2(3sin2(kπh2)+3sin2(lπh2)−2sin2(kπh2)sin2(lπh2)),k,l=1,2,⋯,(N−1).

Moreover, by making use of properties of the Kronecker Product and the spectral decomposition of , it is straightforward to show that has exactly the same eigenvectors as . Applying the inequality for any and the fact we get

 2sin2(kπh2)sin2(lπh2)≤2sin(kπh2)sin(lπh2)≤sin2(kπh2)+sin2(lπh2),

which gives

 8(k2+l2)3≤83h2(sin2(kπh2)+sin2(lπh2))≤μk,l≤λk,l≤8h2.

Following the same proof arguments in part (1), we can obtain (with replaced by and different constants)

 ∥v∥∞=max1≤i,j≤N−1|vi,j|≤C∥Fhv∥,

which proves the part (2) of the conclusion. ∎

## Appendix B Proof of Theorem 2.1: the second-order accuracy of the scheme (4).

###### Proof.

Under the given regularity assumptions, the exact solution triple of (3) satisfies the system

 ⎧⎨⎩−Δhz−u=fh+F,−Δhp+z=gh+G,αu−p=0, (31)

with the truncation errors terms and for some generic positive constant .

Let , , and denote the approximation errors on all mesh grid points. Then the difference between (31) and (4) gives

 ⎧⎪⎨⎪⎩−Δhez−eu=F,−Δhep+ez=G,αeu−ep=0, (32)

The discrete inner product of the first equation in (32) and yields

 (33)

Similarly, the discrete inner product of the second equation in (32) and yields

 (34)

Also, the discrete inner product of the third equation in (32) and yields

 α(∇heu,∇hez)=(∇hep,∇hez). (35)

By using (35), the sum of (33) and (34) leads to

 α∥Δhez∥2+∥Δhep∥2=α(−F,Δhez)+(−G,Δhep)≤α∥F∥∥Δhez∥+∥G∥∥Δhep∥≤αCh2∥Δhez∥+Ch2∥Δhep∥≤α(C2h4+14∥Δhez∥2)+(C2h4+14∥Δhep∥2)=(α+1)C2h4+α4∥Δhez∥2+14∥Δhep∥2. (36)

By combining the like terms, we further obtain

 (37)

Due to the discrete Sobolev embedding inequalities (see Lemma 2.1)

 ∥ez∥∞≤C∥Δhez∥and∥ep∥∞≤C∥Δhep∥,

we get

 ∥ez∥∞≤Ch2, ∥ep∥∞≤Ch2, and ∥eu∥∞=1α∥ep∥∞≤Ch2. (38)

This concludes the second-order accuracy in the infinity norm of the finite difference scheme (4) . ∎

## Appendix C Proof of Theorem 2.2: the fourth-order accuracy of the scheme (5).

###### Proof.

Under the given regularity assumptions, the exact solution triple of (3) satisfies the system

 ⎧⎨⎩Fhz−Rhu=Rhfh+H,Fhp+Rhz=Rhgh+S,αu−p=0, (39)

with the truncation errors terms and for some generic positive constant .

Let , , and denote the approximation errors. Then the difference between (39) and (5) gives

 ⎧⎪⎨⎪⎩Fhez−Rheu=H,F