Optimal control of a class of reaction–diffusion systems

. The optimal control of a system of nonlinear reaction-diﬀusion equations is considered that covers several important equations of mathemat- ical physics. In particular equations are covered that develop traveling wave fronts, spiral waves, scroll rings, or propagating spot solutions. Well-posedness of the system and diﬀerentiability of the control-to-state mapping are proved. Associated optimal control problems with pointwise constraints on the con- trol and the state are discussed. The existence of optimal controls is proved under weaker assumptions than usually expected. Moreover, necessary ﬁrst- order optimality conditions are derived. Several challenging numerical examples are presented that include in particular an application of pointwise state constraints where the latter prevent a moving localized spot from hitting the domain boundary.


Introduction
We consider a class of optimal control problems for a general system of nonlinear reaction-diffusion equations that covers a variety of particular cases with important applications in mathematical physics. Control problems of this type received increasing attention in the recent past, we refer for instance to [16,24,29,33] with respect to control methods of theoretical physics or [3,4,13,12,11] from a more mathematical perspective.
The general system is posed in Q := Ω ×(0, T ), where Ω ⊂ R d is a bounded spatial domain and T > 0 is a finite time. It has the form E ∂y ∂t − D∆y + R(x, t, y) + A(x, t) y = B(x, t) u (1.1) subject to appropriate initial and boundary conditions. Here, y : Q → R n is the state vector function and u : Q → R nc is the control vector function. Moreover, constant diagonal (n, n)-matrices E, D, and matrix-valued functions A, B of suitable dimensions are given, where E has positive diagonal entries and the ones of D are non-negative. The nonlinearity of the system is defined by the function R : Q × R n → R n that we suppose to have a diagonal structure, i.e. the j-th component of the vector function R depends only on (x, t, y j ). The assumptions on A, B, and R are detailed in the next section.
In particular, the system (1.1) includes the following special cases that we recall with increasing order of complexity. Here, we have n = n c = 1, the matrices E = D = B reduce to the real number 1 and A is equal to the real number 0. The system (1.3) is known to exhibit traveling waves as typical solutions. We refer to [28] and to the standard textbooks [14,22,31] on reaction-diffusion equations. The optimal control of (1.3) was investigated in [4] with focus on the numerical analysis and in [16] from a physical point of view; we also refer to various examples in [29].
(ii) The FitzHugh-Nagumo system A standard type of this well-known system is σ y ∂y ∂t − D y ∆y + R(y) + α z = u σ z ∂z ∂t + β z − γ y + δ = 0, (1.4) where σ y , σ z , and D y are positive and α, β, δ, γ are real numbers. This system fits in (1.1) with n = 2, n c = 1; A has the rows (0 α) and (−γ β), B = (1 0) , and R(x, t, y, z) = (R(y), δ) . For spatial dimension d = 1, the system has impulses as typical solutions, while it develops turning spirals or scroll rings for d = 2 and d = 3, respectively. We refer to [14,22,31] for general results on the equation, to [24,33] on applications of control methods in theoretical physics and to [11,12] with respect to the optimal control from a mathematical perspective.
(iv) Coupling of (1.3) with a linear system of PDEs In this more general system, the equation (1.5) is coupled with the linear equations σ j ∂z j ∂t − D j ∆z j + β j z j − γ j y + δ j = 0, j = 1, . . . , m. (1.7) Here, we have almost the same quantities as in (iii), but D has to be defined by D = diag(D y , D 1 , . . . , D m ) with positive numbers D j . This system was discussed in the PhD thesis [23], too.
Other similar reaction-diffusion-equations covered by the general system (1.1) and associated control strategies have been of great interest during recent years. We mention exemplarily [1,20,21,25,27,32].
We will present several numerical examples for the different types of equations listed above. All of them are adopted from the thesis [23]. Our mathematical analysis is developed for the more general system (1.1) that includes all these particular cases for the case of unit matrices E and D.
Our paper contains the following novelties: We extend the analysis on existence, uniqueness, and differentiability of the control-to-state mapping to the general state equation (1.1). Here, we apply a different method of proof than in [23] and avoid the application of analytic semigroups that turns out to be too technical for the general case. We prove the existence of optimal controls under weaker assumptions than usually expected: In the unconstrained case a = −∞ or b = ∞ it seems that the existence of optimal controls cannot be shown. The L 2 -Tikhonov term ensures only the L 2 (Q) nc -boundedness of the set of controls, where the infimum of the problem can be attained. However, with exception of d = 1, the control-to-state mapping G is not continuous from L 2 (Q) nc to the state space. Therefore, this boundedness in L 2 (Q) nc seems to be useless. Nevertheless, we are able to prove the existence of optimal controls. Moreover, we consider problems with pointwise state constraints and prove associated necessary optimality conditions.
The consideration of the state constraints is motivated by an interesting application in theoretical physics, namely the problem of preventing localized moving spots from reaching the boundary of the spatial domain. This issue is discussed in Example 3 of Section 4. Remark 1.1 We will develop our analysis for (n, n)-unit matrices E = D = I d , since this is less technical. If E and D have positive diagonal entries, then the results remain true with E and D substituted for I d in the associated positions. The proofs need only minor modifications. If some of the diagonal entries of D are vanishing as in (1.4) or (1.6), then some equations of (1.1) are ordinary differential equations w.r. to t. In these equations, x plays the role of a parameter and boundary conditions are not given. If they are linear as in (1.6), then the associated vector function z can be represented by the variation-ofconstants formula in terms of y and herafter eliminated. Then a system for y of the form (1.1) is obtained, where D has positive diagonal elements again. If z and y appear nonlinearly in the ordinary differential equations, then the situation is more complicated. This case is not covered by our theory. The matrices E and D will be needed in our numerical examples.

