SECOND ORDER AND STABILITY ANALYSIS FOR OPTIMAL SPARSE CONTROL OF THE FITZHUGH-NAGUMO EQUATION

Abstract. Optimal sparse control problems are considered for the FitzHugh-Nagumo system including the so-called Schlögl model. The non-differentiable objective functional of tracking type includes a quadratic Tikhonov regularization term and the L1-norm of the control that accounts for the sparsity. Though the objective functional is not differentiable, a theory of second order sufficient optimality conditions is established for Tikhonov regularization parameter ν > 0 and also for the case ν = 0. In this context, also local minima are discussed that are strong in the sense of the calculus of variations. The second order conditions are used as main assumption for proving the stability of locally optimal solutions with respect to ν → 0 and with respect to perturbations of the desired state functions. The theory is confirmed by numerical examples that are resolved with high precision to confirm that the optimal solution obeys the system of necessary optimality conditions.

Introduction.This paper contributes to several recent developments in the optimal control of partial differential equations.First, and this is perhaps the main novelty of our work, it improves the theory of second order sufficient optimality conditions for the optimal control of nonlinear evolution equations.We prove results that are not only new for the control of the FitzHugh-Nagumo equations.They are also not yet known for the optimal control theory of standard academic semilinear elliptic or parabolic control problems with smooth objective functionals.We even go beyond this and consider problems with nonsmooth functionals in the context of sparse optimal control, where the L 1 -norm of the control appears in the objective functional.
Our control system of FitzHugh-Nagumo equations plays an important role in physics, chemistry, and mathematical biology.Here, we continue our research in [8] on the optimal control of wave-type solutions such as traveling waves or spiral waves that appear as typical solutions of this system in unbounded spatial domains.Investigating this class of control problems with its interesting applications, we observed numerical effects that we wanted to confirm by a deeper mathematical analysis.Eventually, we arrived at novel second order conditions that we needed to explain our numerical observations that are examplarily presented at the end of this paper.
Let us detail our main achievements in the theory of second order sufficient optimality conditions a bit more: Second order conditions are related to certain critical cones that must be chosen as small as possible.For sparse controls, in Lemma 3.1 (3) we introduce the cone of critical directions C ūν that is smaller than the associated one suggested in [7].In this way, we improve the results of [7] that were developed for a simple class of elliptic equations.Even for this class of simpler problems, our result is new.
Moreover, our paper contains second order conditions that are sufficient for strong local minima in the sense of calculus of variations.So far, almost all contributions to the theory of second order conditions in the optimal control of PDEs addressed weak local minima.To our best knowledge, the first result on strong local minima in PDE control was recently obtained in [1] for the case of semilinear elliptic equations.Although our results are similar to the ones of [1], they are more general, even if they are transferred to the elliptic case discussed there.In particular, we introduce two extended cones of critical directions, C τ ū in (3.14) and E τ ū prior to Theorem 3.13.They can be used to deal with the case of vanishing Tikhonov regularization parameter ν; see the more detailed remarks in the introduction to section 3. We do not know papers on the control of PDEs that deal with strong local minima for the case ν = 0.
The case ν = 0 is important for the numerical application of control methods to the FitzHugh-Nagumo equations that we present in the last section.In our preceding paper [8], we observed that the numerical methods worked fairly stable also for very small parameter ν > 0 so that we became interested in the convergence analysis as ν → 0. To prove the observed stability of optimal controls for ν → 0, we had to develop our new second order conditions that are applicable to the case ν = 0.This is a situation where strong local minima are needed.We also discuss the stability of optimal solutions with respect to perturbations of the given desired state functions.
Compared with standard semilinear elliptic or parabolic equations, the analysis of the FitzHugh-Nagumo system is more difficult.Therefore, we have to prove the fairly technical Lemmata 3.8-3.11.Although these results do not directly extend the optimal control theory, their detailed proofs are unavoidable for building our theory and occupy a major part of our paper.
The analysis of optimal control problems for FitzHugh-Nagumo systems was already considered in [3], which investigated associated problems by the Dubovitskij-Milyutin optimality conditions, and in [13], which concentrated on the analysis of time-optimal control problems for a linear version of the FitzHugh-Nagumo equations.Later, resolving certain obstacles in the analysis of differentiability, in [8] we proved the first order optimality conditions for the (nonlinear) FitzHugh-Nagumo equations by showing the differentiability of the control-to-state mapping of any order.Moreover, we reported on a variety of computational results.In this context, we also mention [4], where the analysis and numerical treatment of optimal control problems for traveling wave fronts is discussed for the so-called Schlögl model (also known as the Nagumo equation).
The discussion of control and feedback control problems of reaction-diffusion equations has a long tradition in the community of theoretical physics.We mention, for instance, contributions in [15], [17], [20], and the references therein.We also refer to the survey volume by [16].Our investigations on the control of FitzHugh-Nagumo equations and systems of related type were initiated by these ideas.
In this system, the partial differential equation for y is said to be the activator equation, while the one for z is called the inhibitor equation.The function y is the state that is to be controlled.In our paper, the inhibitor z has only some auxiliary character with respect to the control.For the choice α = 0, both equations decouple and the state function y has to solve the Schlögl equation.We should mention that this equation is also known as the Nagumo equation and it is also a particular case of the Allan-Cahn equation.Here, the inhibitor equation is meaningless.For α = 1, the FitzHugh-Nagumo system is obtained.These values α ∈ {0, 1} are the values of our interest, but the analysis for the system (1.1) is the same for any arbitrary real α.
Following the usual notation we set The following theorem was proved in [8, Theorem 2.1 and Corollary 2.1].
Following again [8], we introduce the mapping The second derivative (ω v1,v2 , χ v1,v2 is given by the pair (ω, χ) solving the equation where and G (u) can be extended to continuous linear and bilinear mappings from L 2 (Q) to W (0, T ) 2 for any u ∈ L p (Q) with p > 5  2 .Now, we formulate our control problem In this setting, we have given constants , and the set of admissible controls and C Z T are nonnegative.Remark 1.4.We introduced fairly general coefficients in the objective functional F ν .The reason for that is the fairly different applications we have in mind.For Downloaded 11/06/15 to 193.144.185.28.Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.phpinstance, in [8] we considered the problem of extinguishing a moving spiral wave.To do so, we observed only the tip of the desired spiral wave that moved from the center to the boundary of the domain Ω.Here, C Y Q was the characteristic function of a circle around the tip that moved with the tip through Ω.We also mention that our theory works under the slightly weaker assumption y T , z T ∈ L p (Ω), too.However, this would lead to more technicalities, hence we do not use this assumption.
Thanks to Theorems 1.1 and 1.2, it is easy to prove that ∀ ν ≥ 0 the control problem (P ν ) has at least one solution ūν ; see [8,Theorem 3.1].
We finish this section by analyzing the cost functional J ν that consists of two terms having different smoothness.While the first part is smooth, the second part, j : It admits directional derivatives given by the formula (1.4) where Of course, the relations defining these sets are required only almost everywhere.Moreover, we know that a measurable function λ belongs to the subdifferential in the sense of the convex analysis of j at a point u, λ ∈ ∂j(u), if and only if it satisfies for a.a.
For every 0 < ρ < 1 we have These inequalities are valid for any convex and locally Lipschitz functional j; see, for instance, [2].The first part of J ν has the following differentiability properties.Theorem 1.5.The function and we have for the first and second derivatives the expressions where ϕ u along with ψ u are the solutions in W (0, T ) ∩ L ∞ (Q) of the adjoint system The proof of the existence and uniqueness of the solution (ϕ u , ψ u ) ∈ [W (0, T ) ∩ L ∞ (Q)] 2 of the adjoint system and the formula (1.6) can be found in [8, section 3.2].The expression of the second derivative follows from the chain rule, (1.3), and (1.8).
2. First order optimality conditions.Since the control problem (P ν ) is not convex, we must consider local minima.In this section, the goal is to set up the first order necessary optimality conditions satisfied by the local minima and to draw some conclusions from the optimality system.We say that ūν is a local minimum of problem (P ν ) in the L p (Q) sense, 1 ≤ p ≤ +∞, if there exists ε > 0 such that where B ε (ū ν ) denotes the L p (Q)-ball centered at ūν with radius ε.Let us mention that the boundedness of U ad in L ∞ (Q) implies that ūν is a local minimum in the L 2 (Q) sense if and only if it is a local minimum in the L p (Q) sense for any 1 ≤ p < +∞.On the other hand, if ūν is a local minimum in the L ∞ (Q) sense, then it is a local minimum in the L p (Q) sense for any 1 ≤ p ≤ +∞.Hereafter, local minima will always be understood as local minima in the L 2 (Q) sense.
Remark 2.1.Minima of this type are, viewed in the sense of calculus of variations, weak local minima.They ensure optimality in a neighborhood of the locally optimal control ūν .Later, we will also investigate conditions for local minima that are strong in the sense of calculus of variations.In that case, we have no matter how far u is from ūν .Now we state the optimality conditions satisfied by a local minimum of (P ν ).Theorem 2.2.Let ūν be a local minimum of (P ν ).Then, there exist (ȳ ν , zν ) and Then, for any u ∈ U ad , u = ūν , and ∀ 0 < ρ < ε/ u − ūν L 2 (Q) , we have that ūν + ρ(u − ūν ) ∈ U ad ∩ B ε (ū ν ).Hence, using the convexity of j, we get Now, passing to the limit ρ → 0, we infer Taking (ȳ ν , zν ) and ( φν , ψν ) as the solutions of (2.1) and (2.2), respectively, and using the expression (1.6) of the derivative of F ν , we deduce from the above inequality that Therefore, ūν solves the convex optimization problem that can be considered in L 1 (Ω).Finally, from rules of subdifferential calculus we deduce the existence of λν ∈ ∂j(ū ν ) such that (2.3) holds.
From inequality (2.3), we deduce some interesting relations.Corollary 2.3.Let ūν , φν , and λν be as in Theorem 2.2, and assume that ν > 0. Then the following relations hold: Moreover, from the last two representation formulas it follows that λν is unique and Corollary 2.4.Let ūν , φν and λν be as in Theorem 2.2 and assume that ν = 0. Then the following relations are satisfied: The last representation formula yields that λν ∈ L 2 (0, T ; For the proofs of these relations, the reader is referred to [7, Corollary 3.2] and [6, Theorem 3.1], which deal with the cases ν > 0 and ν = 0, respectively, in a similar situation.The relations in the corollary show that local minima ūν can be zero in large regions, the support of ūν being monitored by κ. In particular, a value κ 0 < ∞ exists such that for κ > κ 0 we have only one local minimum, namely, ūν = 0. Indeed, since U ad is bounded in L ∞ (Q), the set {(y u , z u )} u∈U ad is bounded in L ∞ (Q), and hence from (1.8) we deduce that the set {(ϕ u , ψ u )} u∈U ad is also bounded in L ∞ (Q).Then, our statement holds obviously for In the case ν = 0, the control ūν can admit values outside of {a, 0, b} on a set of positive measure only if the set {(x, t) ∈ Q : | φν (x, t)| = κ} has a positive Lebesgue measure.If the measure of this set is zero, then ūν (x, t) must belong to {a, 0, b} for almost every (x, t) ∈ Q.Such optimal controls ūν (x, t) are of "bang-bang-bang" type.

Second order optimality conditions.
In this section, we carry out the second order analysis for (P ν ).First, we establish second order necessary conditions.Distinguishing between the cases ν > 0 and ν = 0, we consider sufficient conditions.The reader might be surprised how different both cases can be.
Let us give at this point a short orientation on the second order conditions.To deal with the different situations depending on ν, we will introduce three different cones of critical directions, namely, the cone C ū and the cones C τ ū and E τ ū .The cone C ū is quite standard and well known.It will be used in the case ν > 0. In Theorem 3.4 we prove that sufficient second order conditions based upon C ū ensure (weak) local optimality.However, we show just as a corollary that then ūν yields even a strong local minimum.
The two other cones are needed for the case ν = 0.The second order conditions based on E τ ū imply that ū is a strong local minimum in the sense of calculus of variations; cf.Theorem 3.13 and relation (3.38).This property of strong local optimality cannot be deduced if we replace E τ ū by C τ ū .Section 4 is devoted to the stability of locally optimal solutions with respect to perburbations.For ν > 0 the cone C τ ū can be used to show Lipschitz stability of optimal controls; cf.Theorem 4.2.Here, weak local minima are useful.
In section 4.1.2,we discuss the stability of locally optimal solutions for ν = 0. Here, we invoke the second order conditions based on E τ ū and prove the stability of the optimal states with respect to perturbations of the desired state.Now, stability can be proved only if ū is a strong minimum.For weak minima we cannot prove this property.
We should mention that this strong minimum property would also be needed if we want to prove error estimates for the optimal states when the control problem is discretized.Indeed, we only have weak convergence of the optimal discrete controls to the continuous ones, but the corresponding states have strong convergence properties, which allows the use of the sufficient second order condition.
The reader should observe that the contribution of the term j(u) in the second order optimality conditions is through the critical cone C ūν , but only F ν (ū ν )v 2 is involved in the second order approximation.In some sense, j is piecewise linear and there is no second order contribution from it.Next we consider the sufficient second order conditions.
Let us now start the second order analysis.The presence of the so-called Tikhonov term ν 2 u 2 L 2 (Q) in the cost functional is extremely important.Due to this term, the second order sufficient conditions for the local optimality of ūν are a straightforward extension of the corresponding conditions for a finite dimensional optimization problem.Actually, it is well known that the equivalence stated in the following theorem is not in general true for infinite dimensional optimization problems.However, it works perfectly in this case because of the presence of the Tikhonov term.
Theorem 3.3.The following statements are equivalent: (2) There exists σ > 0 such that We omit the proof of this result, because it was proved for a certain class of parabolic partial differential equations in [9,Theorem 4.11].For the FitzHugh-Nagumo system, the situation is slightly more involved, but the extension of the proof in [9] to this system is more or less straightforward.
Proof.This proof follows the classical one that is performed by contradiction.The only difference to take into account is the nondifferentiable term j involved in the cost functional and in the definition of the cone of critical directions.Let us sketch some details of the proof.Arguing by contradiction, we assume that there exists a sequence Now, we take the last convergence after possibly selecting a subsequence of {v k } ∞ k=1 .The proof is split into three parts.
(i) v ∈ C ūν .Since every v k satisfies the sign conditions of the definition of C ūν , we deduce that v also satisfies them.Let us prove that F ν (ū ν )v + κj ν (ū ν ; v) = 0. From (1.5) and (2.3) we get that To prove the converse inequality, we apply the mean value theorem and use (3.7), k 2k with some 0 ≤ θ k ≤ 1.Therefore, with (1.5) we obtain from above Now by convexity and continuity of the mapping Here, we have used that Performing a Taylor expansion and using (3.7) again, we find From (1.5), we obtain . Thanks to the expression of F ν given in (1.7), it is easy to pass to the limit with the aid of the compactness of the linear operator k = 0. Now we prove the surprising fact that the local optimality ensured by Theorem 3.4 is even strong.Corollary 3.5.Let ūν satisfy the assumptions of Theorem 3.4.Then, there exist δ > 0 and ε > 0 such that Proof.Let us assume that (3.9) does not hold for any δ and ε .Then, for any integer k ≥ 1, we can find a control We can take a subsequence, denoted in the same way, such that Indeed, from the boundedness of {u k } k≥1 in L ∞ (Q) and Theorem 1.1 we infer the boundedness of {(y k , z k )} k≥1 in W (0, T ).Therefore, taking again a subsequence if necessary, we can assume the weak convergence of {(y k , z k )} k≥1 in W (0, T ) and of {u k } k≥1 in L p (Q) for every p < +∞.The weak limit of {y k } k≥1 is obviously ȳν .From the state equation satisfied by the functions z k we deduce that z k zν in W (0, T ).Finally, passing to the limit in the first equation of the system (1.1) we conclude that u k ūν in L 2 (Q).From y u k → ȳν in L ∞ (Q) and the last two equations of (1.1) we also obtain that z u k → zν in L ∞ (Q).Now, we deduce from (3.10) Passing to the limit, we get Notice that F 0 (u k ) → F 0 (ū ν ) follows from the strong convergence of y u k and the circumstance that F 0 does not explicitly depend on u.

3.2.
The case ν = 0.In this case, we write (P) instead of (P 0 ), J and F instead of J 0 and F 0 , and ū, ȳ, φ, etc., instead of ū0 , ȳ0 , φ0 , etc.In general, for an infinitedimensional optimization problem, the strict positivity of the second derivative of the functional on the critical cone is not enough for local optimality.The reader can find an example for this fact in [12].Therefore, we have to consider an extended cone.Given 0 < τ < κ, we define the cone C τ ū as the set of elements v ∈ L 2 (Q) satisfying Proposition 3.6.The extended cone C τ ū has the following properties: is fulfilled, where Q v denotes the set of points (x, t) ∈ Q such that (3.12) is not satisfied by v(x, t).
From (3.4) and (3.12) we infer that C τ ū is a small relaxation of C ū when ν = 0.In the finite dimensional case, both cones are the same for all sufficiently small τ .Now, one might be tempted to formulate the second order sufficient conditions as (3.14) ∃τ > 0 and ∃σ > 0 such that Downloaded 11/06/15 to 193.144.185.28.Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.phpHowever, for ν = 0, this condition can hold only in exceptional cases, as we will prove at the end of this section.In fact, it is enough to take a look at the expression (1.7) of F (ū)v 2 to understand why the condition (3.14) cannot be expected for ν = 0.The next theorem provides the correct sufficient condition.
Before proving this theorem, we derive some auxiliary results.Lemma 3.8.There exist constants C a,b , C 1 , C 2 , C 3 , and C ∞ such that ∀ u ∈ U ad we have the estimates Proof.The first estimate follows from Theorem 2.1 and Corollary 2.1 of [8].To prove the second and third estimates, we introduce the family of operators for μ ≥ 0. Now, we take w(x, t) = e −μt (y u − ȳ)(x, t).Subtracting the equations satisfied by y u and ȳ and observing that (z u − z)(x, t) = K 0 (y u − ȳ)(x, t), we obtain the following equation for w: The estimation of ϕ u − φ is performed for the version of the adjoint equation where ψ u and ψ are eliminated on using the adjoint integral operator K * μ .We omit the related computations, because they are straightforward.
Finally, we prove the estimate (3.20).From [8, (10)], we know the estimate Multiplying the equation by w and integrating in Q we infer From here we get Writing η v (x, t) = e μt w(x, t) we deduce (3.20).Lemma 3.9.There exist a constant C η and ε > 0 such that where from now on, differing from the notation in Theorem 1.2, η u,v and η v are the first components of G (u)v and G (ū)v, respectively.Moreover, there is a constant C 0 such that Proof.We introduce the operator K μ as in the proof of Lemma 3.8 and define w(x, t) = e −μt (η u,v − η v ).Now, we subtract the equation satisfied by η v from that for η u,v and apply the mean value theorem to get where ỹu (x, t) = ȳ(x, t)+ϑ(x, t)(y u (x, t)− ȳ(x, t)), 0 ≤ ϑ(x, t) ≤ 1, and μ is sufficiently large.Invoking again [8, Lemma 2.3], we obtain (3.22).
To show (3.23), we proceed similarly.We set w = e −μt (y u − ȳ −η u−ū ) and perform the Taylor expansion From the PDEs for y u , ȳ, and η u−ū , we deduce analogously to (3.24)  Proof.This estimate follows easily from the expression (1.7) of F , (3.17), and the fact that the functions C Y Q , C Y T , C Z Q , and C Z T are bounded.In addition, we have to recall that ζ v (x, t) = K 0 (η v )(x, t), where K 0 is the operator introduced in the proof of Lemma 3.8.From this identity we get Next, we prove our last lemma.Lemma 3.11.For every ρ > 0 there exists ε > 0 such that Let us indicate how all these terms can be estimated.First, from (3.22) we deduce that From here, we also get We use (3.28) and (3.30) to estimate I 1 : I 3 is estimated from (3.29) and (3.31) as follows: Now, we proceed in the same way as for the estimates I 1 and I 3 .Finally, we estimate I 2 as follows:

By (3.19), we can estimate the first two terms by C y
for some constant C .The third term is estimated as I 1 .The statement of the lemma is a straightforward consequence of the obtained estimates.
Proof of Theorem 3.7.By Lemma 3.11, we have the existence of Using Lemma 3.8, we deduce that Let us take From here it follows Inserting this inequality in (3.34) and using (3.15) and (3.25) it follows that and, with Young's inequality, Inserting this inequality in (3.35) we obtain Proof.We define w = e −μt (y u − ȳ − η u−ū ); then we know the estimate (3.23), Let now ε be as in Theorem 3.7 and take hence, moving the term 1 2 y u − ȳ L 2 (Q) to the other side, Moreover, we have Therefore, we get from the last two estimates ∀u ∈ U ad such that y u − ȳ L ∞ (Q) < ε.Finally, we take 0 < δ ≤ σ 8C3 and 0 < ε ≤ ε 0 so that y u − ȳ L ∞ (Q) < ε 0 for every u ∈ U ad ∩ B ε (ū).Then, (3.36) follows from (3.16).
The growth condition (3.36) is valid in a ball around ū. Therefore, we analyze if a result similar to Corollary 3.5 can be proved for the case ν = 0.This would be applicable in a ball around ȳ, hence in a possibly larger neighborhood of ū.
We have not been able to establish such a result based on the cone C τ ū .To deal with this problem we introduce a different extended cone defined by From Lemma 3.1(1) and (1.5) we infer that C ū ⊂ E τ ū for every τ > 0. Thus the cone E τ ū can be considered as a small extension of C ū.We are able to prove the following result on second order sufficiency that is based on E τ ū .Theorem 3.13.Let ū ∈ U ad satisfy the optimality system (2.1)-(2.3)along with the state (ȳ, z) and the adjoint state ( φ, ψ).Assume also that ∃τ > 0 and ∃σ > 0 such that (3.37) Downloaded 11/06/15 to 193.144.185.28.Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.phpThen, there exists ε > 0 such that Proof.First, we prove that there is a constant M > 0 such that With Lemma 3.9, (3.23), and (3.17) we obtain Hence, (3.39) holds with where M a,b was introduced in Lemma 3.10.From Lemma 3.11, we deduce the existence of ε 2 > 0 such that ∀u ∈ U ad with Next, we prove that (3.38) holds with ε = min{ε 1 , ε 2 }.We take u ∈ U ad such that Here we have Notice that (3.39) yields ( Corollary 3.14.Let ū satisfy the assumptions of Theorem 3.13.Then, there exist δ > 0 and ε > 0 The proof of this corollary is almost the as that of Corollary 3.12.Corollary 3.12 also holds under the second order condition (3.37).This is a consequence of (3.41) and (3.19).By Corollary 3.14 the local optimality of ū is strong in the sense of calculus of variations.
We finish this section by proving that the second order condition (3.14) can only hold in some exceptional cases.
Proof.We argue again by contradiction.Assume that (3.14) holds.Then, arguing as in the proof of Theorem 3.7, we deduce the existence of δ > 0 and ε > 0 such that Then, ū is a solution of the problem H(x, t, y, z, ϕ, u) = L(x, t, y, z

The Hamiltonian of this control problem is given by
where and H denotes the part of the Hamiltonian involving u.  [5] or [14,Chapter 4] for the proof of Pontryagin's principle for control problems associated with evolution partial differential equations.However, (3.42) contradicts the following inequalities: Indeed, let us prove that the first and third statements contradict (3.42).The other two statements are analyzed similarly.If κ < φ(x, t), then Corollary 2.4 implies that ū(x, t) = a.Now, using that φ(x, t) < κ Let us consider the case 0 < φ(x, t) < κ.Then Corollary 2.4 implies that ū(x, t) = 0. On the other hand, using the inequality Finally, since φ ∈ C(Q), and ū ≡ 0, ū ≡ a, ū ≡ b, φ ≡ κ, and φ ≡ −κ, we deduce from Corollary 2.4 that the sets are strict open subsets of Q, and at least one of them is not empty.Hence, the set of points (x, t) satisfying one of the conditions of (3.43) is open and nonempty.Consequently, the set of points where Pontryagin's principle fails has a positive Lebesgue measure, which is not possible.Remark 3.16.The coercivity assumption (3.15) can hold only if there is a constant δ 0 > 0 such that C Y Q (x, t) ≥ δ 0 holds for a.a.(x, t) ∈ Q and C Y T (x) ≥ δ 0 holds for a.a.x ∈ Ω.Our theory is worked out for this case.

If C Y
T does not obey this assumption, for instance, if C Y T = 0, then the main results must be modified in an obvious way.The reader will confirm that then the theory remains valid under the following changes: The norm squares of η v (T ), η u−ū (T ) in (3.15), (3.16), (3.37), and (3.38) are to be deleted.Analogously, the norm squares of y u (T ) − ȳ(T ) must be deleted in (3.36) and (3.41).

Applications to the stability analysis with respect to perturbations.
Here, we exemplarily show by three case studies how the second order sufficient optimality conditions can be invoked as the main assumption for proving the stability of optimal solutions with respect to perturbations of our optimal control problem.We consider perturbations of the desired state functions and also the case where the Tikhonov regularization parameter ν tends to zero.This can be viewed as a perturbation of the reference parameter ν = 0.
is satisfied.We also assume that {y ε T , z ε T } ε>0 is bounded in L ∞ (Ω).Associated with these perturbed data, we define the perturbed objective functionals, , and the family of perturbed optimal control problems, Prior to proving the next result, for convenience we introduce the notation . These functions contain only the parts of J and J ε that depend on the state functions.The Tikhonov regularization term and κ j are separated from them.
Theorem 4.1.If {ū ε ν } ε is any sequence of optimal controls of problems (P ε ν ) that converges weakly in L 2 (Q) to some ūν , then ūν is optimal for (P ν ) and Conversely, if ūν is a strict locally optimal control of (P ν ), then there exists a sequence {ū ε ν } ε of locally optimal controls of (P ε ν ) converging to ūν .This sequence obeys (4.2).Furthermore, there exists ρ > 0 such that every ūε ν affords a global minimum to J ε ν in U ad ∩ Bρ (ū ν ), where B ρ (ū ν ) denotes the L 2 (Q)-ball.Proof.Let us skip the subscript ν in the proof for convenience.We write ū := ūν and ūε := ūε ν .From ūε ū we get ū ∈ U ad and, in a standard way, that ū is optimal for (P ν ).It remains to show the strong convergence.
Applying a result of [11], we know that ȳε − ȳ is bounded in some Hölder space C λ ( Q), λ ∈ (0, 1).Therefore we obtain by compact embedding in C( Q) that a subsequence ȳε := y ūε converges strongly in L ∞ (Q) toward ȳ as ε → 0. Since the same holds for all subsequences, we even have ȳε − ȳ L ∞ (Q) → 0 as ε → 0 for the whole sequence.From (1.1), we also deduced that zε for short, we find . Downloaded 11/06/15 to 193.144.185.28.Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php In the same way, all other integrals appearing in (4.4) are first simplified and next estimated by the Cauchy-Schwarz inequality.Now the result follows immediately from these estimates and (4.4).

4.1.2.
Case ν = 0. Now we consider the problem (P ν ) for ν = 0 and we write for short (P):= (P 0 ).Let us prove the following result analogous to Theorem 4.1.
Theorem 4.3.If {ū ε } ε is any sequence of optimal controls of problems (P ε ) that converges weakly in L 2 (Q) to some ū, then ū is optimal for (P) and Conversely, if ū is a strict locally optimal control of (P), then there exists a sequence {ū ε } ε of locally optimal controls of (P ε ) converging weakly to ū.Moreover, (4.5) holds.Furthermore, there exists ρ > 0 such that every ūε affords a global minimum to J ε with respect to the elements u ∈ U ad such that y u − ȳ ∞ ≤ ρ.
Proof.The first part of the theorem can be proved arguing as in the proof of Theorem 4.1.For the second part, we take ρ > 0 such that ū is the minimizer of the functional J in the convex and closed set Now we take ū as a global minimizer of J ε in K ρ .Arguing again as in the proof of Theorem 4.1, we deduce that {ū ε } ε converges weakly to ū in L p (Q) and (4.5) holds.
Theorem 4.4 (Lipschitz stability for ν = 0).Let ū be a locally optimal control of (P) that satisfies the second order sufficient optimality condition (3.37) and let {ū ε } be a sequence of locally optimal controls of (P ε ) that converges weakly to ū in L p (Q) as ε → 0 with the properties established in Theorem 4.3.Denote the associated states by (ȳ, z) and (ȳ ε , zε ), respectively.Then there exists C > 0 such that Proof.Since ȳε − ȳ L ∞ (Q) → 0 as ε → 0, for all sufficiently small ε > 0, ȳε belongs to a neighborhood of ȳ, where the quadratic growth condition (3.41) holds by Corollary 3.14.Thanks to this growth condition, we can argue as follows: After rearranging, the last inequality admits the form In the same way, all other integrals appearing in (4.7) are first simplified and next estimated by the Cauchy-Schwarz inequality.Finally, we find Here, we have exploited the estimate and an associated inequality for the time t = T (cf. the proof of Lemma 3.8) to get the second inequality and the equivalence of all norms in R 4 to obtain the third one.Now the result follows immediately.

Tikhonov parameter tending to zero.
In this section, we investigate the behavior of a sequence of optimal controls {ū ν } ν>0 of (P ν ) and the corresponding states {(ȳ ν , zν )} ν>0 as ν → 0. Since U ad is bounded in L ∞ (Q), any sequence of solutions of (P ν ) contains subsequences converging weakly * in L ∞ (Q).Below, we will deduce consequences of this convergence.
Proof.First we observe that the boundedness of {ū ν } ν>0 in L ∞ (Q) and the weak convergence ūν ū for ν ↓ 0 in L 2 (Q) implies that ūν ū for ν ↓ 0 in L p (Q) for any 1 ≤ p < ∞.Moreover, ū ∈ U ad holds.Let (ȳ ν , zν ) = G(ū ν ) and (ȳ, z) = G(ū).From the equation satisfied by zν , we deduce (4.9) zν (x, t) = e −βt z 0 (x, t) We have a similar representation for z in terms of ȳ.Using (3.17) and (3.19) and the compactness of the embedding W (0, T ) ⊂ L 2 (Q), it is easy to pass to the limit in the state equation (1.1) satisfied by (ȳ ν , zν ) and to confirm that (ȳ ν , zν ) (ȳ, z) in W (0, T ).By the continuity of the embedding W (0, T ) → C([0, T ], L 2 (Ω)), the convexity of the cost functional with respect to (y, z, u), and using that ūν is a solution of (P ν ), it follows that Hence, the inequality is fulfilled, where the right-hand side converges to zero.The proof is finished by recalling that We should mention that this estimate is fairly pessimistic.All our numerical tests below show that the convergence has the order ν, i.e., numerically we see a Lipschitz rather than a Hölder estimate.Though this is not a proof, there is a theoretical explanation for this.It is known from similar discussions for problems with linear elliptic equations that a Lipschitz estimate for ν → 0 is true provided that the adjoint state function for ν = 0 satisfies a certain assumption on the behavior in its zeros; cf.[19].This assumption is often fulfilled and the result of Lipschitz stability seems to be satisfied also in our case of a nonlinear system of parabolic equations.However, an associated discussion would go beyond the scope of our paper.

Numerical results.
The goal of this section is to confirm the convergence of ȳν against ȳ for ν → 0. This needs a very precise computation of optimal solutions and therefore we concentrate on the spatial one-dimensional (1D) case.In one dimension, the Schlögl model develops traveling wave fronts, while the FitzHugh-Nagumo system exhibits traveling pulses as typical solutions.For 2D examples, we refer to the various examples in our former paper [8].For instance, we showed the control of moving spiral waves.
We apply a semismooth Newton method since it allows us to determine solutions that satisfy the first order necessary optimality conditions with high accuracy.This method requires the solution of forward-backward parabolic systems that can be efficiently solved in one dimension but would be very demanding in two dimensions.
Let us briefly sketch our numerical approach.The optimal solutions must obey the projection formulas for the optimal control and the associated subgradient.We define active and inactive sets for the control u depending on the adjoint state ϕ u : Almost everywhere on A a , A 0 , and A b , the optimal control must be equal to a, 0, and b, respectively.Almost everywhere in I − , respectively, I + , the equation u = must hold true.By these sets, the system of necessary optimality conditions can be expressed in the form , where χ k , k ∈ {A a , I − , A 0 , I + , A b }, denotes the indicator function of the associated sets.This optimality system, along with the associated initial and boundary conditions, is solved by the semismooth Newton method.To shorten the presentation, we do not explain the steps of this standard method and refer to [18], where the method is presented in the context of sparse controls for elliptic equations.The method is terminated if a suitable norm F is sufficiently small.To avoid certain oscillation effects of the control iterates, we did not take the whole Newton step and added the computed Newton direction (d y , d z , d ϕ , d ψ , d u ) with some step-size s to the last iterate (y k , z k , ϕ k , ψ k , u k ).We invoked a modified Armijo-Goldstein rule, since the standard one slowed down the algorithm considerably.
To have a sufficiently accurate initial iterate for the semismooth Newton method, we applied a few iterations of the nonlinear CG method that we also used in [8].Moreover, the CG method was our method of choice in the case ν = 0.
To set up an optimality indicator for u, we define the function As an optimality indicator, we used where |Q| is the volume of Q.Our optimization algorithms are terminated for sufficiently small Θ. in Figure 5.1 (left) at t = 0, the target is to approach a piecewise linear shape defined by For Ω = (0, 20), partitioned by 201 equidistant nodes, the initial value of our desired wave front satisfies y 0 (10) = 1 2 .Figure 5.1 displays the natural development of this wave profile as well as the desired state y Q for the observed time horizon (0, 10).

Example 1 (Schlögl model). At first, we consider the Schlögl model and fix
As a sparse parameter, we set κ = 1e − 2 and fix the control bounds a = −1 and b = 1.For ν = 1e − 9 the optimal control is displayed in Figure 2 together with its associated state, its adjoint state, and the related sets A a , A b , A 0 that should coincide with the active/inactive sets of ū.As in the other examples, they show a very good coincidence.
Let us explain why this form of the optimal control can be expected.In the unconstrained case, we insert the desired function y Q in the left-hand side of the 1D Schlögl model, The result is a control such that the objective functional is zero for ν = 0. Notice that the first derivative of y Q is a function of Heaviside type with jumps at 9.5 and 10.5, hence the Dirac Delta functions δ 9.5 and −δ 10.5 are obtained as derivatives.The "best" unconstrained control û is a measure and not a measurable function.
With bounded and measurable control functions taking values in [a, b], the desired state y Q cannot be reached.However, the optimal controls will approximate the measure defined above.This explains why ū is of bang-bang-bang type; cf. Figure 2.
To compute the optimal solution for ν = 0, we used the nonlinear CG method and took the optimal control for ν = 1e − 9 as the initial iterate.However, the CG method did not improve the initial iterate; the optimality indicator Θ = 6.21e − 13 from (5.1) was very small.For this reason, we took ȳ := ȳ1e−9 as a reference solution to determine the order of convergence as ν ↓ 0. In both the L 2 (Q)-and L ∞ (Q)-norms, ȳν − ȳ appears to decay linearly for ν ↓ 0; cf.

Example 2 (FitzHugh-Nagumo system).
Here we consider the 1D FitzHugh-Nagumo system for Ω = (0, 75), T = 10, α = 1, β = 0, γ = 0.33, δ = −0.429,and R(y) = y(y − √ 3)(y + √ 3).We take a natural impulse as initial value; this is a snapshot of an (instationary) impulse that develops in the uncontrolled case as a typical solution of the 1D FitzHugh-Nagumo system.For κ = 0.1, a = −10, and b = 10, the goal is to double the velocity of the natural impulse in the y-component of the solution (y, z).To this aim, we fix c Y T = c Z T = c Z Q = 0 and c Y Q = 1.The natural development of the initial impulse and the desired trajectory are shown in Figure 4.
Figure 5 displays the optimal control and the adjoint state for ν = 1e − 10.Graphically, the optimal state does not differ from the desired state displayed in Figure 4.Moreover, also here the active sets of the optimal control and the sets I + , A 0 , I − coincide.We skip these figures for space reasons.The calculated control does not touch the bounds a = −10 or b = 10.Its maximal absolute value is 8.89.
As in Example 1, taking this result as initial control for the case ν = 0, the CGmethod did not further decrease the objective function; the value Θ = 9.81e − 13 of the optimality indicator already is very small.Also in this example, the decay of ȳν − ȳ in the L 2 (Q)-and L ∞ (Q)-norms seems to be almost linear.The deviation of the last value for ν = 1e − 9 is due to numerical limitations.For very small ν > 0, the accuracy cannot be expected to get any better.Notice that the term 1/ν appears in the optimality system.We dispense with an associated graphical representation as in Figure 3, as the decay looks similar.
Let us emphasize that the values of ȳν − y Q L ∞ (Q) and ȳν − y Q L 2 (Q) from Table 1 are quite large.Since the sparse optimal control acts very localized, not all features of the desired trajectory can be achieved.As Figure 6 displays, the profile of the impulse is not conserved perfectly, but the main characteristic is.Therefore, we believe the result is close to global optimality.

Example 3 (FitzHugh-Nagumo system with explicitly known optimal solution
).In the preceding example, for ν = 0 or ν very close to 0, the figures revealed areas with nonempty interior, where the optimal control took nonzero values strictly between a and b.According to the analysis, the absolute value of the adjoint state should be equal to κ in these regions.However, this is a numerical result for a finite dimensional approximation.We cannot be certain about the structure of the real optimal solution.There are two possibilities.The exact set {(x, t) ∈ Q | | φ(x, t)| = κ} might have positive measure as the computation seems to indicate.On the other hand, this might be a numerical artifact and the measure of the set is zero.In that case, the optimal control should touch the bounds ± 10.Perhaps our numerical resolution is not sufficiently fine to see this.To show that also the first explanation might hold true, we construct an associated example as follows.to obtain φ.For the control, we define ū(x, t) := max(0, (1 − (t − 2.5) 2 ) * ū0 (x)), where ū0 (x) := max(0, 1 − ((x − 25)/10) 2 ).This analytically defined optimal control along with its associated state and adjoint state are shown in Figure 7.This approach has the advantage that we have a well-defined stationary point for ν = 0 and can confirm the convergence against this point.Also in this example, we observe a linear decay of ȳν − ȳ .For the norm of L 2 (Q) this is displayed in Figure 7; for the norm of L ∞ (Q) this looks similar.
Remark.Since ū only takes values in [0, 1] and our approach does not depend on the bounds a and b, ū satisfies the first order necessary optimality conditions for any a ≤ 0 and b ≥ 1.In particular, this works for a < 0 and b > 1 so that the control does not touch the bounds.Let us mention that our example is not completely analytic.The states ȳ and y Q are numerically obtained.We emphasize that the first order necessary optimality conditions are satisfied but that global optimality of ū cannot Downloaded 11/06/15 to 193.144.185.28.Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.phpbe deduced from this.Even the local optimality cannot be ensured.For this purpose, we would have to check second order sufficient optimality conditions.This cannot be done numerically.However, even with the quite high penalization from the sparse term, the computed state ȳ is very close to the desired y Q .Most likely, the consistence of the first order optimality conditions together with the quite small objective value indicate global optimality of ū.

, Lemma 2 . 3 ]
directly leads to the estimate for y u − ȳ in(3.18).The estimate of ϕ u − φ in(3.19) is obtained by considering the difference of the adjoint equations for ϕ u and φ.Notice that the right-hand side of the associated adjoint PDE of this difference is just C Y Q (y u − ȳ), while in the terminal condition the term C Y T (y u (T ) − ȳ(T )) appears.Here, we use that y u the same homogeneous boundary and initial conditions as for w above.Then the estimate (3.23) for w = y u − ū − η u−ū follows as the one for w.Lemma 3.10.There exists a constant M a,b such that ∀u ∈ U ad and ∀v 1 , v 2 ∈ L 2 (Q) the following estimate holds:

Fig. 1 .
Fig. 1.Example 1, natural development y nat of the wave profile for u = 0 (left) and the desired state y Q (right).

Fig. 4 .
Fig. 4. Example 2, natural development y nat of the impulse for u = 0 (left) and the desired state y Q (right).