New preconditioners for semi-linear pde-constrained optimal control in annular geometries

Optimization problems that are constrained by partial differential equations (PDEs) often pose significant computational challenges to black-box optimization algorithms. This is particular the case for non-linear problems that typically rely on multi-query solution of large-scale linearized subsystems. In this regard, this paper proposes a customized solution strategy for a class of semi-linear PDE-constrained optimization problems in annular domains. At its core, the strategy relies on a collection of new Poisson-like preconditioners that are based on boundary-adapted spectral Galerkin methods. The preconditioners are matrix-free and scale linearly with the problem size. To establish proof-of-concept, a case study solves a benchmark control problem for various model parameters. The results show that the preconditioners lead to fast solution strategies that outperform conventional direct approaches.


Introduction
Large-scale optimization problems that are constrained by partial differential equations (PDEs) play a key role in various fields of science and engineering [2,10].As a challenge, the size and complexity of the PDE-constraints presents severe computational difficulties that often prevent the use of general-purpose black-box optimizers.As a consequence, cost efficient, specialized solvers become essential [1,3,7,8].As a contribution in this direction, this paper demonstrates how to extend seminal ideas of Shen [15][16][17] to construct fast and memory-efficient optimizers for the class of semi-linear PDE-constrained optimization problems with non-linear reaction kinetics min y, u∈U ad 1 2 Ω (y(x) − y d (x)) 2 dx + ρ 2 Ω u(x) 2 dx, (1a) The paper focuses on the specific cases of either homogeneous (1) Dirichlet or (2) Neumann boundary conditions, where Ω ⊂ R 2 is an annular domain of the type For a given non-linear reaction term, G(•), and Tikhonov regularization parameter, ρ > 0, the control problem (1) aims to determine the optimal state and control variables, (y * , u * ), that minimize the objective (1a).Here the optimal solution must belong to the set of feasible pairs, (y, u), that satisfy the PDE-constraints (1b) and the additional admissibility condition, u ∈ U ad .To be concrete, this paper focuses on the case of bi-lateral point-wise control constraints Point-wise bounds of the type (3) appear in a number of practical applications, where the control must satisfy, e.g., operational limitations that are not naturally captured by the underlying PDE (1b).In the limiting case, where u a := −∞ and u b := ∞, the admissible set becomes U ad = L 2 (Ω d ).This corresponds to the case where the PDE (1b) constitutes the only constraint.

Main Contributions and Outline
This paper contributes to a recent series of efforts by the authors that seek to construct fast, iterative solvers for a range of PDE-constrained optimization problems by exploiting the properties of customized spectral bases [4][5][6].This series of work aims to introduce a high-order alternative to the widely-used constellation of low-order finite-element methods and Schur-complement preconditioners that currently predominates the literature on PDE control [12][13][14].Previous efforts have mainly considered distributed control of elliptic and parabolic non-linear diffusionreaction systems.The main focus has been on problems in rectangular domains, where PDEs constitute the only constraints.As a natural extension, this paper investigates how to modify the existing methods to account for (1) bound constraints of the type (3) and (2) different geometries.For the sake of brevity, the paper restricts attention to annular domains (2).However, with slight modifications, the approach generalizes to cylindrical geometries of the type As the main contribution, this work proposes a collection of Poisson-like preconditioners that are customized for efficient solution of the control problems (1) by a semi-smooth Newton (SSN) strategy [9].Similar to a traditional Newton method, the SSN scheme solves (1) iteratively by finding a locally optimal solution to the non-linear Karuhn-Kush-Tucker (KKT) optimality conditions by solving a sequence of linearized, variable-coefficient subproblems.Direct solution of the subproblems is often time consuming and requires considerable memory-allocation.To this end, the new preconditioners are designed to promote efficient solution of the SSN subproblems by appropriate Krylov subspace (KSP) methods.Following seminal ideas of Shen [16], the preconditioners rely on fast direct solvers for constantcoefficient problems that exploit (1) the structure of boundary-adapted spectral bases and (2) the separable nature of annular domains.As the main feature, inversion of the preconditioners decouples to form to a sequence of independent 2×2 systems.This implies that the preconditioners can be applied matrix-free and scale linearly with the problem size.In addition, the independence of the 2 × 2 systems makes the preconditioners amenable to parallelization.To establish proof-of-concept, a numerical case study solves (1), where G(•) is given by a cubic non-linearity.
The results demonstrate computational efficiency and show that the preconditioners respond well to different problem sizes, boundary conditions, point-wise bound constraints and various choices of the regularization parameter, ρ > 0.
To establish the necessary background, Sect. 2 outlines how to solve the optimal control problem (1) using the SSN scheme.Further, to motivate the contributions of this paper, the section discusses some of the computational challenges that arise from discretization of the associated linearized subproblems.These challenges naturally leads to the construction of the new Poisson-like preconditioners in Sect.3. Section 4 presents numerical results, while Sect. 5 draws overall conclusions and addresses future work.