State equation.
We shall consider the following system of parabolic partial differential equations in Ω, where Ω is an open bounded domain in R d with d ∈ {1, 2, 3}, Γ is the boundary of Ω that we assume to be Lipschitz, ν is the unit outward normal vector Given a control u : Q −→ R nc , we look for the corresponding solution y u : Q −→ R n of (2.1). The analysis of this equation will be done below under the following assumptions on R. We assume that every function R j is of class C 1 with respect to the last variable and satisfies for 1 ≤ j ≤ n R j (·, 0) ∈ Lp(0, T ; Lq(Ω)) withp,q ∈ [2, ∞] and Along this paper, we take A L ∞ (Q,R n×n ) := ess sup (x,t)∈Q A(x, t) , where A(x, t) denotes the matrix norm induced by the Euclidean norm in R n . Analogously, we define B L ∞ (Q,R n×nc ) .
Proof For every M > 0, we set R M (x, t, y) = R(x, t, Proj [−M,+M ] n (y)). Now, given w ∈ L 2 (Q) n , we consider the system in Ω. (2.6) For every η > A L ∞ (Q,R n×n ) , the operator −∆+A+ηI is coercive in H 1 (Ω) n , hence there exists a unique solution y w ∈ L 2 (0, T ; H 1 (Ω)) n of (2.6); see, for instance, [30, page 112]. Further, from the partial differential equation we get that ∂ t y w ∈ L 2 (0, T ; H 1 (Ω) * ) n as well, hence y w ∈ W (0, T ) n . For fixed M > 0, the right hand side of the partial differential equation (2.6) is bounded in L 2 (Q), hence the Schauder fixed point theorem applied to the mapping w ∈ L 2 (Q) n → y w ∈ L 2 (Q) n along with the compactness of the embedding W (0, T ) ⊂ L 2 (Q), implies the existence of a fixed point y M ∈ L 2 (0, T ; H 1 (Ω)) n that satisfies the equation in Ω. (2.7) Next, we derive some estimates for y w . Let us set η = |C R |+ A L ∞ (Q,R n×n ) +1 and take 0 < T ≤ T . Multiplying (2.7) by e −2ηt y M and integrating in Q = Ω × (0, T ), we get For the first term of the left hand side of the above identity we have Moreover, applying the mean value theorem and using (2.3) we obtain We also have that y M A(x, t)y M ≥ − A L ∞ (Q,R n×n ) |y M | 2 . Using these two inequalities, inserting (2.9) in (2.8), and taking into account the definition of η, we infer by the Schwarz inequality that Applying Young's inequality and multiplying the resulting inequality by 2 we get Therefore, dropping the second term, using e −2ηt ≥ e −2ηT with 0 < T ≤ T , and multiplying the inequality by e 2ηT , we conclude (2.10) Since T was arbitrary, we deduce that {y M } M is uniformly bounded in the space L ∞ (0, T ; L 2 (Ω)) n ∩ L 2 (0, T ; H 1 (Ω)) n . It remains to prove the boundedness of y M in L ∞ (Q) n . To this end, consider the functions f j : Q × R −→ R and g j , 1 ≤ j ≤ n, defined by From (2.7), it is obvious that y M,j satisfies the equation y M,j (x, 0) = y 0,j (x) in Ω. (2.11) Due to the regularity of g j and the fact that f j is monotone non decreasing, we can use the methods of [15,§III.7] to infer the existence of a constant C j independent of M such that y M,j L ∞ (Q) ≤ C j y 0;j L ∞ (Ω) + (Bu) j Lp(0,T ;Lq(Ω)) + R j (·, ·, 0) Lp(0,T ;Lq(Ω)) + (1 + |C R |)y M,j − (Ay M ) j L ∞ (0,T ;L 2 (Q)) ≤ C j y 0 L ∞ (Ω) n + u Lp(0,T ;Lq(Ω)) nc + R(·, ·, 0) Lp(0,T ;Lq(Ω)) n + y M L ∞ (0,T ;L 2 (Ω)) n . (2.12) Combining (2.10) and (2.12), we obtain the boundedness in L ∞ (Q) n of every y M by a constant C independent of M . Hence, taking M > C, we see that R M (x, t, y M (x, t)) = R(x, t, y(x, t)). This implies that y M is solution of (2.1). Moreover, (2.5) follows from (2.10) and (2.12). Now, we prove the uniqueness of a solution. Let us assume that y 1 , y 2 ∈ Y are two solutions of (2.1). We set y = y 2 − y 1 . Subtracting the equations satisfied by these functions, we obtain y(x, 0) = 0 in Ω.
and arguing as above, we get This shows that y = 0. Finally, if y 0 ∈ C(Ω) n , then we can apply the results of [15,§III.7] to deduce that y ∈ C(Q) n .

