Table Of Content

Optimization Methods for Dirichlet Control Problems Mariano Mateos 7 Dpto. de Matema´ticas, Universidad de Oviedo, Campus de Gijón, 33203 Gijón, Spain 1 0 2 ARTICLE HISTORY n Compiled January 27, 2017 a J ABSTRACT 6 We discuss several optimization procedures to solve finite element approximations 2 of linear-quadratic Dirichlet optimal control problems governed by an elliptic par- tial differential equation posed on a 2D or 3D Lipschitz domain. The control is ] discretized explicitely using continuous piecewise linear approximations. Uncon- C strained, control-constrained, state-constrained and control-and-state constrained O problems are analyzed. A preconditioned conjugate method for a reduced prob- leminthecontrolvariableisproposedtosolvetheunconstrainedproblem,whereas . h semismoothNewtonmethodsarediscussedforthesolutionofconstrainedproblems. t StateconstraintsaretreatedviaaMoreau-Yosidapenalization.Convergenceisstud- a iedfor boththecontinuousproblemsandthefinitedimensional approximations.In m thefinitedimensionalcase,weareabletoshowconvergenceoftheoptimizationpro- [ cedures even in the absence of Tikhonov regularization parameter. Computational aspects are also treated and several numerical examples are included to illustrate 1 thetheoretical results. v 9 KEYWORDS 1 Dirichlet optimal control, discretization, constrained optimization, preconditionate 6 conjugate gradient, semismooth Newton methods 7 0 . AMS CLASSIFICATION 1 49J20, 49M05, 49M15, 65N30 0 7 1 : v Xi 1. Introduction r a OnlyinthelasttenyearsithasbeenpossibletodevelopasystematicstudyofDirichlet optimalcontrolproblemsgovernedbyellipticequations,whichstartedwiththeseminal paper [1]. In that work a 2D control-constrained problem governed by a semilinear equation posed on a convex polygonal domain is studied. Several other works have been published about numerical error estimates; see [2] or [3, Ch 2.1] for a variational approach to control-constrained problems posed on smooth 2D or 3D domains, [4] for a superconvergence result on the state approximation for unconstrained 2D problems, [5] for control-and-state constrained problems posed on convex polygonal domains. In [6] the authors study the control-constrained problem in a 2D smooth convex domain taking into account the problems derived by the boundary approximation and in [7] an apparent paradox between [2] and [6] is explained. The regularity of the solution in Email:[email protected]. 1 Optimization for Dirichlet Control Problems 2 possiblynonconvex polygonal planedomainsis studied in[8];see alsotheintroduction of that paper for further references about related problems. In the recent publication [9], error estimates for a Dirichlet control problem governed by a parabolic equation are obtained. In this work, the spatial discretization of the control is studied in both the cases of continuous piecewise linear finite elements and variational approach. Just before that, several works dealing with efficient optimization methods for control or state constrained problems had appeared; in [10–16] the semismooth Newton method is thoroughly studied. Since all these papers about optimization are previous to the papers about the numerical analysis of Dirichlet control problems, very little or no reference is made in them to their applicability to the problems we are going to deal with. Only in [16] two different Dirichlet control problems with more regular solutions than the ones we treat here are studied in the infinite-dimensional case. Let us also mention that in [17] the Dirichlet boundary condition is replaced by a Robin penalization. For this kind of penalization the methods developed in the aforementioned references are directly applicable. In [18], the semismooth Newton method is studied inthe context of thevariational approach of the control. Although theauthors only exemplify their results through distributed and Robin boundary control, a com- bination of their results and some of the results we present in Section 4 can be applied to Dirichlet control. See Remark 4 below. InthisworkwedescribeoptimizationmethodsforDirichletcontrolproblemsinboth the infinite and finite dimensional cases. Convergence proofs, examples and practical implementation details are discussed for all the algorithms through the paper. In Section 2 we state the problem and prove that it is well posed in Lipschitz domains;seeLemma2.1.Sofar,onlysmooth,convex,polygonalorpolyhedraldomains had been studied. Next, we discretize the problem. As is usual in control problems, we have three choices to discretize the control: the first option is not to discretize the control, using a variational discretization as introduced in [19] for distributed problems; as a second option we may use piecewise constant functions; and finally we can use continuous piecewise linear functions. ThechoiceisnottrivialbecausefirstorderoptimalityconditionsforDirichletcontrol problems involve the normal derivative of the adjoint state. If we discretize directly the optimality system using e.g. continuous piecewise linear finite elements, we would obtain two different approximations of the control: the trace of the discrete state would be continuous piecewise linear and the normal derivative of the discrete adjoint state would be piecewise constant. In [20, Ch. 3.2.7.3] and in [3, Example 2.1.11] this difficulty is solved using a mixed formulation of the state equation which is discretized with the lowest order Raviart-Thomas element. In [21] a convergence proof for this kind of discretization is given. In this way, both the trace of the discrete state and the normal derivative of the discrete adjoint state are piecewise constantfunctions,sotheidentification ofbothwiththediscretecontrolismeaningful. The discretization of the control in the presence of control constraints is carried out in these papers using a variational discretization. This kind of approximation may be convenient whenthegradient of thestate is thevariable of interest, sincethis quantity is treated as one of the variables of the problem. Another approach, in finite dimension, is to use the variational discrete normal derivative of the discrete adjoint state introduced in [1], which is a continuous piecewise linear function. Doing so, both the trace of the discrete state and the normal derivative of the discrete adjoint state are continuous piecewise linear functions, so the identification of both with the discrete control is meaningful. Following this idea, Optimization for Dirichlet Control Problems 3 in [2] the control is not discretized for the control-constrained case, but a variational approach is followed. In this work we investigate optimization methods for the case of discretizing the control explicitely using continuous piecewise linear functions. The end of Section 2 is devoted to describe with some detail some computational aspects that will be important in the rest of the work. WenextprovideinSection3anefficientmethodtosolvetheunconstrainedproblem. We propose in Subsection 3.2 a preconditioned conjugate gradient (pcg in the rest of the work) method for a reduced problem in the spirit of [22, Section 5]. We are able to prove convergence of the conjugate gradient even for the case where the Tikhonov regularization parameter ν vanishes. In sections 4, 5 and 6 we study the convergence of the semismooth Newton method fortheconstrainedproblemsandwritepracticalalgorithmsforthesolutionofthefinite dimensional approximations. In Section 4 we deal with the control-constrained problem. In many of the aforementioned references about the semismooth Newton method forcontrol-constrainedproblemstheauthorsstudyconvergenceofanabstractproblem oraninfinitedimensionalproblem.FordistributedandNeumanncontrolproblemsthe results are immediately applicable to the finite dimensional approximations because the controls can be discretized using piecewise constant functions, and therefore the variational inequality arising from first order necessary optimality conditions can be written in an element-wise form. The same idea applies when we are dealing with a variational approach as proposed in [2, 3, 19] or a mixed formulation, as studied in [3, 20, 21]. When we use continuous piecewise linear elements,the variational inequality cannot be written in a point-wise or elementwise form; see (36d) and Remark 5. We include the analysis of Newton methods for both the original problem and the discretized one. In Section 5 we study the state-constrained problem using a Moreau-Yosida penalization. Since there are no control constraints, the analysis of the semismooth Newton method for the infinite dimensional problem is applicable to the finite dimensional one, so we do not need to repeat it. We prove that the finite-element approximation of the penalized problems converge to the solution of the penalized problem. This result cannot be deduced straightforward from the ones in the literature since the penalized functional is not of class C2. A continuation strategy as proposed in [16] is developed. Finally, in Section 6 we discuss the problem with both control and state constraints. It is well known that the main difficulty with Dirichlet control problems is the low regularity of the solutions. This regularity and related error estimates are, so far, well established in 2D polygonal domains [1, 5, 8] and 2D or 3D smooth domains [2] but there is not, up to our best knowledge, a general study in 3D polyhedral domains. Although the main focus of this work is on optimization methods, we also study the regularity of the solution and error estimates of the approximations in one example case in a polyhedron; see Example 3.6. 2. Statement of the problem and first properties Let Ω Rd, d = 2 or d = 3, be a bounded domain with Lipschitz boundary Γ and in ⊂ this domain consider a target state y L2(Ω). Consider also the continuous linear Ω ∈ operator S : L2(Γ) L2(Ω) such that y = Su if and only if y is the solution in the → Optimization for Dirichlet Control Problems 4 transposition sense (see Definition 2.2 below) of ∆y = 0 in Ω, y = u on Γ. (1) − Let ω Ω bea domainsuch that ω¯ Ω anddefinetwo pairs of functionsα, β C(Γ) ⊂ ⊂ ∈ and a, b C(ω¯) such that α(x) < β(x) for all x Γ and a(x) < b(x) for all x ω¯. ∈ ∈ ∈ For some fixed regularization parameter ν 0, define ≥ 1 ν J(u) = Su y 2 + u 2 2k − ΩkL2(Ω) 2k kL2(Γ) and consider the sets U = u L2(Γ) : α(x) u(x) β(x) for a.e. x Γ , α,β { ∈ ≤ ≤ ∈ } and K = y L2(Ω) C(ω¯) : a(x) y(x) β(x) for a.e. x ω¯ . a,b { ∈ ∩ ≤ ≤ ∈ } In this work we will study optimization procedures for the following four Dirichlet control problems: (PU) min J(u), (PC) min J(u), (PS) min J(u), (PCS) min J(u), u∈L2(Γ) u∈Uα,β Su∈Ka,b u∈U α,β Su∈K a,b namely the unconstrained, control-constrained, state-constrained and control-and- state constrained problems. Remark 1. Almost all the literature related to these problems is written using the state equation (1). There would be no problem in taking into account an equation of the form d Ay := ∂ (a ∂ y)+a y = F in Ω, y = u+G on Γ i i,j j 0 − Xi=1 with regular enough F, G, a 0 and a = a satysfing an uniform ellipticity 0 i,j j,i ≥ condition. Letusstate precisely whatwemean bysolution inthetranspositionsense.Consider the space Φ = φ : φ H1(Ω) and ∆φ L2(Ω) . { ∈ 0 ∈ } This is a Banach space with the graph norm φ = φ + ∆φ . Further, Φ H1(Ω) L2(Ω) k k k k k k the functions in this space satisfy ∂ φ L2(Γ). This is known to be true for smooth n ∈ domains;convexdomains,see[17,LemmaA.1];planepolygonaldomains,see[8,Corol- lary 2.3]; or polyhedral domains, see [23, Theorem 2.6.7] and the usual trace theorem. We have not been able to find a proof of this fact for general Lipschitz domains. In [24, Theorem B.2] the regularity φ H3/2(Ω) is proved. Nevertheless, as the authors ∈ Optimization for Dirichlet Control Problems 5 notice in page 165 of this reference, the trace theorem is not valid neither in H3/2(Ω) nor in H1/2(Ω), so we cannot deduce immediately that ∂ φ L2(Γ). The results in n ∈ [24,25] imply that the usualtrace result can beextended to harmonic functions in the limit cases. We show next how to take advantage of this to prove that ∂ φ L2(Γ) n ∈ in Lipschitz domains. Regarding the analysis of semismooth Newton methods –see Lemma 4.1 below– we also prove Lq(Γ) regularity for some q > 2. Lemma 2.1. Let Ω Rd, d = 2 or d = 3, be a bounded domain with Lipschitz ⊂ boundary Γ and consider φ Φ. Then, there exists q > 2 depending on the domain 0 ∈ such that, for 2 q < q , we have ∂ φ Lq(Γ) and 0 n ≤ ∈ ∂ φ C ∆φ . (2) n Lq(Γ) L2(Ω) k k ≤ k k If, further, Ω is smooth or convex or polygonal or polyhedral, then there exists t > 0 such that ∂ φ Ht(Γ) and n ∈ ∂ φ C ∆φ . (3) n Ht(Γ) L2(Ω) k k ≤ k k Proof. Denote z = ∆φ, extend z by zero to Rd, consider the Newtonian potential − centered at the origin N(x) and definew = z N, theconvolution productof z and N. ∗ Then w H2(Ω) and w H1(Ω)d, so it is clear that Tr( w) H1/2(Γ)d ֒ Lq(Γ)d ∈ ∇ ∈ ∇ ∈ → for all q < if d = 2 and all q 4 if d = 3, since the dimension of Γ is d 1. This ∞ ≤ − implies that: (a) ∂ w = w n Lq(Γ) because the boundary Γ is Lipschitz and therefore it has n ∇ · ∈ a unit normal vector defined almost everywhere; and there exists C > 0 such that ∂ w C ∆φ ; (4) n Lq(Γ) L2(Ω) k k ≤ k k (b) g = Tr(w) W1,q(Γ) dueprecisely to the definition of W1,q(Γ), and there exists ∈ C > 0 such that g C ∆φ . (5) W1,q(Γ) L2(Ω) k k ≤ k k Define now v H1(Ω) the unique solution of ∆v = 0 in Ω, v = g on Γ. Using ∈ − [24, Theorem 5.6], we have that there exists q > 2 such that if 2 q < q , then 0 0 ≤ the nontangential maximal function of the gradient of v satisfies M( v) Lq(Γ) and ∇ ∈ there exists C > 0 such that M( v) C g . (6) Lq(Γ) W1,q(Γ) k ∇ k ≤ k k As is pointed out in [26, p. 438], this implies that v has nontangential limit a.e. on ∇ Γ and we can define the normal derivative of v at a point s Γ as the nontangential ∈ limit as x s of v(x) n(s). For a precise definition of nontangential limit and → ∇ · nontangential maximal function see, e.g., the introduction of the work [26]. This, together with inequalities (6) and (5) imply that ∂ v Lq(Ω) and n ∈ ∂ v C ∆φ . (7) n Lq(Γ) L2(Ω) k k ≤ k k So we have that φ =v w and we can define in a natural way ∂ φ= ∂ v ∂ w n n n − − ∈ Lq(Γ). The estimate (2) follows from (4) and (7). Optimization for Dirichlet Control Problems 6 Forsmooth,convex,polygonalorpolyhedraldomains,thesecondresultfollowsfrom theregularity φ H3/2+t(Ω)andtheusualtracetheoremfor φ.See[17,LemmaA.1] ∈ ∇ for convex domains, [8, Corollary 2.3] for plane polygonal domains and [23, Theorem 2.6.7] for polyhedral domains. Definition 2.2. We will say that y L2(Ω) is the solution in the transposition sense ∈ of (1) if (y, ∆φ) = (u,∂ φ) for all φ Φ. (8) Ω n Γ − − ∈ Here and in the rest of the work (, ) stands for the standard inner product in X · · L2(X). The adjoint operator of S is S∗ : L2(Ω) L2(Γ) defined by S∗z = ∂ φ, where φ n → − is the unique weak solution of ∆φ= z in Ω, φ= 0 on Γ. (9) − We can write now 1 ν 1 J(u) = (Su y ,Su y ) + (u,u) = (S∗Su+νu,u) (S∗y ,u) +c Ω Ω Ω Γ Γ Ω Γ Ω 2 − − 2 2 − where c = 0.5 y 2 is a constant. Ω k ΩkL2(Ω) Fromthisexpression,wecaneasily computethederivativeofJ atapointu L2(Γ) ∈ in the direction v L2(Γ): ∈ J′(u)v = (S∗Su+νu,v) (S∗y ,v) . Γ Ω Γ − For later use, we will define now for every u L2(Γ), y = Su H1/2(Γ) and u ∈ ∈ ϕ H1(Ω) the unique solution of u ∈ 0 ∆ϕ = y y in Ω, ϕ = 0 on Γ. u u Ω u − − Discretization. To discretize the problems. consider a regular family of triangulations of Ω¯. To h h {T } simplifythenotation,wewillsupposethatΓispolygonalorpolyhedral.Relatedtothe mesh,wewillcallN thenumberofnodesand x N thenodesofthemesh.Wedefine { j}j=1 the sets of interior indexes, boundary indexes and indexes in ω¯ as I = j : x Ω , j { ∈ } B = j : x Γ and J = j : x ω¯ . For the discretization of the state and the j j { ∈ } { ∈ } adjoint state we use the space of linear finite elements Y H1(Ω), h ⊂ Y = y C(Ω¯): y P1(T) T . h h h h { ∈ ∈ ∀ ∈ T } As usual, we will abbreviate Y = Y H1(Ω). For the control we use the space U h0 h ∩ 0 h of continuous piecewise linear functions on Γ U = u C(Γ): u P1(T Γ) T . h h h h { ∈ ∈ ∩ ∀ ∈ T } Optimization for Dirichlet Control Problems 7 Notice that the elements of U are the traces of the elements of Y . We will denote h h I : C(Ω¯) Y or I : C(Γ) U the interpolation operator related to these spaces h h h h → → and Π : L2(Ω) Y or Π : L2(Γ) U the projection onto this spaces in the L2 h h h h → → sense. For all y L2(Ω) and all u L2(Γ): ∈ ∈ (Π y,y ) = (y,y ) y Y and (Π u,u ) = (u,u ) u U . h h Ω h Ω h h h h Γ h Γ h h ∀ ∈ ∀ ∈ We discretize the state equation following the work by Berggren [27]: define S : h L2(Γ) L2(Ω) such that for u L2(Γ), y = S u if and only if y Y is the unique h h h h → ∈ ∈ solution of a(y ,z )= 0 z Y , (y ,v ) = (u,v ) v U , (10) h h h h0 h h Γ h Γ h h ∀ ∈ ∀ ∈ where a(, ) is the bilinear form associated to the operator in the PDE. In the case of · · the Laplace operator a(y,z) = Ty zdx. The discrete functional is thus defined Ω∇ ∇ as R 1 ν J (u) = (S u y ,S u y ) + (u,u) . h h Ω h Ω Ω Γ 2 − − 2 We define now Uh = u U : α(x ) u (x ) β(x ) j B , α,β { h ∈ h j ≤ h j ≤ j ∀ ∈ } and Kh = y Y : a(x ) y (x ) b(x ) j J . a,b { h ∈ h j ≤ h j ≤ j ∀ ∈ } The discrete problems reads thus as (PU) min J (u ), (PC) min J (u ), (PS) min J (u ), (PCS) min J (u ). h h h h h h h h h h h h uh∈Uh u∈Uαh,β Shuh∈Kah,b uh∈Uαh,β S u ∈Kh h h a,b The adjoint operator of S is given by the discrete variational normal derivative. h See [1]. S∗ : L2(Ω) L2(Γ) and w = S∗y if w U satisfies h → h h h ∈ h (w ,z ) = (y,z ) a(z ,φ ) for all z Y , (11) h h Γ h Ω h h h h − ∈ where φ Y is the unique solution of h h0 ∈ a(z ,φ ) = (y,z ) for all z Y . (12) h h h Ω h h0 ∈ It is customary to write S∗y = ∂hφ . For u L2(Γ) and y L2(Ω), we have h − n h ∈ ∈ (S u,y) = (u,S∗y) . h Ω h Γ We can then write 1 J (u) = (S∗S u+νu,u) (S∗y ,u) +c . h 2 h h Γ− h Ω Γ Ω Optimization for Dirichlet Control Problems 8 ThecomputationofthederivativeofJ atapointu L2(Γ)inthedirectionv L2(Γ) h ∈ ∈ is then obtained with J′(u)v = (S∗S u+νu,v) (S∗y ,v) . h h h Γ − h Ω Γ Since our final objective is to optimize in U , let us see in more detail how to make h the computations in this space. Consider e N the nodal basis in Y , where N is the dimension of Y and satisfies { j}1 h h N = NI+NB,thelatter beingrespectively thenumberofinteriorandboundarynodes. With an abuse of notation, we will also denote e the restriction of e to Γ. Define the j j usual finite element stress, mass and boundary mass matrices by = a(e ,e ), = (e ,e ) , = (e ,e ) for 1 i,j N. i,j i j i,j i j Ω i,j i j Γ K M B ≤ ≤ We will also use RN×N for the identity matrix, RN×N for the zero matrix I ∈ O ∈ and usually refer to submatrices as I,I or I,: defined by taking the rows or columns K K designed by the sets of indexes in the subscripts, the semicolon meaning “all the indexes”. For instace, the usual boundary mass matrix is B,B. B Given u = u e , we denote u RNB×1 the vector (u ,...,u )T and for h j∈B j j ∈ 1 NB y = N y e Pwe denote y RN×1 the vector (y ,...,y )T. Using (10), we have h j=1 j j ∈ 1 N that yPh = Shuh iff KI,I KI,B yI = 0 . (13) (cid:20) B,I B,B (cid:21)(cid:20) yB (cid:21) (cid:20) B,Bu (cid:21) O B B Since B,B is nonsingular, we can write this as B KI,IyI = −KI,Bu, (14) yB = u. If we define RNB×N as S ∈ −1 = KI,I KI,B :,B S (cid:20) B,I B,B (cid:21) I O I we have that y = S u if and only if y = u. h h h S Given y Y , let φ = φ e be the solution of (12) for y = y . Denoting h ∈ h h j∈I j j h φ RNI×1 the correspondingPvector, it can be computed as the solution of ∈ I,Iφ = I,:y. (15) K M Then we could compute w RNB×1, the vector whose components are the coefficients ∈ of the (minus) discrete normal derivative ∂hφ = S∗y = w = w e solving − n h h h h j∈B j j the system P B,Bw = B,:y B,Iφ. (16) B M −K Formally, we can also write that w = −1 T y. To finish this section we also define BB,BS M Optimization for Dirichlet Control Problems 9 the matrix RNB,×NB and the vector f RNB×1 by A∈ ∈ = (S∗S e +νe ,e ) and f = (S∗y ,e ) . Ai,j h h i i j Γ i h Ω i Γ We have that 1 J (u )= uT u fTu+c . (17) h h Ω 2 A − We also note that = T +ν BB. (18) A S MS B To compute the vector f, we consider the projection of y onto Y in the L2(Ω) sense, Ω h y = Π y and denote y RN×1, the vector whose j-th component is y (x ), Ω,h h Ω Ω Ω,h j ∈ and f = T y . (19) Ω S M So for u ,v U , the latter represented by the vector v, we can compute h h h ∈ J′(uh)vh = vT B,Bw+νvT B,Bu vTf. B B − Notice that applying the equality (16), the explicit computation of w is not necessary, and we can write J′(uh)vh = vT( B,:y B,Iφ+ν B,Bu f). M −K B − 3. Unconstrained problem Problem (PU) has a unique solution u¯ L2(Γ) that satisfies J′(u¯) = 0. For every ∈ 0 < h < h , problem (PU) has also a unique solution u¯ that satisfies J′(u¯ ) = 0. The 0 h h h h problems being convex, these conditions are also sufficient. Moreover it is known that u¯ u¯ strongly in L2(Γ) and also error estimates are available in terms of the mesh h → sizehifν > 0andthedomainiseither2Dandpolygonalorsmooth.See[1,2,4–6,28]. LetuscommenttwodifferentapproachesforthesolutionoftheFEMapproximation of (PU). 3.1. Solve the optimality system for the state, the adjoint state and the control. We write first order optimality conditions for (PU). There exists some h > 0 such h 0 that for every 0 < h < h , there exist uniqueu¯ U , y¯ Y and ϕ¯ Y satisfying 0 h h h h h h0 ∈ ∈ ∈ the optimality system: a(y¯ ,z ) = 0 for all z Y , (20a) h h h h0 ∈ y¯ u¯ on Γ, (20b) h h ≡ a(z ,ϕ¯ ) = (y¯ y ,z ) for all z Y , (20c) h h h Ω h Ω h h0 − ∈ νu¯ ∂hϕ¯ on Γ. (20d) h ≡ n h Optimization for Dirichlet Control Problems 10 Taking into account the definition of discrete variational normal derivative and rela- tions (14)–(16), we can write this optimality system as I,I I,B I,B I,I yI 0I K O K O  OB,I IB,B −IB,B OB,I  yB  =  0B . u y I,I I,B I,B I,I I,: Ω  −M −M O K    −M   B,I B,B ν B,B B,I  ϕI   B,:yΩ   M M B −K    M  We may eliminate u with the boundary condition, and reorder the equations to get the linear system M+νB −K:,I y = MyΩ . (21) (cid:20) I,: I,I (cid:21)(cid:20) ϕI (cid:21) (cid:20) 0 (cid:21) −K O Notice that system (21) is solvable even for ν = 0. Solving this equation would com- pletely solve the problem for the unconstrained case. When the discretization is very fine, the number of unknowns may make the solution of the system by direct methods too difficult. The preconditioned conjugate gradient method for this kind of linear systems is studied by Scho¨berl and Zulehner in [29] and Herzog and Sachs in [30]. A preconditionercan bebuiltusingmatrices representingscalars productsinY andY . h h0 Following Algorithm 1 in the last-mentioned reference, at each iteration three linear systems must be solved: two of size N and one of size NI. In [30, Algorithm 3], the systems are solved using a multigrid method. Reduced problems in the adjoint state variable and related pcg have also been studied in [31]. Nevertheless, the structure of the matrix we have in (21) is different to theonetreated inthatreferenceandapplication oftheir resultsisnotstraightforward. 3.2. Use an iterative method to solve a reduced problem in the control variable. Let us see another way of solving (PU). Using (17), first order optimality conditions h can also be written as u= f. (22) A Lemma 3.1. The matrix is symmetric and positive definite for all ν 0 and there A ≥ exists C > 0 independent of h and ν such that its condition number is bounded by κ( ) CλNλB1((MBB,)B)λN(M)+νλNB(BB,B), (23) A ≤ λ1( )+νλ1( B,B) M B where 0 < λ1(BB,B) < λNB(BB,B) and 0 < λ1(M) < λN(M) are respectively the smallest and greatest eigenvalues of the matrices B,B and . B M Proof. It is clear from (18) that is symmetric. The mass matrices and B,B are A M B symmetric and positive definite and therefore λ1( ) > 0 and λ1( B,B) > 0. Since the M B boundarycomponents of Su are exactly u, we have that u u and hence kS kRN ≥ k kRNB uTAu = uT(STMS+νBB,B)u = (Su)TM(Su)+νuTBB,Bu) ≥ (λ1(M)+νλ1(BB,B))kuk2RNB (24)