Motivation: A Semi-smooth Newton Method
This paper solves the control problem (1) by a semi-smooth Newton strategy [9].The SSN scheme seeks to generate a locally optimal solution, (y, u), by solving the first-order necessary optimality system Here the boundary conditions of the original problem (1) are preserved, G y denotes the Fréchet derivative of G with respect to the state variable, y, and the optimal control satisfies u = H (p) = max(u a , min(ρ −1 p(x), u b )).In the special case In the concrete case of annular domains (2), the system (5) can be recast to polar coordinates.To this end, define the functions where The optimality system then reads where To solve the KKT conditions (7), the SSN scheme considers the system as an operator equation F (y, p) = 0 and solves it by generating a recursive sequence of iterates, x i := (Y i , P i ), 1 ≤ i ≤ k, where the next iterate, x k+1 := (Y, P ), is found by solution of the linearized optimality conditions: Here where H p denotes the generalized Newton derivative of H with respect to the adjoint variable, P , i.e.,

Numerical Challenges: Discretization of the SSN Subproblems
As a numerical challenge, the SSN scheme relies on successive solution of coupled PDEs in the form (8). Upon discretization, this leads to repeated solution of large saddle-point problems.To illustrate the associated difficulties, consider a spectral-Galerkin discretization of the linear subproblems (8).To this end, define the boundary-adapted approximation spaces Let K := N • M and define S K := V N × F M .The discrete Galerkin approximation of (8) then seeks to find Y, P ∈ S K such that where v, w := and P N,M , consider the truncated series expansions where l(k) := k + M 2 .Now, define the (N − 1) × (N − 1) matrices associated with the basis {ψ k } N−2 k=0 : Note that appropriate choices of the basis functions {ψ k } N−2 k=0 ∈ V n will be constructed in Sect.3. Further, let Γ and Ξ denote the M × M diagonal matrices defined by γ mn = e in (•) , e im(•) = 2πδ mn , ξ mn = mn e in (•) , e im(•) = 2πnmδ mn , (16) where δ mn denotes the Kronecker delta.Finally, consider the (MN × 1) vectors The discretized linear subproblem (8) can then be written in matrix form where Here the matrices M C , = 1, 2, 3 are defined by the elements where i, j satisfy that

New Poisson-Like Preconditioners
As a significant challenge to the numerical solution of ( 7), the SSN scheme relies on repeated solution of saddle-point problems (21) of dimension 2(N − 1)M × 2(N − 1)M.Consequently, direct solution strategies often become computational intractable.As a cost efficient alternative, the following introduces new preconditioners that seek to accelerate the inner SSN subproblems (8) by using appropriate Krylov subspace methods to solve the associated preconditioned linear systems Concretely, this paper proposes approximative constraint preconditioners of the type Following ideas of traditional Poisson preconditioners, the new preconditioners are constructed by approximating each block of the SSN subproblem (21) by the matrices, B and M c , = 0, 1, 2, that come from a spectral Galerkin discretization of the corresponding constant-coefficient problem that determines Y, P ∈ S K such that where To be efficient, the new preconditioners crucially rely on carefully chosen basis functions {ψ k } N−2 k=0 for the discrete approximation space, V N (11).To this end, this paper uses Fourier-like (FL) bases that were originally introduced by Shen and Wang in the context of traditional initial-boundary-value problems [17].As a key property to construction of the preconditioners, the FL bases lead to diagonal massand stiffness matrices, i.e., The FL bases can be constructed as part of an offline preprocessing stage in two steps: 1. Let {L k (•)} N k=0 be the Legendre polynomials.Then there exists a unique set of coefficients {a k , b k } N−2 k=0 such that Furthermore, the mass matrix, M A = ( φ j , φ i ) ij , is penta-diagonal and symmetric positive definite, whereas the stiffness matrix, S A = ( ∂ x φ j , ∂ x φ i ) ij , becomes diagonal [15].In the concrete cases of Dirichlet and Neumann boundary conditions, the coefficients, {a k , b k } N−2 k=0 are given by respectively 2. The second step computes the diagonalization Λ = Q T M A Q, where Q = (q ij ) denotes the matrix of eigenvectors and {λ i } N−2 i=1 are the associated eigenvalues.Using the matrix Q, the FL basis can be constructed by the linear combinations: (29)