Remark 2.2
Of course the previous theorem remains valid if we consider a more general elliptic operator with L ∞ (Q) coefficients. Moreover, the space Y = W (0, T ) n ∩ L ∞ (Q) n introduced previously can be substituted by Y = W (0, T ) n ∩ C(Q) n , provided that y 0 ∈ C(Ω) n . Then, Theorem 2.1 and the next results are true under this new definition of Y . We also mention that the diagonal structure of R could be avoided under some additional assumptions. For instance, Theorem 2.1 is valid assuming that R is of polynomial order with respect to y and certain kind of monotonicity is satisfied: there exist constants C i > 0, i = 1, 2, 3, such that for almost all (x, t) ∈ Q and all y, The difficulty of dealing with more general functions is due to the lack of L ∞ -estimates, which are crucial in the above analysis.

Remark 2.3
As usual, assuming more regularity of y 0 and Ω, we get extra regularity of y: If y 0 ∈ H 1 (Ω) n , then y belongs to H 1 (Q) n , and an associated estimate is valid. If, in addition to this regularity of y 0 , Γ is of class C 1,1 or Ω is convex, then y ∈ H 2,1 (Q) n .
Due to Theorem 2.1, we can define a mapping G : Lp(0, T ; Lq(Ω)) nc −→ Y by G(u) = y u , where y u is the solution of (2.1) associated with u. The next theorem states the differentiability of this mapping.

