Preconditioning of discrete state- and control-constrained optimal control convection-diffusion problems

We consider the iterative solution of algebraic systems, arising in optimal control problems constrained by a partial differential equation with additional box constraints on the state and the control variables, and sparsity imposed on the control. A nonsymmetric two-by-two block preconditioner is analysed and tested for a wide range of problem, regularization and discretization parameters. The constraint equation characterizes convection-diffusion processes.


Introduction and problem setting
Optimal control problems with constraints given by partial differential equations (PDE) arise in numerous important applications, where we pose the task to control processes, governed by PDEs.In addition to the state equation constraint, in practice it is often required that the solution obeys also other types of restrictions, such as inequality or sparsity constraints.Due to their complexity, such problems can in general only be analysed via computer simulations.During the simulation process we face the necessity to solve large scale nonlinear and/or linear systems of algebraic equations, making the use of iterative solution methods and preconditioning a must.
In this paper we analyse and illustrate numerically the performance of a particular nonsymmetric two-by-two block preconditioner, suitable for PDE-constrained optimal control problems with additional inequality constraints on the state and the control variables.
The optimal control problem we consider in this work reads as follows: find a solution (y, u) that minimizes the cost functional , (1) subject to the following constraints.
(C.1)The stationary convection-diffusion equation with boundary conditions y = g D on ∂ D , ∂ y ∂n = g N on ∂ N .(C.2) Inequality constraints for the control, referred to as box-constraints, namely, u ≤ u ≤ u, u < 0 < u (cf., e.g., [1,2]) together with the requirement that the control variable is sparse, thus, it is zero in some parts of the domain.The latter is achieved by introducing the term (T2) (the L 1 term) in (1).(C.3) Inequality constraints on the state, requiring that y ≤ y ≤ y holds in .State constraints can be handled in various ways, cf., e.g., [7,8].As in [9], we choose to use the Moreau-Yosida penalty function method, expressed in term (T3).We mention that in practice the box constraints could be imposed only on a part of the domain, 0 ⊂ , however we assume here that 0 = .
Above, y is the state variable, y d is some given desired (target) state we want to achieve, reflected in term (T0), and u is the control variable, imposed as a source term, to attain the target; b is a given vector field and 0 < ≤ 1 is a problem parameter, which determines the strength of convection related to diffusion.The space domain is assumed to be convex polyhedral.The constant α in term (T1) is the Tikhonov regularization parameter, which is positive and small.The L 1 norm constraint introduces one more regularization parameter, β, again positive and small.For details, see [3,4] and with emphasis on solution techniques for such problems, for instance, [5,6].The parameter γ , 0 < γ < 1 is yet another small regularization parameter.
The nonlinear systems of algebraic equations, arising from discrete formulation of the first order necessary conditions (9), referred to also as Karush-Kuhn-Tucker (KKT) conditions, and how to efficiently compute their solution are the main target of this work.
The area of optimal problems constrained by PDEs has been of intense research interest for decades.The published results are extensive.In this paper we refer to some articles and the references therein, particularly related to the presented work.
As optimization framework we choose the first order necessary KKT conditions, entailing large scale nonlinear and linear systems of algebraic equations to be solved.The latter triggers the necessity to use iterative methods and robust and computationally efficient preconditioners to accelerate the convergence of those when solving the arising linear systems.
Designing highly efficient solution techniques for the distributed control problem with no additional constraints on the state and the control is well studied.Optimal and near optimal preconditioning methods are available, cf., e.g., [5,[10][11][12][13][14][15][16][17].Applying box constraints on the state is discussed in, e.g., [18] and numerical methods for that formulation of the problem are considered in, e.g., [19] for both state and control constraints, and later in [9,20].Imposing sparsity on the control variable is discussed, e.g, in [3][4][5].To our knowledge, the formulation with both state and control boxconstraints combined with L 1 -regularization and preconditioners for the case when all these constraints are active have not been studied by other authors in the semismooth Newton framework.Preconditioning is considered in [6] but using another optimization framework, namely, interior point methods.The systems to solve in the latter setting are of dimension larger than those arising in the KKT framework.Initial testing of the approach, discussed in this paper for the case when the state constraint is given by the Poisson equation is presented in [21], see also [32].
The optimal control problem, where the constraint is given by the convectiondiffusion equation is considered in, for instance in [22,23], however without all the three additional constraints being active.
The paper is organized as follows.In Sect. 2 we define the discrete form of the optimization problem, corresponding to additional state and control constraints, and its compressed two-by-two block form.In Sect. 3 we define the preconditioner and analyse its properties.The interplay between problem, discretization and method parameters is discussed in Sect. 4. In Sect. 5 we present numerical examples to illustrate the performance of the method, as well as a comparison against a relevant preconditioning technique.Concluding remarks are given in Sect.6.

Discrete algebraic form of the target optimization problem
We proceed with considering the numerical solution of minimizing (1).
For simplicity, we choose to work in two space dimensions, assume that the vector field b is divergence-free, apply only Dirichlet boundary conditions and also assume that f = 0. To discretize the problem, we use the Finite Element method (FEM).
To refer to the discretized variables in vector form, we use lower-case bold letters.As in many related studies, we assume also that all variables are discretized using the same finite element spaces, thus, the mass matrices, corresponding to y, p and q (to appear in the sequel) are the same.This is utilized in order to compress the arising algebraic systems and save computations.
The discrete constraint equation becomes then K y = M u + g.Here K represents the stiffness matrix, M is the mass matrix and g contains the state boundary conditions.To be able to easily handle the box constraints for the state and the control, we work with the lumped mass matrix, thus, in the sequel, M is diagonal, as in [40], for instance.
To solve the optimization problem numerically, we construct the discrete form of the first order necessary KKT conditions.To better explain the preconditioner analysed in this work and its relation to existing efficient preconditioning techniques we consider three cases.For notational simplicity we denote the arising system matrix as A even when we reduce its dimension.Case 1: We consider (1) with terms (T0) and (T1) when the only constraint is the state equation (C.1).The optimality system in strong form reads The variable p is the solution of the dual problem, referred also as the Lagrangian multiplier to handle the state equation.We discretize (2s * ) as K y − M u = 0 where K = L + F, where L is the discretized Laplace operator − , which is symmetric (L = L T ) and F is the discretized convection term b • ∇ y which, as ∇ • b = 0, is skew-symmetric (F = −F T ).Equation (1s * ) can thus be written as M y + K T p = M y d and the first order necessary conditions in matrix-vector form read This is a well-studied case and we include it here for completeness, as well as to illustrate the gradual increase of the complexity in the arising algebraic systems.From the third equation in (4) we express the control via the Lagrange multiplier, u = 1 α p and obtain the reduced system In order to ensure that the computed state is as close as possible to the desired state y d , the regularization parameter 0 < α < 1 is supposed to be small.The inverse of it introduces a poor scaling in the matrix, which would make the solution of systems with it difficult for the iterative methods.Therefore, even before constructing a preconditioner, we weaken the poor scaling by multiplying the second equation by √ α and replacing p by − √ α p.The system to be solved then becomes Case 2: Consider next box constraints on the control of the form u ≤ u ≤ u, understood pointwise.As shown, for instance, in [1,2], the optimality system then is as in (3), but with the following additional equation to handle the box constraints, − q +max(0, q +c(u −u))+min(0, q + c(u − u)) = 0 for any constant c > 0, (6) where q is another Lagrange multiplier.As already mentioned, the presence of conditional terms in (6) results in KKT conditions that constitute a nonlinear problem with a non-smooth Jacobian.Further, as stated, relation ( 6) is understood pointwise, thus, after discretization the conditions are checked for u in every discretization point.To facilitate writing the KKT conditions in matrix form, we introduce the sets a and I (u)  i , denoting correspondingly the points, where any of the box control constraints is violated (active points) and the rest of the points, where the control satisfies the constraints (inactive points).We introduce also the notations M (u)  a and M (u) i , denoting diagonal matrices, equal to M in the points I (u)  a and zero in I (u) i for M (u) a , and vice versa for M (u)  i .With the above assumptions the linearized KKT conditions in matrixvector form, to be solved iteratively in general until the active sets stop changing, become (e.g., [24]) where f 1 = M y d , f 2 contains the Dirichlet boundary conditions for the state y and f 3 contains the conditional terms in (6), that is, it depends on the current solution and − u.Clearly, the matrix and the right-hand-side vector in (7) depend on the solution and a and f 3 are recomputed in every nonlinear iteration.The derivations below show that A can be compressed and the constant c can be eliminated.It plays a role implicitly in determining the active and inactive sets in the nonlinear process.
To include the L 1 norm constraint, relation (6) has to be changed as follows, cf., e.g., [3,21], for any constants c 1 > 0, c 2 > 0. Thus, the conditions how to classify an entry from u into I (u) i change, however, the structure of the KKT system remains the same, as in (7) with a change in f 3 , corresponding to (8).The constants c 1,2 do not appear in the compressed KKT-system yet they affect the non-linear iterations in determining the active sets.Some choices of c 1,2 have been observed to lead to non-convergence for the Newton iteration, see e.g., [21].
Case 3: Finally we add to the KKT system the box constraints for the state variable y, i.e., we now include all terms in (1).As solving the PDE-constrained problem with all additional constraints is the main focus of this work, we write out the optimality system in detail.Following [3,19] and e.g., [41], the first order necessary conditions to find the minimum of the cost functional J(y, u) in (1) in strong form read as follows: An observation shows that from (3s) we can eliminate u = 1 α ( p − q), also seen from the corresponding discrete system (11).Equation (4 s) then reads Further, for clarity, we separate the possible different cases in Equation (4 s), and what they mean for the control variable u, as follows, From (10), we clearly see why the values of c 1,2 can be chosen freely without changing the solution.The values of c 1,2 are important as we use the min/max terms in (4s) in Sect. 2 when we define the so-called active and inactive points to detect where the control constraints are violated of fulfilled and, hence, affect the behavior of the Newton iteration.
To handle the state box constraints, analogously to the control box constraints, we introduce two more sets of points -I (y) a and I (y) i , denoting correspondingly the points, where any of the state box constraints in Equation (1 s) are violated (active points) and the rest of the points, where the state satisfies the constraints (inactive points).Similarly as for the control, we introduce also the notations M where From now on we discuss the solution of the linearized system (11), for which the matrix blocks M We first scale and compress the system (11) by performing the following steps: α ( p − q) from the third equation in (11) and then multiply the second equation by √ α, the third equation by α and scale p = − √ α p.The system becomes • Observe, that the Schur complement A of A with respect to the top pivot two-bytwo block can be computed exactly and is of the form In detail, are diagonal with entries equal to one in the active and inactive points for the control u, correspondingly.In the derivation we use the fact that . It is a positive definite diagonal matrix and we can use both its square root and inverse.
To moderate the scaling due to γ −1 we introduce y = I 1 2 y y and multiply the top block-row by y we obtain the form of the system, used in the sequel, We note that for the implementation of the action of A on a vector we do not have to form K explicitly.
Thus, the main aim of this work reduces to solve systems with the matrix A in ( 14) in an efficient and robust way, ensuring fast convergence for a broad range of problem ( ), discretization (h) and method parameters (α, β, γ ).The system is nonlinear and the Jacobian is not differentiable in the classical sense.Therefore, we use the semismooth Newton method, see [2,25].Furthermore, there are two characteristics of A that make some of the established and well-studied preconditioning techniques inefficient.Namely, earlier theoretical results on preconditioning matrices of the form ( 13) are not applicable here, because (i) the matrix K is nonsymmetric due to the convectiondiffusion state equation and (ii) the block M (u) i is singular.
Remark 1 As we are dealing with convection-diffusion state constraint equation, a note on stabilizing the problem is due.As is well-known, when solving convectiondiffusion problems, some stabilization should be added, in particular when the local Peclet number is larger than one.Otherwise the quality of the computed approximate solution cannot be guaranteed.As discussed, for instance in [40], in the context of PDEconstrained optimal control problems, additional difficulties may arise.Namely, the matrix systems that we obtain when we follow the optimize-then-discretize approach would not be strongly consistent with those when following the discretize-thenoptimize approach, meaning that the solutions to the optimization problem would not necessarily satisfy all the optimality conditions.This is the case when applying the broadly used streamline upwind/Petrov Galerkin (SUPG) stabilization.Again in [40] it is shown that for another stabilization technique, the Local Projection Stabilization (LPS) method, the matrix systems obtained from 'discretize-then-optimize' and from 'optimize-then-discretize' are the same.Another argument to apply LPS instead of SUPG is that the accuracy of the obtained approximate solution when using the former is higher.However, as reported in [40] and also observed via our own numerical experiments, not reported here, the solution methods work well whether SUPG or LPS stabilization is used, and also with no stabilization at all when such an approach is plausible.Therefore in this work we do not use stabilization.Adding, for instance LPS, entails some changes in the stiffness matrix but would not affect the performance of some standard methods, such as algebraic multigrid, shown to be efficient for such types of discrete problems.For the effect of using LPS stabilization in the framework of convection-diffusion optimal control problems we refer to [26].

Preconditioning techniques
To simplify the notations further, denote y .The assumptions regarding the matrices are that M and M 1 are symmetric and positive definite (spd).We assume also that M 2 = 0, thus, the control constraints are satisfied at least in some points.The matrix K may or may not be symmetric.We consider the matrix where K = K T , and seek a preconditioner of the form where PSSB stands for Preconditioner for Square and Singular Blocks.
The preconditioner P P SSB possesses the factorization which displays the computational cost to solve a system with P P SS B , amounting to solutions with M 1 + K T and K + M 2 , one multiplication with K T and some vector updates.The computational steps to solve systems with P P SS B are given in Algorithm 1.
Algorithm 1 Solution of a system with the PSSB preconditioner Algorithm 2 Solution of a system with the PRESB 11 preconditioner Before studying the spectral properties of the generalized eigenvalue problem Az = λP P SSB z (18) we first give some details to relate P P SSB to another efficient preconditioner in earlier studies on the problem in focus, the PRESB preconditioner.

No additional constraints on y and u
Imposing no box or sparsity constraints means that M (y) a = 0 (e.g., I = M (e.g., I u a = 0) we obtain the classical distributed control problem and the resulting matrix in compressed two-by-two block form (cf., e.g., [26]) The properties below hold for any spd matrix M. A preconditioner that is shown to have significant advantages over other techniques, in terms of convergence and computational complexity is the PRESB preconditioner (cf., e.g., [20,26,27] and the references therein) It is shown that for both forms of P P RE SB all eigenvalues of P −1 P RE SB A 0 are real and belong to the interval [ 1  2 , 1].Furthermore, the PRESB preconditioner possesses convenient factorized forms, which means that the solution of a system with P P RE SB requires to solve two half-sized systems with M +K and M +K T and to perform one matrix multiplication with K .For completeness we include Algorithm 2 to solve a system In [21] the preconditioner in ( 16) is introduced, which in the case M 1 = M 2 = M reads as follows As the factorization in (22) suggests, the solution with P 0 P SSB has nearly the same computational cost as P P RE SB .We show below that the spectrum of (P 0 P SSB ) −1 A 0 is also real and has the same eigenvalue bounds as P −1 P RE SB A 0 .To this end, consider and follow the well-known approach to simplify the generalized eigenvalue problem, with μ = 1 λ − 1.From the first equation in (23) we see that for z 1 = M −1 K T z 2 we have that μ = 0 and λ = 1.For μ = 0 and z 1 = M −1 K T z 2 we obtain the relation Clearly, μ is real and positive, which shows that λ < 1.We next multiply by M − 1 2 from left and right, denote and see that μ is an eigenvalue of Utilising the fact that the eigenvalues of a matrix and its transpose are the same, μ is of the form 2a 1+a 2 .Using the assumption that K + K T is positive definite, then a is positive and the algebraic expression is bounded independently of the value of a, namely, We prove in this way the following result.
Theorem 1 Let M and K be real matrices of size n, M be symmetric positive definite, K + K T be positive definite, N(M) ∩ N(K ) = {0} and let

General case
Consider now the general case (18), where y , and analyse Here, K = K T , M 2 is singular, and λ, z 1 and z 2 may be complex.Clearly, λ is not zero.
We transform (24) to the equivalent problem where μ = 1 λ − 1.From (25) we see straightforwardly, as for P 0 P SSB , that for z 1 and z 2 such that z 1 = M −1 K T z 2 we have μ = 0, thus, λ = 1.In other words, 1 is an eigenvalue of multiplicity at least n, the size of M and K .Assume next that μ = 0, i.e., z 1 = M −1 K T z 2 and substitute in the second equation in (25), In this way we obtain another generalized eigenvalue problem with real matrices, where S is spd and Q is nonsymmetric.Thus, to analyse and bound the eigenvalues λ = 1 we analyse μ by considering the generalized eigenvalue problem Then, a straightforward computation shows that λ 2 = − λ 1 μ 2 μ 1 +1 and then (iii) If μ 2 = 0 then λ is real.The undesired case would be to have μ 1 < 0 and The following lemma is useful.
Then for the eigenvalues μ of the generalized eigenvalue problem X w = μY w, μ ∈ C, Here i denotes the complex unit, and * H -the complex conjugate of * .
In our case X = Q and Y = S.We note, that M and M 2 are diagonal and they commute.For the symmetric part of Q we obtain Analogously we obtain is the skew-symmetric part of Q.The matrix W + is diagonal with entries equal to either 1 or 2, thus it is spd.The matrix W − is also diagonal with entries equal to either 1 or 0, thus, it is symmetric positive semidefinite (spsd).To recall, I is the identity matrix, I a ∩I (u) i is not necessarily empty.In these terms, in our notations, relation (28) becomes The consequence of Lemma 1 is that the eigenvalues of Q + v = μ + Sv define lower and upper bounds for μ 1 and, analogously, the eigenvalues of To estimate the eigenvalue Re(μ), i.e., μ 1 we make use of the following corollary of Courant-Fischer's lemma.Below, unless stated otherwise, λ i (A) denotes the ith eigenvalue of a generic square matrix A.
Lemma 2 [29,Corollary 3.14] Let A and B be symmetric matrices and let λ i (A), λ i (B), λ i (A + B), λ i (AB) denote the ith eigenvalue, where we order the eigenvalues in an increasing order.Then Clearly, Q + is symmetric and all its eigenvalues are real.However, there is no guarantee that the smallest eigenvalue μ + min is positive.In fact, there may be some negative eigenvalues, which are the cause some of the eigenvalues λ to have real part larger than 1. (Fig. 8 indicates that possible loss of efficiency for P P SS B can be due to very small eigenvalues of the preconditioned system.)On the other hand, the largest eigenvalue μ + max is positive because otherwise the matrix Q + would be negative definite, which is not the case.With μ + max > 0 then condition (b) from Lemma 2 can be used to obtain an upper bound of μ + max which, in turn, is an upper bound for Re(μ).To analyse the case, we provide a lower bound of the real part λ 1 , respectively, a upper bound for μ 1 .Thus, the next task is to estimate μ + max , which is shown in Theorem 2. For clarity we keep the currently used notations.

Theorem 2 Consider the generalized eigenvalue problem
, K is a real nonsymmetric matrix, originating from discretizing the convection-diffusion operator − + b • ∇ using piece-wise linear finite elements on a regular triangular mesh with characteristic mesh-size parameter h.Further, M is diagonal positive definite (a lumped mass matrix) and M 2 = M I where C is a generic constant, not depending on α, γ , or h.

Proof
We trail the following chain of arguments.
Observe that , which is diagonal with positive entries, i.e., it is spd.Next we see that with K = W −1/2 K W 1/2 .Thus, K and K are similar and have the same spectrum.Since (32) defines a congruence transformation, we also see that Q + and Q have the same number of positive and negative eigenvalues.However, as already mentioned, there is no guarantee that all eigenvalues are positive.
Without loss of we assume that the largest eigenvalue of Q + , and, respectively, of Q is positive.From (32), using Lemma 2(b) we see that Next, we need to estimate λ max ( Q).Even though K and K are similar, we cannot use this directly in estimating the spectrum of Q.Therefore, we use Gershgorin type of arguments to obtain an upper bound of the eigenvalues of Q, based on the same type of arguments for K .Denote K = {k i j } and W = {w i }, i, j = 1, . . ., n.Then for the entries of Then, the maximum absolute row sum of K , r K can be related to that of K , r K , as follows, For linear FEM discretization on a regular triangular mesh we have where |b| = max For the spectral bounds of the other matrices involved the following eigenvalue bounds hold, Next we obtain eigenvalue estimates for S. Taking into account that diagonal matrices commute, we have where w = K T v.
As the eigenvalues of K K T are equal to the square of the singular values σ 1 ≥ . . .≥ σ n of K , then, following [30] and the references therein, to estimate the term ( * ) in (39) we can use the following Gershgorin-type of bounds, where r i (K ) = n j=1 |k i, j | and c j (K ) = n i=1 |k i, j | are the absolute row and column sums, correspondingly.Also, r i (K ) = r i (K ) − |k ii | and analogously for the other terms.As b is bounded, the bounds ξ 1 and ξ 2 decrease with decreasing .
To quantify ξ 1 we use that k ii = 4 and also and after applying Lemma 2 combined with the estimates in (34) and (35) we get If we in addition assume that > |b|h, thus, that the local Péclet number is less than one, we can simplify (42) to obtain where c, c , c , c, c, C are some generic positive constants, not depending on h, α, or γ .
We note, that estimates of the type k ii = 4 can be done for any irregular mesh, where instead k ii ≤ d , where d is the maximum number of edges per mesh-node.
Based on the estimate (42) and using the form of λ 1 in ( 27) we obtain a lower bound of it in (24), Thus, the larger the value of the factor ≡ (1 + γ )h 2 √ α γ is, the closer to zero the bound for λ 1 becomes.
To summarize, n of the eigenvalues in ( 24) are real and equal to 1.The imaginary parts of all the rest are by modulus less than 1.The lower bound of the real parts is bounded from below as in (43).The bound approaches zero if α, γ or approach 0. However, for fixed values of those parameters and h → 0 the bound approaches 1/2.Figures 3, 6 and 7 illustrate typical cases of parameter value combinations, where the lower bound is around 0.2. Figure 8 shows some specific parameter combinations, where the smallest eigenvalue of the preconditioned system becomes close to zero.The factor also indicates how the interplay between the problem and discretization parameters influence the quality of the preconditioner.This question is discussed in the next section.

Interplay between problem, discretization and method parameters
The target numerical problem involves five parameters.The problem parameter determines how singularly perturbed is the given convection-diffusion problem.The method parameters are the regularization variables α, β and γ and h is the discretization parameter.As is shown below, the convergence of the PSSB-preconditioned compressed KKT system depends explicitly on , α, γ, h and implicitly on β, which determines the number of active/inactive points for the additional state and control box-and sparsity constraints.Due to the complexity of the problem and the existing mutual parameter dependencies, a convergence study cannot be done by letting these parameters individually approach their limit values.There are several considerations to take into account.
1. Since no stabilization is included in the discretization, h should be chosen small enough with respect to .2. In [7,20] it has been shown that when we impose box constraints on the state using the Moreau-Yosida approach, keeping h fixed and decreasing γ leads to a stagnation of the achieved error.In [7] some optimal relations between h and γ have been discussed, of the type γ ≈ h 2 in 2D. 3. A certain relation between and α can also be foreseen.To illustrate this interplay, consider two mathematically equivalent forms of the state equation, We can also regard (44) (right) to be a dimensionless form of (44) (left).
Then the corresponding optimization problem becomes The second form is obtained by multiplying the right side of the constraining PDE by 1/ while leaving the objective functional unchanged.For = 1, since the second formulation does not take the scaling into account in the target functional, the solutions to the optimization problems will differ.To take the scaling into account we can modify the target functional as In (47) the minimum of relates to the minimum of J 1 as y 1 = y 2 and u 2 = 1 u 1 .That is to say, we have equivalence up to a scaling on the control for the optimization problem and given the solution to either formulation we can trivially reconstruct the other.In this sense the solutions to both formulations contain the same information and this observation is useful when predicting solver performance.The form (47) also indicates that for small values of we have an effective L 2 -stabilization parameter, defined as α e f f := 2 α, suggesting that when is small, we do not have to choose α very small to ensure good quality of the state.Clearly, scaling u in the equation entails scaling of all terms in the cost functional that contain the control, revealing further relations between the problem and the regularization parameters.These relations warrant special attention in a separate study.

Numerical illustrations
We illustrate the performance of the PSSB preconditioner for various values of the parameters h, α, β and γ .

Discretization
The optimal control test problems in this study are discretized using a triangular mesh and linear finite element basis functions.Also, as already commented in Remark 1 and in order to avoid excessive complications in the analysis of PSSB, we do not consider any stabilization in the discrete convection-diffusion problem.

Stopping criteria and solver details
The solution procedure comprises a nonlinear solver, the semismooth Newton method.The KKT system is solved until the number of the points in the active sets is minimized (ceases changing).To mitigate cases when the active set changes with only very few points that are counted active or inactive due to numerical fluctuations of the values of the current state (present in the compressed system) we consider the solution to be converged when all changes leading to a set change are less than the discretization error h 2 .
As is well known, the convergence of Newton-type of methods depends on the quality of the initial guess.The numerical experiments are performed on a sequence of nested meshes and the initial guess for the semismooth Newton method on each next finer mesh is the prolonged solution obtained on the previous coarse mesh.
In addition, the initial guess for the linear solver at each nonlinear iteration is updated to correspond to the current active sets, as described in Remark 3. The procedure on a given mesh is outlined in Algorithm 3.
Algorithm 3 Solution of the non-linear problem on a given mesh 1: Given initial guess ( y, p, q) 2: while no convergence achieved do 3: Use y to construct the active sets I − y).Use p and q to construct the active sets i and the right hand side term f u using (8) and u = 1 α ( p − q).

4:
Solve (14) using GCR to obtain ( y, p), using initial guess y = I Use active sets, the max/min conditions from Step 3 and ( y, p) to recover q 7: Check convergence After initialization, Steps 6, 7, 3 can be done in a single for-loop.8: end while Remark 3 When carrying over the current solutions from one Newton step to the next on a given mesh, a subtle but important detail is due to be mentioned.In each Newton step, after solving the linear system in (14), we obtain the solution [ y, p] T .As the relations between [ y, p] T and [ y, p] T depend on the currently active points, we recover [ y, p] T and, using the new active set, generate the initial guess for the next [ y, p[ T .These transformations (Steps 3, 4 and 5 in Algorithm 3) entail a small additional computational cost, however they result in significantly reducing the number of linear iterations, thus we gain in efficiency.For reference, the iteration counts presented in the tables in [31] correspond to the case when [ y, p] T is used in the next Newton step without the latter transformation.
The linear system with the Jacobian is solved iteratively with a relative stopping criterion chosen as P −1 (Ax − f ) / P −1 f < 10 −6 .The stopping criterion for the inner solvers with In order to avoid explicitly computing K at each non-linear step, instead of solving the systems with y , and similarly for AGMG2.The constants c 1 and c 2 in ( 8) are chosen to be c 1 = c 2 = 1/α in all tests.As a side comment we note that for very small values of α the matrix y , which is a diagonal matrix.In such cases one could possibly replace the AMG solver by some simple solver, e.g., use only the diagonal term and solve it exactly.However, this case is not discussed any further as no problems in solving the inner systems with AGMG were observed for the parameter ranges tested.
Implementation All experiments are performed on a Dell Latitude 7490 laptop with an Intel ® Core ™ i7-8650U CPU and 16 GB of RAM.The generation of the meshes and the FEM matrices is done in FEniCS ( [33]).The outer linear solver is the PSSB-preconditioned generalized Conjugate Residual (GCR) method with a maximum of 40 iterations allowed, cf.[34].As an inner solver for the block systems we use the AGMG-preconditioned GCR ( [35][36][37][38]).
The numerical tests are performed in Julia 1. 6.3 ( [39]).We provide elapsed run time for the solvers on each mesh level (i.e., excluding mesh operations and FEM assembly).
We minimize the cost functional J (y, u) as in (1), subject to (2) in two space dimensions using the test problems below.The problem domain is the unit square Problem P1 Consider the optimal control problem with We assume Dirichlet boundary conditions along the entire boundary with values defined by the desired state.
Problem P2 Let the desired state be y d = sin(2π x 1 ) sin(2π x 2 ) and the convection field be defined as b = 1 0 T .We assume homogeneous Neumann boundary condition at x 1 = 1 and Dirichlet boundary conditions on the rest of the boundary, defined by the desired state.1], the performance of the preconditioners P P RE SB and P P SSB is compared using Problem P1 without non-linear constraints, i.e., the standard distributed control problem (1).The computational steps to solve systems with P P RE SB and P P SSB are nearly the same and, as expected, both perform very similarly although the iteration counts are not identical.We note that PRESB cannot be straightforwardly generalized for systems of the form (15). Performance of the equivalent forms of the state equation (44) Consider Problem P1 and denote the vector field as b e f f = δb with b as in (48).Thus, δ = 1 corresponds to the standard formulation of the state equation as in (44)(left) and δ = 1/ corresponds to the scaled formulation (44)(right).In Table 2 in [31] we illustrate the equivalence of the formulations on pairs of parameter choices.The observation is that the iteration counts are identical across all pairs and the elapsed time for the solvers are comparable.The obtained solutions are also equivalent up to a scaling of the control.In all tests in this work we use the standard formulation.

PSSB performance with all constraints active
To illustrate the performance of the PSSB preconditioner in the case when all additional constraints are imposed, we consider Problem P1.For = 0.1 we impose control box constraints as u ∈ [5, −2.5] and state box constraints y ≤ −1.3 and β = 10 −4 .As in [21] γ = α 1/4 .For = 0.01 we impose the control box constraints as u ∈ [−2, 2] and β = 10 −3 .The iteration counts are shown in Table 1 and the control and the state are plotted in Figs. 1 and 2. In Fig. 3 we plot the eigenvalues of the preconditioned linear systems for one particular discretization.
The performance of the PSSB preconditioner without state-constraints is shown in Table 2.
Performance of PSSB with all constraints active, alternate field We consider Problem P2 with state box constraints −0.5 ≤ y ≤ 0.5 and β = 10 −3 .For = 0.1 we set the control box constraints as u ∈ [10, −10].For = 0.01 we set the control box constraints as u ∈ [3, −3].We fix γ = 0.05 for both values of .Iteration counts are listed in Table 3. Plots of the obtained control and state are shown in Figs. 4 and  5, plots of the eigenvalues of the preconditioned systems are shown in Fig. 7.

Assessment of the results of the numerical tests
The target problem, considered in this work, being a combination of a convection-diffusion equation as a general state constraint and additional constraints, imposed on the state and the control, poses significant difficulties when solving the arising nonlinear and linear algebraic systems.
For all experiments, the nonlinear iterations have been very few, which is due to several factors.On one side, the solution on a previous coarser discretization mesh provides a good initial guess for the nonlinear iteration on the next finer mesh as well as for the related (outer) linear iterations.This approach is not new and has been utilized, cf., e.g., [20] and the references therein.The addition of the constant c 2 in (8) is crucial in aiding the convergence of the Newton iteration, especially for small values of α.Another factor is the efficient PSSB preconditioner, which allows in a relatively few iterations to achieve sufficient accuracy of the linear solver for the related Jacobian system, so that the performance of the nonlinear iterations is not degraded.
As seen in the results presented in Tables 1, 2, 3, the performance of the PSSB preconditioner is very satisfactory across problem sizes, exhibiting nearly independent convergence with respect to h.The outer convergence is in general affected by the   As PSSB requires solutions of inner systems with some block matrices, we show that these can be solved very efficiently using an algebraic multigrid preconditioner.We also show that the relative stopping criterion of the inner systems can be chosen relatively large, 10 −2 , which speeds up the solution with PSSB without sacrificing the convergence of the outer iterative solution method.Cases of unfavourable spectrum of the PSSB-preconditioned matrix In Tables 1, 2, 3 we observe some increase in the number of linear iterations when α decreases, which is also aligned with the analysis of the spectrum of the preconditioned matrix.However, the increase in the iteration counts for the outer linear solver observed for the tested parameters is not significant.As we see from (43), choosing very small values of α, γ and can be unfavourable.To illustrate this effect, we consider Problem P2 with α = 10 −5 , = 10 −4 and plot the eigenvalues of the preconditioned system in Fig. 8a.
To illustrate the case with a small value of γ , we revert to = 0.1 and choose γ = 5 × 10 −5 .The eigenvalue plot is depicted in Fig. 8b.However, the choices of the parameters, in particular a very small γ disobey the mutual interplay dependencies, discussed in Sect. 4 and should thus be avoided when solving problems of practical interest.Comparisons with a Schur complement based preconditioner In this study, when solving the linear systems arising in the semismooth Newton method, we use the novel preconditioner P P SS B .Although we have no knowledge of previous research that considers the given problem formulation, there are numerous studies dealing with matrices with similar algebraic structure.For completeness we include a comparison against a preconditioner, used when solving the problem with L 1 -regularization and control box-constraints, not including box-constraints on the state.The preconditioner in question is the Schur complement based preconditioner from [5], of the form where S is an approximation of the Schur complement S = K M −1 K T + M 2 of the system.In the current setting there are a few viable choices of S, the direct extension from [5] to the form with re-scaled K -blocks being S = ( K + M 2 )M −1 ( K T + M 2 ).Numerical tests, not included here, show that for our problem the choice S = ( K + M)M −1 ( K T + M 2 ) results in faster convergence and thus we use this form in  the comparisons.Note that, applying P Schur requires a solution of systems with the inner blocks ( K + M) and ( K T + M 2 ), which is done by AGMG1/2, respectively, implemented as described for the blocks when using P P SS B .Table 4 show the iteration counts and timings when using P Schur for Problem P1 for the case = 0.01.We include iteration counts for two different stopping tolerances of the inner AGMG block solver.All other problem parameters are identical to the tests with P P SS B presented in Table 1.
We point out that, when using P Schur , the number of outer linear iterations depends quite significantly on the tolerances in the AGMG block solvers, which is especially evident for the examples with larger α.When using P P SS B , such a sensitivity is not  observed and the inner tolerance is set to 10 −2 for all included tests.For the cases where α is small, the number of outer iterations for both preconditioners is similar, but P P SS B is computationally cheaper, resulting in faster overall execution time.Notably, the iteration counts for P Schur are obtained for the scaled version of the reduced KKT system (14), which differs from the form used in [5].In [5] the nonlinear iteration counts are much larger than those reported here.The difference can be explained by a multitude of factors.For instance, we employ a mesh hierarchy which provides a good initial guess on the finer levels, utilized in both the non-linear and the linear solvers.The solutions from one non-linear iteration to the next are updated to match the currently determined active sets of points, as described in Remark 3. Furthermore, the inclusion of the constant c 2 in ( 8) is crucial and reduces significantly the number of non-linear iterations, especially when α is small.The latter improvement is independent of the use of mesh hierarchies.

Conclusions
In this work we consider the numerical solution of discrete optimal control problems, constrained by the convection-diffusion equation with three additional constraints, namely, two-sided constraints on the state and the control, and required sparsity for the control.The algebraic problems arising from the KKT system are nonlinear and large-scale, in two-by-two block form with nonsymmetric off-diagonal blocks and one of the diagonal blocks being singular.
The solution framework consists of nonlinear semismooth Newton iterations and an outer linear solver for the Jacobian, which changes at each nonlinear iteration.The nonlinearity is of a particular nature as it reflects the fact that the set of active points can differ for different mesh and problem parameters.
To solve the related linear systems iteratively we consider a particular two-by-two block preconditioner, PSSB.Since applying PSSB requires the solution of two systems with block matrices, we use an inner iterative solver, preconditioned by an efficient algebraic multigrid scheme, AGMG.
The properties and the performance of PSSB are analysed as functions of problem, discretization and three regularization parameters.The good performance of PSSB is confirmed by extensive numerical experiments.The timing results confirm a nearly linear computational complexity of the overall nonlinear solver, which includes the inner-outer solver for the corresponding Jacobian.
We note that some of the regularization parameters, namely, α and γ counteract.Due to α and h being small and 1/γ being large, the arising matrices are in general very ill-conditioned, also due to poor scaling of the blocks.To mitigate the latter effect we introduce a particular additional scaling of the system matrix.
Provided the complexity of the considered optimization problem, the solution procedure exhibits nearly optimal performance for a broad range of values of the involved parameters, that are of practical interest.
Assuming for simplicity that c 1 = 1 and renaming c 2 = c, the matrix form of the linearized KKT conditions becomes as follows,
matrices with zero elements on the diagonal for indices in the set I (y) a and I(u)  i , correspondingly, and the intersection I (y u) = I (y)

a
= 0 and I is the identity matrix.The parameters α and γ are positive constants in the interval (0, 1).Assume that b is bounded.Then, for the maximum eigenvalue μ + max of the above generalized eigenvalue problem, the following upper bound holds true,

l b 2
x,l + b 2 y,l , b x,l and b y,l are the x and y coordinates of b at node l, l = 1, . . ., n 2 .Thus,

i
and the right hand side term f y = γ −1 (I (y) + y + I (y)

Fig. 8
Fig. 8 Problem P2 with state constraints: Eigenvalues of the preconditioned system, α = 10 −5 n all eigenvalues of the generalized eigenvalue problem A 0 z = λP 0 P SSB z.Then all λ i are real and lie in the interval [ 1 2 , 1].We note another relation between P 22 P RE SB and P 0 P SSB , namely, there holds

Table 1
Problem P1: All constraints active

Table 4
Problem P1, = 0.01: P Schur , All constraints active.(P P SS B iterations for the same problem in the right column of