Efficient Inversion of the Preconditioners
As the main feature of the preconditioners, P k , the following describes an efficient inversion procedure that exploits the orthogonal structures of the FL bases (27).To this end, consider the following preconditioning problem that is solved during each iteration of the KSP method: Note that (30) corresponds to the discrete first-order necessary optimality conditions associated with the constant-coefficient optimal control problem (26).Hence, by definition ( 22), it follows that where S ij = ( ∂ t ψ j , ∂ t ψ i ) ij and M ij = ( ψ j , ψ i ) ij .Further, by the orthogonal properties of the Fourier bases ( 16), the matrices, Γ and Ξ , are diagonal.Therefore, using the notation, it follows that the preconditioning problem (30) can be written as M independent linear systems where In addition, the properties of the FL basis, {ψ k } N−2 k=0 , implies that S and M become diagonal (29).Hence, the system (32) reduces to M(N − 1) independent 2 × 2 linear systems in the form where σ lm := C A + (C B k(l) 2 + 2πC 0 )λ m .By (33), it follows that the original preconditioning problem (30) decouples into (N − 1)M independent 2 × 2 subsystems.As a consequence, the Poisson-like preconditioners (25) scale linearly with the problem size and can be applied matrix-free.

Numerical Results
To investigate the potential of the Poisson-like preconditioners, the following case study solves the control problem (1), where the reaction term is given by the cubic non-linearity G(y) := y 3 .The corresponding problem serves as a recurring example in the control literature [18].In this case study, the goal is to track the desired state of the type where a ≤ α < β ≤ b.The following example uses the parameters, Z = 4, a = 30, α = 40 and β = b = 60.The main purpose of the study is to investigate efficiency and robustness of the preconditioners (25).To this end, the study solves (1) for different choices of (1) problem size, (2) boundary conditions, (3) regularization parameter, and (4) point-wise bound constraints of the type (3). 1s a benchmark reference, the results are compared to MATLABs state-of-the-art direct solver.All computations are carried out in [11] on a 2.9 GHz Intel processor.
The SSN scheme is said to have converged when the 2-norm difference between successive iterates is below η = 10 −4 .The KSP iterations are performed using the MATLAB function GMRES with a tolerance of = 10 −9 .The direct solver relies on MATLABs backslash command.Table 1 lists the results, where KSP iter denotes the average number of KSP iterations required for each SSN step.Note also that DOF denotes the number of degrees of freedom for each individual SSN subproblem.Hence, the total degrees of freedom, DOF T , is therefore given by #SSN steps × DOF.The results reflect some overall tendencies that generalize to other choices of the parameters, Z, a, α, β and b.Firstly, the preconditioners provide significant reductions in CPU-time compared to the direct strategy.In particular, the results show that the non-linear control problem with up to DOF T = 875,000 unknowns can be solved in less than a minute using modest hardware.Secondly, the preconditioners prove robust with respect to the problem size and the choice of boundary conditions.Thirdly, as a drawback, the number of SSN steps and KSP iterations increase as the point-wise bounds become more strict.The authors suspect that these increases in SSN steps and KSP iterations are caused by the combination of a decrease in regularity of the solution and an increase in nonlinearity of the KKT system (Fig. 1).Inspired by [16], the new preconditioners exploit the orthogonal properties of customized, boundary-adapted spectral bases.This leads to matrix-free preconditioners that scale linearly with the problem size.Numerical results have demonstrated that the preconditioners lead to fast solution of large-scale optimization problems with significant computational benefits compared to MATLABs state-of-the-art direct methods.Furthermore, the preconditioners have proven to be robust with respect to the problem size for both homogeneous Dirichlet and Neumann boundary conditions.As a challenge, numerical experiments indicated that the non-linearity of the problem increases as the point-wise bound constraints become more strict.In turn, this leads to an increase in the number of SSN steps and KSP iterations that are required to reach convergence.A future study seeks to improve this situation by providing the SSN scheme with an educated starting guess that uses a coarse-grid solution to a similar control problem with less restrictive constraints.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/),which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
To represent the approximate solutions, Y N,M

Fig. 1
Fig. 1 The computed states for (1) Dirichlet boundary conditions, (2) Neumann boundary conditions and (3) the desired state for u a = −35, u b = 35, ρ = 10 −4 .Note that both solutions manage to approximate the desired state well, despite of the bound constraints