Theorem 2.2
The mapping G is of class C 1 and the derivatives z v = G (u)v are the solutions of the system Proof We define the space By the definition of V and assumption (2.2)-(2.4), it is clear that F is well defined and is of class C 1 . Moreover, the identity F(G(u), u) = 0 holds for all u ∈ Lp(0, T ; Lq(Ω)) n . We want to apply the implicit function theorem to deduce the differentiability of G and to obtain z v = G (u)v as the solution of (2.13). To this end, we only need to prove that the linear operator has a unique solution in V and the solution z depends continuously on the data (f, g, z 0 ). This existence, uniqueness and continuity follows from Theorem 2.1. Therefore, G is of class C 1 . Finally, (2.13) follows from the identity The additional regularity of z v is a consequence of the initial condition z v (0) = 0, because z v (0) is continuous and belongs to H 1 (Ω).

Optimal control problem.
In this section, associated with the state equation (2.1), we consider the following control problem (P) min n Ω , and the admissible set Of course, the relations a ≤ u(x, t) ≤ b are to be understood componentwise. Obviously, this assumption implies that U a,b = ∅.
As an immediate consequence of Theorem 2.2, we deduce the following corollary just by an application of the chain rule. (2.14) where ϕ u ∈ L 2 (0, T ; H 1 (Ω)) n is the unique solution of the adjoint state equation Using this result, we can derive the optimality condition for a local solution to (P). Theorem 2.3 Assume that a, b ∈ R nc . Then there exists at least one solution of the control problem (P). Any local solutionū satisfies the variational inequality whereφ = ϕū is the adjoint state associated withū.
Since U a,b was assumed to be bounded in L ∞ (Q) nc , the proof of the existence of an optimal control is standard by taking a minimizing sequence that converges weakly * in L ∞ (Q) to someū ∈ U a,b . Moreover, we mention that u k * u weakly * in L ∞ (Q) implies that y u k → y u strongly in L 2 (Q). This follows from Theorem 2.1 and the Aubin-Lions Theorem. Finally, we takep =q = +∞ and apply Corollary 2.1 to deduce the optimality conditions.
Let us now analyze the case where U a,b is not bounded, which in particular includes the case without control constraints. Here, we slightly extend ideas of [8]. First we observe that Theorem 2.1 cannot be applied to deduce the existence of a solution of the state equation for arbitrary elements in U a,b because the L ∞ -estimates for the states fail. We need to assume more regularity for the elements of U a,b , the assumed L 2 (Q) n C regularity in the definition of U a,b is not enough. Therefore, we have to work with a control space of type Lp(0, T ; Lq(Ω)) nc withp,q ∈ [2, +∞] and 1 p + d 2q < 1. The reader can easily identify the difficulty of proving a solution of the corresponding control problem because of the lack of coercivity of the cost functional in these spaces. Notice that the Tikhonov regularization term yields only coercivity in the space L 2 (Q) nc but not in the control space introduced above. Let us redefine Theorem 2.4 Let λ be strictly positive. Then, the optimal control problem (P) has at least one solutionū. Any local solution in this space satisfies the variational inequality (2.16).
Proof For every 1 ≤ j ≤ n c , we select an element ξ j ∈ R such that a j < ξ j < b j . We set ξ = (ξ j ) nc j=1 ∈ R nc and M 0 = |ξ| |Ω|, where |Ω| denotes the Lebesgue measure of Ω. For every M ≥ M 0 , we define Then we have that u 0 ∈ U M ∀M ≥ M 0 and, consequently, U M = ∅. According to Theorem 2.1 we know that (2.1) has a unique solution y u ∈ Y for every u ∈ L ∞ (0, T ; L 2 (Ω)) nc . We formulate the following optimal control problem For every M ≥ M 0 , this problem has at least one solution u M . Indeed, any minimizing sequence {u k } is bounded in L ∞ (0, T ; L 2 (Ω)) nc , hence the sequence of corresponding states {y k } is bounded in Y . Moreover, from (2.1) we also obtain the boundedness of {y k } in W (0, T ) n . Then we can select subsequences, denoted in the same way, such that u k * u M in L ∞ (0, T ; L 2 (Ω)) nc and y k → y M strongly in L 2 (Q) n due to the compactness of the embedding W (0, T ) ⊂ L 2 (Q). Now, it is easy to pass to the limit in the state equation and to deduce that y M is the state associated with u M . Finally, we have that Therefore, u M is a solution of (P M ). Then from (2.14) we infer that where ϕ M is the adjoint state associated with u M . This relation implies that where the projection is taken in the L 2 (Q) nc norm. Let us prove that there Since y M ∈ L ∞ (Q) n , we can multiply the state equations (2.1) by e −2ηt y M with η = 1 + |C R | + A L ∞ (Q,R n×n ) and argue as in the proof of Theorem 2.1, without any truncation of R (replace R M by R in that proof), to obtain (2.10). Then, using (2.18), we deduce from (2.10) the existence of a constant C 1 such that y M L ∞ (0,T ;L 2 (Ω)) n ≤ C 1 ∀M ≥ M 0 .
Since y M ∈ W (0, T ) n ⊂ C(0, T ; L 2 (Ω)) n , we infer from the above inequality Now, from the adjoint state equation satisfied by ϕ M , it follows the existence of a constant C 2 such that Combining this inequality and (2.19), we conclude that Notice that I a or I b or both can be empty. In any of these cases, the corresponding maximum is taken as 0. If we defineũ ∈ U a,b bỹ This, along with (2.17), implies In this way, we have proved that u M L ∞ (0,T ;L 2 (Ω)) nc ≤M ∀M ≥M . Now, we show that, for every M ≥M , u M is a solution of (P). Let us take u in the space L ∞ (0, T ; L 2 (Ω)) nc and set M = u L ∞ (0,T ;L 2 (Ω)) nc . If Finally, an application of the Corollary 2.1 withp = +∞ andq = 2 leads to (2.16).

Remark 2.4
Let us notice that the optimal controls, whose existence is proved in Theorems 2.3 and 2.4, belong to C([0, T ]; L 2 (Ω)) nc ∩ L 2 (0, T ; H 1 (Ω)) nc assuming that b ij ∈ C([0, T ]; W 1,∞ (Ω)) for every entry of B. Indeed, this regularity is well known for the adjoint state. Since the associated optimal control is taken as a projection on U a,b , this regularity is transferred to the control. The projection formula in the first case follows from (2.6) and it is given in (2.21) in the second case.
3 State constrained control problem

Optimal control problem and existence of a solution
In this section we analyze the following state constrained control problem We impose the following assumptions on the data of the control problem.
(A1) The control u is related to the state y u through the system (2.1). We assume that (2.2)-(2.4) hold and y 0 ∈ C(Ω). We set Y = L 2 (0, T ; H 1 (Ω)) n ∩ C(Q) n . Notice that according to Theorem 2.1, the states y u ∈ Y due to the continuity of y 0 ; see Remark 2.2.
(A2) The cost functional is given as in §2.2 with λ ≥ 0 and the same regularity for matrix functions C Q , C Ω , and functions y Q , and y Ω .
(A3) We assume that a, b ∈ R nc with a j < b j for 1 ≤ j ≤ n c .
(A4) K is a compact subset ofQ and the function g : K × R n −→ R is continuous, together with its partial derivatives ∂g/∂y j ∈ C(K × R n ), j = 1, . . . , n. We also assume that g 1 , g 2 : K −→ R are continuous functions with Under the above assumptions, we have that U a,b is a closed, convex and bounded subset of L ∞ (Q) nc and the mapping G : L ∞ (Q) nc −→ Y , defined as in §2.1 by G(u) = y u , is of class C 1 . The derivative z v = G (u)v ∈ Y is the solution of (2.13). From Corollary 2.1 we know that J : L ∞ (Q) nc −→ R is of class C 1 and its derivative J (u)v is given by (2.14).
By using the above assumptions we can prove the following theorem on existence of an optimal control. Proof The proof follows classical arguments. Indeed, since U ad = ∅ there exists a minimizing sequence {u k } ∞ k=1 ⊂ U ad of (PS). This sequence is bounded in L ∞ (Q) nc and the associated states {y k } ∞ k=1 are bounded in Y . Additionally, this boundedness along with the partial differential equation (2.1) implies that Therefore, we can take subsequences, denoted in the same way, such that u k * ū in L ∞ (Q) nc withū ∈ U a,b , and y k →ȳ strongly in L 2 (Q) n and y k (x, t) →ȳ(x, t) for almost all (x, t) ∈ Q. Now, it is easy to pass to the limit in the state equation (2.1) as k → ∞ and to prove thatȳ ∈ Y is the state associated toū. Moreover, using the continuity of g we can pass to the limit in the inequality g 1 (x, t) ≤ g(x, t, y k (x, t)) ≤ g 2 (x, t) to obtain g 1 (x, t) ≤ g(x, t,ȳ(x, t)) ≤ g 2 (x, t) for almost all (x, t) ∈ K. But the continuity of the functionsȳ, g 1 , g 2 and g implies that the above inequality holds for all (x, t) ∈ K. Hence, we have thatū ∈ U ad . Finally it is obvious that which proves thatū is solution of (PS).

Optimality system
Hereafter,ū will denote a local minimum of (PS) with associated stateȳ. In order to get the optimality conditions satisfied byū in a qualified form we make the following linearized Slater condition.

(3.2)
This section is devoted to the proof of the following theorem.
Before proving this theorem, we recall what a solution of (3.4) is. We will do it in an abstract framework. We consider a vector of real and regular Borel measures µ ∈ M(Q) n such that |µ j |(Ω × {0}) = 0 for 1 ≤ j ≤ n. We decompose µ = µ Q + µ Σ + µ Ω by taking the restrictions of µ to Q, Σ and Ω × {T }, respectively. Now, we consider the system where D ∈ L ∞ (Q, R n×n ).
The reader is referred to [6] or [10] for the proof of this theorem in the case of a scalar equation. The arguments are identical for the above system. In the proof of theorem some regularity results for the adjoint system to (3.7) are required, they are deduced from Theorem 2.1 just taking R ≡ 0. From the regularity of φ established in the above theorem and using a density argument, it is easy to prove that for every φ ∈ [H 1 (Q) ∩ C(Q)] n such that ∂ ν φ ∈ L ∞ (Σ) n and φ(·, 0) = 0 in Ω. Now, we can apply Theorem 3.3 to the system (3.4) with taking into account thatȳ ∈ C(Q) n and using (2.4). For the right hand side of the equations we consider the measures µ Ω (C) = Theorem 3.4 Let U and Z be two Banach spaces and K ⊂ U and C ⊂ Z two convex subsets, C having a nonempty interior. Letū ∈ K be a solution of the optimization problem: where J : U −→ (−∞, +∞] and F : U −→ Z are two Gâteaux differentiable mappings atū. Then there exist a real numberμ 0 ≥ 0 and an element µ ∈ Z * such thatμ Moreoverμ 0 can be taken equal to 1 if the following linearized Slater condition is satisfied: Now, Theorem 3.3 follows by taking U = L ∞ (Q) nc , K = U a,b , Z = C(K), C = Y g , J is the cost functional of (PS) and F (u) = g(·, ·, G(u)). Then, we have that (3.9) and (3.10) coincide with (3.3) and (3.5), respectively. Let us prove that (3.11) is the same as (3.6). To this end, we first introduceφ as the solution of the system (3.4). We observe that the chain rule along with the expression of z v = G (ū)v provided in (2.13) leads to From here we get, once again with the chain rule, where µ = µ Q +µ Σ +µ Ω is the measure above introduced. Now, we notice that Theorem 2.2 implies z v ∈ Φ and ∂ ν z v = 0. Then from Definition 3.1 applied to the system (3.4) with φ = z v , and recalling that D = ∂g ∂y + A we obtain with (2.13) The last two identities lead to Hence, taking v = u −ū with u ∈ U a,b , we conclude that (3.11) and (3.6) are identical. It is obvious that (3.1) and (3.12) are equal, hence the possibility of takingμ 0 = 1 under the Slater assumption follows from Theorem 3.4. After having proved Theorem 3.2, let us draw some conclusion, namely some information onμ that follows from (3.5).

A regularity result for local solutions
As in the previous section,ū will denote a local minimum of (PS) with associate state and adjoint stateȳ andφ, respectively. In this section we impose the following additional assumption on the problem (PS).
Under this assumption, the well known regularityū ∈ L r (0, T ; W 1,s (Ω)) nc for all r, s ∈ [1, 2) with 2 r + d s > d + 1 follows from the projection formulā which is well known to be deduced from (3.6). However, this projection formula leads to higher regularity, namelyū ∈ L 2 (0, T ; H 1 (Ω)) nc . The next lemma, proved in the appendix, is the key tool to establish this regularity.
Then, from Lemma 3.1 we know that ϕ M ∈ L 2 (0, T ; H 1 (Ω)) n . Due to the regularity of the functions b j we also have that b j ϕ M,j ∈ L 2 (0, T ; H 1 (Ω)) for every 1 ≤ j ≤ n c . Now, from (3.14) we havē This implies thatū j ∈ L 2 (0, T ; H 1 (Ω)) n as well for 1 ≤ j ≤ n c . It remains to prove that both projections coincide. This is obviously the case if |φ j (x, t)| ≤ M . If |φ j (x, t)| > M , the reader can easily confirm the following facts In either case the equality of the projections follows.

Examples
As applications, we consider systems of equations in two-dimensional spatial domains (d = 2) that develop spiral waves or moving localized spots as solutions. Spiral waves appear for the FitzHugh-Nagumo equations, a system of 2 equations, while localized spots arise for a system of 3 equations. In all the examples, the aim is to move the appearing state function in a prescribed way. All examples are numerically very challenging but show, on the other hand, the geometrical beauty of solutions to the selected reaction-diffusion equations.
As pointed out in Remark 1.1, the system (1.4) does not directly fit to (1.1), since the second diagonal element of D is zero. However, our theory remains true with obvious modifications. For the necessary optimality conditions and the adjoint equation we refer to [12], where this example was not considered. Here, our control task is to translate a naturally developed spiral wave pattern along a given circle. By a standard method that is explained in [23, p. 48], a rotation of the states is triggered: Spiral waves y 0 and z 0 are computed in Ω as initial states for the system (1.4); they are depicted in Fig.1.
The area of the center of the spiral, the so-called "core", is located around the position (3/4 L Ω , 1/2 L Ω ) = (56.25, 37.5). The desired trajectory y Q equals the uncontrolled natural state y evolving from (y 0 , z 0 ) that is translated in counter-clockwise direction along a circular shape with radius 1/4 L Ω around the center of the domain (1/2 L Ω , 1/2 L Ω ). Due to the Neumann boundary conditions, this is delicate issue.
However, the position of a spiral pattern is basically determined by the location of its core. Translating the core, the arm of the spiral follows accordingly with some delay. Hence, we consider the desired trajectory only in a circle-shaped area of radius 15 around the desired center X(t) of the spiral, given by X(t) := (1/2 L Ω + 1/4 L Ω cos(2πt/T ), 1/2 L Ω + 1/4 L Ω sin(2πt/T )). We set Fig. 2 displays C Q at t = 0 as well as the product C Q y Q for some times t. The remaining parameters of the optimal control problem (P) are set to u a = −1, u b = 1, and C Ω = y Ω ≡ 0.
As optimization algorithm, a projected gradient-method with non-linear CG-step was chosen. Due to the large time-horizon with 3001 time-steps, and the circumstance that an entire desired trajectory y Q is given, We employed Model Predictive Control with 4 time-steps of length τ in each sub-problem. This means that we solve a sequence of (discretized) short time optimal control problems in a time horizon of length 4τ starting with the interval [0, 4τ ]. From the 4 values computed for the discretized optimal control in this short time interval, we keep the first value for the final suboptimal control. Next we move the time horizon one time step to the right, compute the next optimal control for this shifted time horizon and keep again its first value. After having solved a finite number of small optimal control problems, we arrive at a suboptimal control on [0, T ]. The short time control problems were solved by a nonlinear CG-step.
Moreover, since the chosen discretization of 101 × 101 grid points in space still leads to fairly high computation times, only a semi-implicit Euler-scheme for solving the discrete systems was applied. Yet, 13.42 hours for computing the suboptimal controlū was quite large. Fig. 3 illustrates the computed (sub-)optimal controlū with associated activator stateȳ at several times t. As shown, the control task is satisfied; the suboptimal value of the objective functional was f (ū,ȳ) = 9.899 · 10 −3 .
One should expect that the optimal control concentrates at the support of the function C Q . However, the highest amplitudes of the control appear at the boundary of the circumcircle of the support. The reason is that the profile of the given desired trajectory is the one of an uncontrolled "standing" spiral wave. A translation of the pattern naturally leads to some deformations of the profile and the control aims to suppress those deformations where C Q is positive.
Similarly to Example 1, the control task is to translate a naturally developed pattern along a circle-shaped curve in Ω in counter-clockwise direction. For this purpose, we proceed as follows. First, we construct a natural developed spot profile as follows: We take as auxiliary initial states for u ≡ 0 with initial data (ỹ 0 ,z 1 0 ,z 2 0 ) subject to periodic boundary conditions. Eventually, after a finite time, a stable spot profile is generated. As soon as the center of mass of the pattern is in the center of the domain Ω, we replace the boundary conditions by homogeneous Neumann-type. The violation of those conditions is negligible and disappears after a few further time-steps. In fact, we even let enough time pass by to have the center of mass of the spot profile in the activator state y situated in (3/4 L Ω , 1/2 L Ω ) = (0.45, 0.3). This state is taken as initial state, cf. Fig. 4. Analogously to Example 1, we define the support function where X(t) is defined as in Example 1. The desired state y Q is defined by and we fix C Ω = y Ω ≡ 0, T = 250.3, as well as u a = −1 and u b = 1.
The problem is solved numerically in the same way as for Example 1, this time with 2504 time-steps. The computed suboptimal objective value is f (ū,ȳ) = 2.49 · 10 −6 . The associated (sub-)optimal computed controlū is shown in Fig. 5 along with the stateȳ for various times t.
An interesting property of the spot can be observed in this example. The pattern is oriented and hence, a natural translation of the spot in positive x 1direction would occur in the uncontrolled case. However, the desired trajectory is determined by a spot profile in which orientation and natural movement do not comply. They only comply in t ≈ 3/4 T which causes the computed control to have much lower amplitudes during these times. In all other times, the control does not only force the spot to change its position but also to keep its "non-complying" profile unchanged.
For control tasks of this type, namely translating a natural developed pattern in a reaction-diffusion equation without changing the profile, we also refer to [19,17,18,24].
To complete the setup of the optimal control problem, we set C Ω = y Ω = C Q = y Q ≡ 0, u a = −1, and u b = 1.
The numerical solution of this example is based on a finite difference discretization with h = 1/200, τ = 0.1 as step sizes in space and time and is performed as in the preceding 2 examples.
Again, we applied Model Predictive Control with 4 time-steps and a nonlinear CG optimization method for solving each subproblem. The state constraints have been included in an associated penalty term.