Control of traveling localized spots

Traveling localized spots represent an important class of self-organized two-dimensional patterns in reaction-diffusion systems. We study open-loop control intended to guide a stable spot along a desired trajectory with desired velocity. Simultaneously, the spot's concentration profile does not change under control. For a given protocol of motion, we first express the control signal analytically in terms of the Goldstone modes and the propagation velocity of the uncontrolled spot. Thus, detailed information about the underlying nonlinear reaction kinetics is unnecessary. Then, we confirm the optimality of this solution by demonstrating numerically its equivalence to the solution of a regularized, optimal control problem. To solve the latter, the analytical expressions for the control are excellent initial guesses speeding-up substantially the otherwise time-consuming calculations.


Introduction
Localized spots, sometimes referred to as auto-solitons [22], dissipative solitons [39], or bumps [24], are a subclass of traveling waves that spontaneously evolve in two-dimensional (2D) dissipative nonlinear systems driven far from thermodynamic equilibrium. In a co-moving reference frame, spots are stationary solutions to coupled nonlinear partial differential equations (PDE), such as reaction-diffusion (RD) or neural field equations, for example. The characteristic length and time scales of the spots, i.e., their wave profile, propagation velocity, etc., are selected by the experimental conditions or the parameters of the model.
Experimentally, localized spots have been observed as current filaments in gas-discharge [40], as bright intensity spots in nonlinear optics and laser physics [1], as well as moving localized regions of increased concentration in chemical reactions [51] or coverage of adsorbed species in heterogeneous catalysis [53]. Further examples include temperature spots in fixed-bed catalytic reactors [52], actin conformation in dictyostelium discoideum [25], neural activity in head-direction cells [47], vegetation patterns [15], and many others.
Although control of self-organized patterns attracted considerable attention over the last decades, compare [33,50] and references therein, it is still a challenging problem in applied nonlinear dynamics. In the context of localized spots, an illustrative applicative example are temperature spots in catalytic reactors. Transient and sustained moving hot-spot activity spontaneously formed near the wall of catalytic packed-bead and flow reversal reactors, respectively, may pose severe safety hazard problems. Another example is the control of localized neural activity including so-called bump solutions to neural field equations describing ensembles of synaptically coupled neurons [11,24,56]. Catalytic reactors and neural networks are only two fields where spontaneously emerging traveling spots may play a role. Other areas include gas-discharge, polymerization fronts, and solid fuel combustion.
Often, one distinguishes between open-loop, closed-loop, and optimal control. Open-loop control is independent of the instantaneous state of the system. As a consequence, it is inherently susceptible to perturbations in the initial conditions as well as to parameter uncertainty. Thus, detailed knowledge of the system's dynamics and in-depth stability analysis are pre-conditions for reliable open-loop schemes. Typical examples of openloop control are space-time dependent external forcing [10,58] or control by imposed geometric constraints [19,57]. On the other hand, in closed-loop or feedback control, the controlled state is permanently monitored to adjust the control signal accordingly [23,30,38]. Particularly, time-delayed feedback can induce pattern forming instabilities in addition to the pattern to be controlled [16].
Optimal control reformulates control problems in terms of the minimization of a cost functional [5,20]. Often, as in our paper, the cost functional measures the distance in function space between a desired target state and the actual controlled state of the system. If a control signal is the unique solution to an optimal control problem, then no other control, be it open-or closed-loop, will be able to enforce a controlled state closer to the target state. For nonlinear models, the uniqueness of optimal controls is a difficult and mathematically widely open issue. If an optimal control satisfies a second-order sufficient optimality condition, then it is at least locally unique, i.e., it delivers an isolated minimum; we refer to [2,9]. Optimal control of self-organized patterns requires complete knowledge of the PDEs governing the system's evolution in time and space. Numerical solutions to optimal control of PDEs often base on computationally expensive iterative algorithms. The convergence to the target state can sensitively depend on an appropriate initial guess for the control signal.
For traveling wave patterns, a general control task is position control aimed at guiding the pattern according to a given protocol of motion (POM), i.e., moving it with desired velocity along a desired trajectory through a spatial domain. In some technical applications like catalytic reactors, it is necessary to avoid the collision of high-temperature spots with the reactor walls or their pinning at heterogeneities of the catalyst's support to maintain operational safety [52]. Another example of open-loop position control is the enhancement of the CO 2 production rate during the low-pressure catalytic oxidation of CO on Pt(110) single crystal surfaces by dragging reaction pulses and fronts using a focused laser beam with a speed differing from their natural propagation velocity in the absence of control [41,54]. In a photosensitive Belousov-Zhabotinsky (BZ) medium, periodic variation of the applied light intensity can force a spiral wave tip to describe a wide range of hypocycloidal and epicycloidal trajectories [45,46]. In optical bistable media like dye-doped liquid crystals and Kerr cavities, interface dynamics can be controlled by spatially inhomogeneous forcing [37]. Position control of traveling wave patterns can be tackled by feedback control as well. For example, the spiral wave core in a photosensitive BZ medium was steered around obstacles using feedback signals obtained from wave activity measured at a point detector, from tangential crossing of wavefronts with detector lines, or from a spatially extended control domain [44,58,59]. Two feedback loops were used to stabilize and guide unstable traveling wave segments along pre-given trajectories [43]. Furthermore, feedback-mediated control loops were employed to prevent transversal instabilities in reaction-diffusion waves [48].
Recently, we proposed an open-loop control method that acts solely via the Goldstone modes of wave patterns [26] and, therefore, can be referred to as Goldstone mode control. The method provides analytical expressions for the amplitude of the control signal to be applied for a given POM. We demonstrated that this control is able to accelerate or decelerate 1D traveling front and pulse solutions to RD equations [4,26,27] without changing their spatial profile. The stability of the control loop with respect to small changes in the initial conditions was discussed in [28]. Goldstone mode control also applies to move the core of a spiral wave at desired velocity along a pre-given trajectory through a 2D spatial domain, or to shape iso-concentration lines of 2D traveling pulses [29]. Interestingly enough, the control turned out to be equivalent to the solution of an appropriately formulated optimal control problem [26,42].
In this paper, we extend Goldstone mode control to spatially localized moving spots. We introduce a threecomponent RD model supporting stable traveling spot solutions in Section 2 and derive analytical expressions for position and orientation control of traveling spots in the fully-actuated case in Section 3.1. The corresponding optimal control problem with an objective functional involving a Tikhonov regularization term is formulated explicitly in Section 3.2. Here, we discuss the relation between Goldstone mode control derived in Section 3.1 and the solutions to the optimal control problem. In Section 4, after a brief description of the numerical methods being used, we discuss examples for fully-actuated position and orientation control of spots in Sections 4.1 and 4.2, respectively, as well as for under-actuated position control by a single control signal 4.3. Summarizing, in Section 5 we conclude that our control method will be applicable if symmetry induced Goldstone modes exist as a consequence of translational and rotational invariance of the underlying evolution equation. Additionally, in the spectrum of the linear stability operator of the uncontrolled stable solution a sufficiently large gap should exist between the symmetry-induced neutral eigenvalues on the imaginary axis and the remaining eigenvalues with negative real part.

Three-component spot model
Throughout this work, we consider the following three-component RD system exhibiting immobile and traveling stable spot solutions in 2D [13,39,49] Here, ∆ = ∂ 2 x + ∂ 2 y represents the Laplacian in Cartesian coordinates, r is the position vector in the spatial domain Ω, r = (x, y) T ∈ Ω ⊂ R 2 , and t indicates time. D u , D v , and D w denote the diffusion coefficients of components u, v, and w while τ and θ set the time scales for the v and w kinetics, respectively. Beside spots, the model (2.1) is capable to support peanut patterns [36], breathing solitons [17], and jumping oscillons [55], for example.
Model (2.1) is an extension of the well-known FitzHugh-Nagumo model [14]. It was first introduced by Purwins and co-workers to describe the dynamics of current filaments in planar gas-discharge [39]. In this context, activator u and inhibitor v represent the current density and the voltage drop over the high-ohmic electrode, respectively. The second inhibitor w is linked to the surface charge, and the additive parameter κ 1 is related to the supply voltage.
In a more general interpretation, the model (2.1) represents a RD system of activator-inhibitor type with cubic nonlinearity. Two inhibiting components couple linearly with linear kinetics into the equation of the activator component. In a propagating spot, the inhibitor v is delayed with respect to the activator u, i.e., the center of the u and v concentration fields do not lay on top of each other. In contrast, because of its small time constant θ and large diffusion coefficient, inhibitor w quickly follows the activator distribution and surrounds it entirely. This second inhibitor stabilizes the moving spot. Altogether, (2.1) define a minimal model for traveling 2-D spot solutions.

Analytical expression for control amplitudes in position control
Let us consider a controlled RD system according to (3.1a) Here, U(r, t) = (u 1 (r, t), . . . , u n (r, t)) T is the vector of n ∈ N state components defined in the two-dimensional spatial domain Ω ⊂ R 2 with r = (x, y) T . Assuming an isotropic medium, the n × n matrix of diffusion coefficients D is diagonal and constant, D = diag(D 1 , . . . , D n ). The vector R(U) = (R 1 (U), . . . , R n (U)) T describes the reaction kinetics of the components. In general, R i (U) are nonlinear functions of the state. For the RD system (2.1), U, D, and R are given by and appropriate boundary conditions. We consider a rectangular domain Ω = (x a , x b ] × (y a , y b ] with periodic boundary conditions such that U as well as its derivatives in the direction normal to the boundary are periodic, (3.1c) The space-time dependent control signals f (r, t) = (f 1 (r, t), . . . , f m (r, t)) T , m ∈ N, on the right hand side of (3.1a) are assumed to act for all times t everywhere within Ω. The constant n × m matrix B determines which components are directly affected by the control signals. A system with strictly less independent control signals than components, m < n, is underactuated. For m = n and B invertible, the system is fully actuated. In what follows, we focus on fully actuated systems and set B equal to the identity matrix 1. The limiting case of single component control, i.e., Bf (r, t) ∝ (f 1 (r, t), 0, . . . , 0) T , we consider in subsection 4.3.
The partial differential equations (3.1a) describe the evolution of the components U(r, t) in the presence of spatio-temporal perturbations f (r, t) that break the translation and rotation invariance of the unperturbed equations. In this interpretation, the response of the unperturbed solution to a given small input f can be calculated perturbatively, see references [3,34,35], and Appendix C.2.
In this paper, however, following [26] we perceive (3.1a) for given desired spot dynamics as a conditional equation for the perturbations which now are considered as control inputs. The goal of the control f is to enforce a state U to follow a given desired distribution U d (r, t) = (u 1,d (r, t), . . . , u n,d (r, t)) T as closely as possible everywhere in the spatial domain Ω and for all times 0 ≤ t ≤ T . We call a desired distribution U d exactly realizable if there exists a control f such that the controlled state U equals U d everywhere in the space-time Inserting U d for U in (3.1a) yields for the control For U d to be exactly realizable, three more conditions must be satisfied: First, the initial condition for the controlled state, (3.1b), must coincide with the initial state of the desired distribution, U(r, 0) = U d (r, 0). Second, all boundary conditions for the desired distribution U d have to comply with the boundary conditions for U, (3.1c). Third, U d must be sufficiently smooth in the space-time cylinder Q = Ω × [0, T ] such that the derivatives ∂ t U d and ∆U d are continuous. The construction of such desired distributions is difficult and hence this method of control is not practicable. It is more convenient to control certain patterns that are asymptotic solutions of the uncontrolled system (3.1) for t → ∞ in the space R 2 , for instance spot solutions. Therefore, we formulate the control goal for spot solutions to the uncontrolled RD equations (3.1). These solutions propagate with constant velocity v 0 = (v x 0 , v y 0 ) T and wave profile U c through the spatial domain. In  Table 1. In both panels, identical domain size and colormap are used.
where, ∇ ξ = (∂ ξx , ∂ ξy ) T and ∆ ξ = ∂ 2 ξx + ∂ 2 ξy denote the component-wise gradient and Laplacian, respectively. We emphasize that resting localized spots, v 0 = 0, are rotationally-symmetric solutions while traveling localized spots are axis-symmetric with the symmetry axis directed tangentially to the trajectory of motion, cf. Figure 1a and b, respectively. We characterize the current position of a spot by the x-and y-coordinates of the maximum value of the activator concentration along its symmetry axis at a given time, Φ(t) = (Φ x (t), Φ y (t)) T , and its orientation by the angle Φ ϕ (t) between the spot's symmetry axis and the x-axis, compare Figure A.1.
Remarkably, any reference to the nonlinear functions R drops out from the result (3.6). This is of great advantage in all applications where the details of the underlying reaction kinetics R are largely unknown or difficult to identify. Once propagation velocity v 0 and wave profile U c of the uncontrolled spot are measured with an accuracy sufficient to calculate the Goldstone modes ∂ ξx U c , ∂ ξy U c , and ∂ ϕ U c , the control signals can be computed in advance for the complete time interval [0, T ]. Consequently, in contrast to feedback control, a continuous recording of the system is not required.
One notices that (3.6) equals the sum of Goldstone modes with time-dependent prefactors, f Gold (r, t) = P 1 (t) ∂ ξx U c (ξ) + P 2 (t) ∂ ξy U c (ξ) + P 3 (t) ∂ ϕ U c (ξ). The Goldstone modes are the right eigenvectors to the linear stability operator L of (3.3) to the eigenvalue zero. They are associated with the translational and rotational invariance of equation (3.1a) in R 2 for f (r, t) = 0. Clearly, the prefactors' magnitudes are proportional to the difference between the intrinsic velocity, v 0 , and the current prescribed spot velocity projected onto the x-and y-axes. If the prescribed POM Ξ(t) coincides with the spot's natural motion, then all prefactors vanish identically and f Gold disappears everywhere in Q. Importantly, the control signal is localized around the spot position and vanishes far from it because the spatial derivatives of its profile decay sufficiently fast, lim ξ →∞ ∇ ξ U c = 0. Alongside with these advantages, limitations in the applicability of Goldstone mode control (3.6) exist as well. For instance, the magnitude of the applied control may locally attain values that are unfeasible to realize physically because f Gold is proportional to the slope of the wave profile U c . Additionally, the stability of the control scheme depends sensitively on how precise the Goldstone modes can be calculated. Further, the complete spatial domain Ω accessible by the spot has to be available for the control as well. As already mentioned above, f Gold cannot be applied to desired trajectories U d which do not comply with initial as well as boundary conditions or which are non-smooth. In contrary, optimal control can deal with many of these complications. An arbitrary (even non-smooth) desired distribution U d is approximated as close as possible by a control f that has to be found by numerical optimization. However, optimal control is computationally much more expensive than Goldstone mode control.
If the desired U d is not realizable by some control f , then the method of this subsection is not applicable. In this case, an optimal control problem can be solved that determines the best approximation of Ud in the L 2 (Q)-norm. This is the issue of the next subsection.

Optimal control
An optimal control f minimizes a so-called objective functional J that is in our case of tracking-type U satisfies the controlled state equation associated to f with respect to given initial and boundary conditions, cf. equations (3.1). The first term appearing in J measures the distance between the actual and the desired solution U and U d up to the terminal time T in an L 2 (Q)-sense. In the second, so-called Tikhonov regularization term, a small positive number ν guarantees the existence of an optimal control f opt that minimizes the objective functional (3.8) for Ω ⊂ R q , q = 1, 2, 3, see reference [6]. If the desired state U d is not exactly realizable and ν = 0, then the functional (3.8) would not have a minimum. Then the optimal control problem is unsolvable. However, for any exactly realizable desired state U d and ν = 0, the solution f to (3.2) is the optimal control for (3.8). If U d is not exactly realizable, then an approximation of U d is obtained as part of the solution to the optimal control problem. The minimization of J must be performed with respect to state U and control f . Expressing U in terms of S(f ), where S : f → U is the solution operator to (3.1) in Q, justifies the definition of a reduced objective functional J(f ) := J(S(f ), f ). In order to minimize J by a descent method, its first directional derivative with respect to f in the direction h is needed that can be determined by the chain rule The state U is constrained to satisfy the controlled state equation together with given initial and boundary conditions, cf. equations (3.1). By a Lagrange multiplier P(r, t) = (p 1 (r, t), . . . , p n (r, t)) T , also called the adjoint state, the state equation can be "eliminated" to simplify (3.9), (3.10) The adjoint state P is the solution of the adjoint equation subject to terminal condition P(·, T ) = 0 in Ω and periodic boundary conditions in ∂Ω, where DR T denotes the transposed Jacobian matrix of R with respect to U. In the minimum, i.e., for an optimal control f = f opt with state U = U opt and adjoint state P = P opt , the derivative J (f opt )h has to be zero in all directions h. In view of (3.10), it is rather obvious that this amounts to Table 1. Parameter values used in the numerical simulations. The parameters θ = 1, κ 2 = 2, κ 3 = 1, and κ 4 = 8.5 are the same for set 1 and set 2.
This is nothing more than the well-known condition that, in a minimum, the gradient of the function to be minimized is zero.
Due to the mixed initial and terminal conditions for U and P, it is rarely possible to find numerical solutions to optimal control by a direct integration method. To reduce numerical costs, we employ Model Predictive Control and split our optimal control problem in subproblems with a 4 time-step small time-horizon [42]. Thereby, each subproblem is solved with a gradient or conjugate gradient optimization method that employs formula (3.10) for gradient computations. We refer to [4,7] for the optimal control of traveling wave fronts or spiral waves and to [8] for optimal control of localized spots. Details on the iteration scheme are discussed in the supplementary information (SI), paragraph S1.

Examples
In the following, we discuss three examples for position control of traveling spot solutions to the threecomponent RD model (2.1). Mainly, we compare Goldstone mode control f Gold with optimal control f opt . If not stated otherwise, the state equation (3.1a) and the adjoint equation

Translational position control of spots
In our first example, we aim to shift the spot's position along a Lissajous curve without controlling its orientation, i.e., the spot's symmetry axis is kept frozen to the x-axis. Thus, the POM Obviously, the control is localized at the current spot position Φ(t) and vanishes far away from it. Despite that the average speed v = L curve /T ≈ 6v x 0 along the studied Lissajous curve (4.1) with arc length L curve is almost five times larger than the propagation velocity of the uncontrolled spot, the magnitude of f u,Gold is of the same order as the 2) between f u,Gold (3.6), and optimal activator control signals f u,opt (3.8) during t ∈ [0, T /2]. We select set 1 in Table 1 for the parameters to (2.1) and set the Tikhonov parameter to ν = 10 −7 .
local reaction terms (2.1). The control signals applied to the inhibitors v and w are one and two magnitudes smaller [SI video1] than the activator's control, respectively.
Graphically, there is no distinguishable difference between f Gold and f opt , compare [SI video1]. Both are always localized close to the current spot position, and their magnitudes change proportional to |Φ(t)|. For a quantitative comparison, we compute the relative errors between f u,Gold and f u,opt measured by the L 1 (Ω)-norm Here, |h (r, t) | indicates the absolute value of h at position r and time t.
In Figure 3c, we depict solely the normalized error for the first half of the protocol because it starts to repeat after T /2, Φ y (t) = −Φ y (t + T /2). The relative error between f u,Gold and optimal control f u,opt (solid line) ranges between 2% and 8%. As reported in Appendix A.1, the limiting error is dominated by the time step chosen in the implicit Euler scheme. Albeit the scheme is A-stable, the error at a specific time t is of the order of O(dt). Consequently, we observe that f u,Gold (t) − f u,opt (t) L 1 (Ω) is bounded from above by dt; dt = 0.1 in the studied example.
The dashed line in Figure 3c shows the relative error between the activator distribution obtained by Goldstone mode control and the one calculated under optimal control, u Gold (t) − u opt (t) L 1 (Ω) / u Gold (t) L 1 (Ω) . At any time, this error is less than 10 −3 , i.e., both controlled states agree remarkably well, despite that f u,Gold (t) − f u,opt (t) L 1 (Ω) is of the order 10 −1 . Additionally, the relative errors between the desired distribution U d and the state solutions U Gold and U opt are less than 10 −7 in both cases (not shown explicitly). This confirms that the Goldstone mode control (3.6), within numerical accuracy, is indeed the solution to the corresponding unregularized optimal control problem. Similar conclusions had been obtained in our previous study of position control of front solutions in one spatial dimension, see [42].
The gradient-type method, used to solve the optimal control problem, relies on an initial guess for the control signal. The closer the starting guess is to the final solution, the fewer iteration steps are necessary to converge for most established optimization methods. Starting every iteration with an initial zero control, it takes on averagē n iter 23 iterations per time step for position control along the Lissajous curve (4.1). Using the control solution of the previously solved subproblem as initial guess reduces the average number of iterations ton iter 14. Taking advantage of the similarity between f Gold and f opt , see Figure 3c, the computational costs reduce even further. The most substantial computational speed-up is obtained by initiating every optimization subproblem with (3.2). Then, the iteration stops on average after the first step,n iter 1.

Stability of position control
Any open-loop controls is sensitive against perturbations of the initial conditions, data uncertainty, or numerical roundoff errors. To test the stability of our Goldstone mode control for position control f Gold , we accelerate or decelerate a single spot from its initial, intrinsic velocity v 0 to a final velocity v 1 using a translational POM for i ∈ {x, y}. Note that both the protocol's velocityΦ(t) and accelerationΦ(t) are continuous functions within the interval [0, T i ]. T i denotes the duration of the protocol. The maximum acceleration π v 1 i and inversely proportional to T i . A sketch of the protocol is depicted in Figure 4b. Since the proposed control scheme is an open-loop control, deviations between the current spot position Φ curr (t) and the POM Φ(t) will grow unbounded in time if the difference between them exceeds a critical value [28]. A specific protocol is called stable and marked by green boxes in Figure 4 if and only if the Euclidean distance is bounded as Φ curr (t) − Φ(t) < L/2 for all times t ∈ [0, t end ]. Otherwise, it is called unstable (red boxes). Note that a protocol is also considered to be unstable if the control leads to the nucleation of additional spots. In order to make the results comparable for different protocol durations, we adjust the terminal simulation time t end according to t end = max 10 t drift , T i + 10 L/|v i 1 | with drift time t drift = L/v x 0 . We stress that all simulation results presented in Figure 4 have been computed for sufficiently long time intervals and do not alter upon an increase of the total simulation time. Figure 4a depicts the numerically evaluated region of stable position control (green boxes) in x-direction as a function of the ratio of terminal spot velocity v x 1 to the initial one v x 0 and the ratio of the control duration T x to the drift time t drift . The translational POM in y-direction is set to zero, Φ y (t) = 0. As expected, the numerical algorithm is stable in the absence of control, v x 1 /v x 0 = 1. Further, it turns out that the control scheme is mostly stable for rapid, T x t drift , to moderately slow POMs, T x 10 t drift , regardless of the velocity change, |v x 1 − v x 0 |. The stability regions exhibit an asymmetry with respect to the sign of the velocity change. Weakly accelerating protocols, 1 < v x 1 /v x 0 2, are unstable (red colored region) while decelerating ones, v x 1 < v x 0 , are always stable  Table 1 and thus the drift time is given by for T x 10 t drift . This finding is in agreement with [28]. The instability for v x 0 < v x 1 2v x 0 is caused by an undesired rotation of the spot induced by numerical truncation errors. These accumulate during the simulation and eventually result in an asymmetric perturbation (with respect to y) acting on the spot pattern. Once the spot starts to rotate and eventually drifts away from the centerline y = 0, the proposed open-loop control f Gold can neither respond nor correct the undesired rotation. The impact of the numerical truncation error becomes more pronounced with growing protocol's duration T x and results in a broad unstable region for long protocols, T x > 10 t drift .
The situation changes if one aims to move the spot pattern perpendicular to its intrinsic direction of propagation, here in y-direction. In Figure 4c, we keep the motion in x unchanged, Φ x (t) = v x 0 t, and accelerate the spot according to (4.3) along the y-direction. Because the controlled spot solution is symmetric with respect to the centerline y = 0, position control in y might be inherently unstable [28]. One notices immediately that regions with unstable position control are much larger compared to Figure 4a. Nevertheless, the control is stable for weak acceleration, v y 1 0.1v x 0 , independent of the protocol's duration. Increasing the terminal velocity further, Goldstone mode control starts to fail. Once a certain deviation between the current spots' position and the proposed POM is attained, the pattern cannot follow the applied control anymore and starts to move freely. With further growing terminal velocity v y 1 , the control's magnitude increases as well and thus successful position control can be re-stabilized. Longer protocols T y result in an accumulating of numerical truncation errors.

Orientation control with speed adjustment
In the previous paragraph, we have demonstrated that the stability of position control can be enhanced if in any current position of the spot its symmetry axis, given by Φ ϕ (t), points tangentially to the direction of motion. Therefore, in our next example, we propose to shift the spot pattern along a circular trajectory by simultaneously controlling its orientation Here, r denotes the radius of the circle and T the protocol's duration. For experimental realization compare [41], for example. In Figure 5, we present the temporal evolution of the activator distribution u (a) controlled by f Gold (b). In line with the POM, the spot always keeps its symmetry axis at the tangent to the desired trajectory of motion. The control f Gold remains localized and is dominated by the translational Goldstone mode ∂ ξx U c due to the acceleration along Ξ(t) (4.4); the average speed isv 2.4v x 0 . Notably, the maximum value of the control magnitude is half as strong compared to position control without adjusting the orientation, Φ ϕ (t) = 0 (not explicitly shown). In panel (c), the temporal behavior of the relative error f u,Gold (t) − f u,opt (t) L 1 (Ω) / f u,Gold (t) L 1 (Ω) measured by the L 1 (Ω) norm is shown (solid line). They are large compared to pure position control along a Lissajous curve, cf. Figure 3. Stronger deviations are caused by interpolation errors arising during numerical rotation of spot patterns by Φ ϕ (t). The relative error attains a maximum at Φ ϕ (t) = m 45 • , m odd. At these angles, the distance between the nodes of the rotated grid and the underlying one is the largest, viz., dx/ √ 2, and, hence, numerical interpolation errors become significant. Contrarily, the relative error minimizes at Φ ϕ (t) = m 90 • , m ∈ Z, at which both grids coincide. Remarkably, the normalized error u Gold (t) − u opt (t) L 1 (Ω) / u Gold (t) L 1 (Ω) (dashed line), is still less than 10 −3 at any instants of time despite that the deviation of the associated controls rises up to ∼ 25%.

Orientation control
If the uncontrolled spot propagates at non-zero velocity v 0 = 0, the simplest way to navigate it through a spatial domain is to control exclusively its current orientation Φ ϕ (t). If so, the translational components of the POM Ξ(t) are determined byΦ .  We select set 1 in Table 1 for the kinetic parameters to (2.1) and set Tikhonov parameter to ν = 10 −7 . Now, we pick up the problem formulated in Section 2, namely, how to prevent pinning of a spot at a local heterogeneity in the domain. The heterogeneity is viewed as circular region where the parameter κ 1 jumps from a background value κ back 1 , ∀r / ∈ Ω • to a defect value κ het 1 , ∀r ∈ Ω • whereby Ω • = {(x, y) ∈ R 2 : (x + R) 2 + y 2 < R 2 } with radius R = 0.05. The orientational POM for avoiding the heterogeneity is set to with duration T = L x /v x 0 . Note that the corresponding prescribed positions (Φ x (t), Φ y (t)) T have to be calculated numerically.
In Figure 6, we present the temporal evolution of the activator distribution u in panel (a) and the corresponding control signal f u,Gold in panel (b). The prescribed translational POM is indicated by the dashed lines. At first glance, the control signal possesses a more complicated shape and its magnitude is significantly reduced, |f u,Gold | 10 −3 , as compared to |f u,Gold | 10 0 and |f u,Gold | 10 −1 in the previous examples, cf. Figure 3 and Figure 5. Thus, orientation control is less invasive than position control. In return, we lose the ability for  Table 1. The circular defect with radius R = 0.05 is modeled by a jump in κ 1 from its background value of κ back 1 = −7.30 to the value inside the heterogeneity κ het 1 = −7.50.
fast intervention into spot dynamics as well as for mayor increase in the speed of the spot. Additionally, orientation control is much more susceptible to fail. The small control magnitudes are too weak to suppress the impact of numerical round-off errors which may result in undesired spot rotation, cf. Section 4.1.1. Caused by the small propagation velocity, the duration T of the POM grows as compared to position control, see Section 4.1, and therefore the probability of failure increases as well.

Position control by a single control signal
So far, we have discussed examples of fully actuated systems for which the number of state components equals the number of independent control signals. If the coupling matrix B is not invertible, expression (3.6) for f Gold cannot be used. The question arises how to extend our approach to underactuated systems [26,31]. In the following example we assume a control acting on the activator u only while inhibitors v and w remain uncontrolled, i.e., f v,Gold (r, t) = f w,Gold (r, t) = 0. Control via an inhibitor has been discussed in detail for the Hodgkin-Huxley model and the three-component Oregonator model for photosensitive BZ reaction, compare supplemental information to [26].
To derive an expression forf u,Gold (r, t), we start with the fully actuated system  Table 1 for the calculations.
Equations (4.8b)-(4.8c) are linear, inhomogeneous PDEs with initial conditions v(r, t 0 ) = v 0 (r) and w(r, t 0 ) = w 0 (r), respectively. Their solutions can be written as where K 0 i and K i , i ∈ {v, w}, are integral operators involving Green's functions to the homogeneous equations corresponding to (4.8b)-(4.8c) with associated initial conditions and to the inhomogeneous equations with zero initial conditions. Plugging (4.9) into (4.8a) gives From the last line of (4.10) we identify the expression for f u,Gold (r, t) to bẽ whereby the component of f Gold are determined by (3.6).
As an example for position control by a single control signal, we guide a spot along the Lissajous curve given by (4.1) with radius r = 0.2 and protocol duration T = 200. The spot's orientation Φ ϕ (t) = 0 remains uncontrolled. The relative errors between desired and controlled states are shown in Figure 7 (bottom). All states are obtained from numerical simulation of (2.1)-(3.1a) with control f Gold = (f u,Gold , 0, 0) T given by (4.11). The relative error for the activator u (solid line) is less than 10 −3 at any time t and thus the controlled activator pattern agrees satisfactorily well with the desired distribution. This finding is corroborated by snapshots of u at different instants of time in Figure 7 (top left). In contrast to the activator, the profile of the inhibitor v is not preserved under control but deformed considerably, see Figure 7 (top right). In particular, an elongated region of activity becomes apparent along the Lissajous curve due to time scale separation in the RDS (2.1). The concentration of the slow inhibitor v, produced in the wake of the activator, decays exponentially to the rest state on a time scale τ = 48 ≈ T /4. Consequently, the relative error of v (dashed line) attains relatively large values of the order 10 −2 . On the other hand, the fast inhibitor w and the activator u vary on the same characteristic time scale as θ = 1 was chosen in the considered example. Thus, we expect only small changes in both profiles in the presence of the control. In fact, the values of the relative error for w turn out to be less than 10 −4 which is even one magnitude smaller than the relative error of u, cf. the dash-dotted line in Figure 7 (bottom).

Conclusion
Localized traveling chemical, chemo-mechanical, electrical or neural activity is ubiquitous in spatially extended nonlinear systems driven far from thermodynamic equilibrium. Hence, to control the current position, orientation and velocity of a traveling spot is a key challenge not only under general aspects but particularly from the perspective of numerous applications. We have demonstrated that the control signals, which one has to apply to solve these tasks, can be analytically expressed by the Goldstone modes of the uncontrolled spot. The control acts locally as long as the Goldstone modes are localized, and preserves the spatial profile of the spot as long as the spectral gap between deformation modes and Goldstone modes in the linear stability operator of the uncontrolled spot solution is sufficiently large.
We would like to stress that a complete knowledge about the underlying nonlinear kinetics R(U) responsible for the spontaneous formation of the spot is not required in order to determine the control f Gold for a given protocol of motion. According to (3.6), we can express the control amplitude through the derivative of the uncontrolled spot profile U c , i.e., R(U) does not enter explicitly into the expression for f Gold . Thus, if R(U) is known, as in case of the RD model (2.1), we numerically solve the nonlinear eigenvalue problem (3.3) and calculate f Gold from (3.6). However, often in practical applications the complete scheme of kinetic steps specifying R(U) is known only approximately because some kinetic steps are unclear or hidden. Nevertheless, if in those cases the stationary profile and the propagation velocity of the uncontrolled spot can be measured with sufficient accuracy, then the method works. In our view, this aspect might extend the applicability of Goldstone mode control considerably.
Strictly speaking, localized spots are solutions on an infinite spatial domain. Our results are correct as long as the spot weakly interacts with the domain boundaries. This will be the case if the domain size is much larger than the characteristic diameter of the spot, l. When we compare the analytical results with numerical simulations, we use periodic boundary conditions with spatial periods much larger than l. The interaction range of a spot with a Neumann boundary is determined by the so-called response functions. These are the eigenfunctions of the adjoint linear stability operator of the uncontrolled spot. Usually, they are well localized and, therefore, the effective interaction between spots or between a spot and a Neumann boundary is short-ranged.
Goldstone mode control is realized by external spatio-temporal forcing, i.e., it is an open-loop control. Contrary to closed-loop or feedback control, continuous monitoring of the system is not required. On the downside, as any open-loop control, the method is sensitive to perturbations. Therefore, the range of applicability has been checked by a stability analysis.
Based on general symmetry considerations, Goldstone mode control is widely applicable. So far, the method has been successfully used to guide traveling interfaces and excitation pulses in 1D [26,42] and spiral waves [42] as well as to shape iso-concentration lines of traveling wave patterns [29] in 2D. Recently, we successfully applied Goldstone mode control to spot solutions of neural field equations [56] that phenomenologically describe the dynamics of synaptically coupled neurons [11]. Remarkably, in all examples considered so far, Goldstone mode control is, within numerical accuracy, equal to solutions of an equivalent, non-regularized optimal control problem. Consequently, our control turns out to be optimal, i.e., no other control enforces the system closer to the desired target state according to the protocol of motion. Furthermore, these control signals have been proven to be excellent initial conditions for regularized optimal control problems achieving a substantial computational speed-up. Generally, Goldstone mode control might serve as consistency check for numerical optimal control algorithms as well.
We emphasize that optimal control is not only computationally demanding but requires full knowledge of the nonlinear kinetics. On the other hand, the scope of optimal control can be extended in various ways: More general objective functionals than stated above can be studied, cf. [21]. Optimal control can also easily deal with control signals confined to prescribed spatial regions or with upper and lower bounds for the control amplitudes in the form of inequalities, i.e., −∞ < f a ≤ f i (r, t) ≤ f b < ∞, ∀i = 1, . . . , n. Moreover, by adding a multiple of the L 1 -norm of the control to the objective functional J, sparse optimal controls are obtained. Sparsity means that the optimal control vanishes in certain regions. The level of sparsity depends on a so-called sparsity parameter. The larger this parameter is the larger is the region where the optimal control is exactly zero. In this way, the optimal control is concentrated on regions that are most important for the minimization. However, this is a more intuitive interpretation of sparsity that cannot be exactly quantified. For technical details as well as examples, we refer to [7-9, 21, 42].

Appendix A. Details on numerical methods
In our numerical simulations, the state equation is solved on a rectangular domain Ω = (x a , x b ] × (y a , y b ] with periodic boundary conditions Here, U(r, t) = (u 1 (r, t), . . . , u n (r, t)) T is the vector of n ∈ N state components defined in the two-dimensional spatial domain Ω ⊂ R 2 with r = (x, y) T . Further, R(U) = (R 1 (U), . . . , R n (U)) T describes the nonlinear kinetics of the components and the n × n matrix D of constant diffusion coefficients is assumed to be diagonal D = diag(D 1 , . . . , D n ) (isotropic medium). On the right hand side of (A.1a), the space-time dependent control signals are denoted by f (r, t) = (f 1 (r, t), . . . , f m (r, t)) T , m ∈ N. Without loss of generality, we fix the spot's direction of motion to coincide with the x-axis, i.e., v x 0 = 0 and v y 0 = 0. Any numerical simulation to equation (A.1a) are initialized with the profile U c . The latter and the natural velocity v 0 are obtained by solving the nonlinear eigenvalue problem with adequate accuracy.

A.1 Simulations based on the Goldstone mode control f Gold
In all position control examples based on Goldstone mode control, we numerically solve (A.1) with spectral methods [12]. Transforming (A.1a) into the Fourier space, one gets for the n-th state component whereû n denotes the Fourier transform F{·} of the n-th component,û n (k, t) = F{u n (r, t)}, and k = (k x , k y ) T is the wave vector. Multiplying (A.3) by e −cnt with c n = D n k 2 and integrating the equation over a single time step from t m to t m+1 = t m + dt, one derives the exponential time differencing (ETD) method This formula is exact. The essence of ETD methods is in approximating the integral in this expression. We use ETD2 implying that Finally, one arrives at the numerical scheme for ETD2 which has a local truncation error of O(5 dt 3F {·}/12) [12]. To avoid interpolation errors due to discretization of the spatial domain, both the desired solution , are computed with spectral methods 1 . First, the uncontrolled spot profile U c (r) is rotated around (0, 0) T in the spatial domain by Φ ϕ (t), see Figure A.1(b). In a second step, the shift as given by the translational POM  Figure A.1(c), is computed in Fourier space by multiplying the Fourier transform of the previously rotated spot solution with exp (−ik · Φ(t)), followed by the inverse Fourier transform F −1 {·}.
Analogously, we calculate f Gold : First, the derivatives with respect to x and y are obtained numerically in Fourier space using ∂ x U c = F −1 {ik x F{U c }} and ∂ y U c = F −1 {ik y F{U c }}, respectively. The angular derivative is given as the linear combination ∂ ϕ U c = −y∂ x U c + x∂ y U c . Second, all derivatives are rotated according to the orientational POM Φ ϕ (t), followed by a shift computed in Fourier space.
In all simulations using ETD2, the squared domain Ω = (−0.5, 0.5] × (−0.5, 0.5] is discretized into 256 × 256 spatial grid points and Fourier modes, respectively. We emphasize that the domain size is chosen sufficiently large to avoid self-interaction of the spots in the periodic simulation domain. The numerical time step dt is adjusted such that the local truncation error is less than 10 −4 for any protocol of motion. But dt is always less or equal to 0.01.

A.2 Numerical calculation of optimal control
For optimal control, both the state equation (A.1) and the adjoint equation have to be solved for any iterate f . In (A.7) U is the state associated to f . It has to be computed prior to solving (A.7). The adjoint state P has to obey the terminal condition P(·, T ) = 0 in Ω at final time T and periodic boundary conditions in ∂Ω. In equation (A.7), DR T denotes the transposed Jacobian matrix of R with respect to U. The solution to the optimal control problem is determined by solving the system for f . Due to the mixed initial and terminal conditions for U and P it is rarely possible to find numerical solutions to optimal control by a direct integration method. To reduce numerical costs, we employ Model Predictive Control and divide our optimal control problem into subproblems with a 4 time-step small time-horizon [42]. Thereby, each subproblem is solved with a gradient-type method which proceeds as follows: Starting for k = 1 with an initial guess f k−1 ≡ f 0 for the control, we compute U k as the solution to (A.1) as well as P k as the solution to (A.7) with U k substituted for U. The gradient of the non-negative tracking-type functional J with respect to f , d k , is calculated as d k = P k + ν f k−1 whereby −d k defines the direction of steepest descent in function space. A new control f k is iteratively obtained by where s denotes a suitable step size [20,21]. If d k satisfies an appropriately chosen termination condition, in our work this is set to d k < 10 −8 , then the algorithm stops and f k is taken as f opt . Otherwise, the steps are repeated iteratively for k := k + 1. Because the convergence of a gradient-type method can be fairly slow, the maximum number of iterations has been limited to N iter = 50. Moreover, we keep the Tikhonov-regularization parameter fixed at ν = 10 −7 . Many modifications of the basic algorithm sketched above exist and often lead to a better performance, as, e.g., the nonlinear conjugate gradient method [4]. However, all iterative algorithms require multiple solutions of state and adjoint state equation, resulting in a rapidly increasing computational cost of numerical optimal control with the number of spatial dimensions and the length of the time interval. Therefore, the domain Ω is discretized by 256 grid points in each direction with spatial step sizes of dx = dy = 1/256 and the Laplacian ∆ is approximated by a 5-point stencil. Due to significant computational cost of the gradient-type algorithm, we choose a moderate temporal resolution of dt = 0.1. For maximum stability the time evolution is computed with an implicit Euler method.

A.3 Error estimation
While the numerical evaluation of f Gold is limited by the accuracy of the first spatial derivatives∇ ξ U c , the numerical computations of f opt is affected by errors arising both in the discretization of space and time.
We stress that f Gold is calculated using spectral differentiation, while we have to use finite difference stencils for the simulation of both the state (A.1) and adjoint equation (A.7) in order to reduce computational costs. Nevertheless, the limiting error arises from the time step chosen in the implicit Euler-scheme. Despite that the latter is A-stable, the local truncation error is O(dt 2 ) and the error at a specific time t is of the order of O(dt). Consequently, one can expect that the relative error f u,Gold − f u,opt is in the L 2 -norm bounded from above by O(dt).

A.4 Simulation times
As already mentioned, optimal control algorithm requires much more computational resources compared to the analytic expression f Gold . For instance, the computation of optimal position control along the Lissajous curve, see Figure 4, with a short period of T = 200, i.e., 2000 time iterations steps with dt = 0.1, takes roughly 50 hours by using 3 cores (Intel i5 − 6500 with 3.20GHz) in parallel. Contrary, the same numerical simulation of (A.1) based on f Gold runs approximately one hour on the same PC. The relatively short simulation times enable us not only to use integration schemes with smaller numerical truncation error, e.g., spectral methods with ETD2 [12], but also to move the spot along longer, more complex trajectories [SI video2,SI video3].

Appendix B. Traveling spot interacting with a circular heterogeneity of jump type
The interaction of traveling spots with different types of parameter heterogeneities in 1D and 2D has been studied in detail by many authors, see [36] and references therein. Penetration, rebound, annihilation, oscillation, as well as stationary or oscillatory pinning of spots were observed. Figure   well as spots orbiting both inside and outside of the defect boundary. If the activator describes the temperature in a catalytic packed-bed reactor, resting hot spots [32] or those pinned to local heterogeneities can damage the catalyst support. In particular, collision of hot spots with the reactor walls must be prevented for safety reasons. Consequently, guidance of a traveling spot with given velocity along a desired trajectory through a bounded spatial domain might be a particular challenge in chemical engineering applications.
In the numerical simulations we realize a local heterogeneity in the medium by a circular defect of jump type (step-function) in the additive parameter κ 1 : where the circular domain with radius R is centered at (0, 0) and is defined by To start a numerical collision experiment, we initiate one spot moving rightwards (v x 0 > 0) at the far left edge of the simulation domain and measure the response to a circular defect with radius R = 0.1 located at (0, 0).
The outcomes of collisions are classified in Figure  is smaller than the background value, the single initiated spot splits into two when reaching the left edge of the heterogeneity. Thereby, the initial velocity vector is conserved. Increasing the value of κ het 1 and thus decreasing the difference , the spot is able to pass the defect. With growing difference , spots first enter the defect and split into multiple solutions when leaving the obstacle at the right edge. Further increasing results in multiple creation of spots inside the defect. For κ het 1 κ back 1 , the spot is able to enter the defect, gets trapped inside, and initiates the creation of multiple spots.
The observed scenarios of spot-defect-interactions are different for the second parameter set in table 2 and κ back 1 = −7.30. There, one discovers three qualitatively different regimes: trapping, transmission, and reflection, Figure B.1(b). For κ het 1 κ back 1 , the spot pins at the outer edge of the defect and starts to move around the heterogeneity. With shrinking difference , the spot gets first reflected for κ het 1 < κ back 1 , then is able to pass the defect for ≈ 0, and eventually gets trapped inside the heterogeneity for ∈ [0.09, 0.26]. In the last scenario, the trapped spot eventually performs circular motion at the inside of the defect. A further increase of κ het 1 results in various reflection scenarios.
Finally, we remark that parameter gradients and the geometrical shape of a heterogeneity will impact the interaction between traveling spots.

Appendix C. Asymptotic perturbation analysis -projection method
We presume that the uncontrolled reaction-diffusion system (A.1a) possesses a stable traveling wave (TW) solution U c which either propagates smoothly with velocity v 0 = (v x 0 , v y 0 ) T or rotates around its center of rotation with angular velocity ω 0 . Such a pattern satisfies in the frame of reference ξ ≡ (ξ x , ξ y ) T = A(−ω 0 t)(r − v 0 t) co-moving with the velocity v 0 and co-rotating with angular frequency ω 0 . Here, A(α) = [cos(α) , − sin(α); sin(α) , cos(α)] is the clockwise rotation matrix in 2D. In the co-moving coordinates ξ, the angular derivative is given by ∂φ = −ξ y ∂ ξx + ξ x ∂ ξy . Obviously, (C.1) depends explicitly on time t if the TW pattern simultaneously moves and rotates, ω 0 v 0 = 0. Consequently, we have to claim that the considered pattern either moves, solely rotates, or is at rest, i.e., |v 0 | = ω 0 = 0. In what follows, we keep both characteristic quantities in our derivations but always have in mind that at least one of them must be zero.

C.1 Goldstone modes and response functions
The stability of the TW pattern U c is determined by the eigenvalues of the linear stability operator L. The latter arises by expanding (C.1) around U c in the frame of reference ξ. Thereby, DR (U c ) denotes the Jacobian matrix of R evaluated at U c . Since we presume that U c (ξ) is stable, the eigenvalue of L with the largest real part is λ 0 = 0 and the corresponding eigenfunctions W i (ξ), i ∈ {ξ x , ξ y ,φ}, also called the propagator modes, can be expressed by the Goldstone modes ∂ i U c (ξ), i ∈ {ξ x , ξ y ,φ} . Calculating the derivative of (C.1) with respect to ξ x , ξ y , andφ, one gets These relations can be re-written in a compact form The adjoint operator L † to L is defined with respect to the standard inner product in two-dimensional function space where u denotes complex conjugation; yielding Next, we introduce the so-called response functions W † i (ξ), i ∈ {ξ x , ξ y ,φ}. Although the latter do not coincide with the eigenfunctions of L † to eigenvalue zero, they are closely related to these. The response functions are also not identical to W i , i ∈ {ξ x , ξ y ,φ} because L is in general not self-adjoint. Regarding the Goldstone modes and response functions, one can distinguish three cases corresponding to (i) stationary solutions, (ii) traveling localized patterns like spots, and (iii) rotating solutions as e.g. spiral waves.
(i) For stationary, immobile solutions with v x 0 = v y 0 = ω 0 = 0, all Goldstone-modes ∂ i U c , i ∈ {ξ x , ξ y ,φ}, are ordinary eigenfunctions of L to eigenvalue λ 0 = 0. Additionally, the response functions are also ordinary eigenfunctions of L † to λ 0 = 0. For stationary patterns, L can be expressed by the product of a diagonal square matrix M and a self-adjoint matrix L, L = M L. Then, the eigenfunction W † i are determined by the relation W † i = M −1 ∂ i U c , i ∈ {ξ x , ξ y ,φ}, for a given Goldstone mode ∂ i U c [18]. (ii) For traveling localized patterns with ω 0 = 0, both translational Goldstone modes ∂ ξx U c and ∂ ξy U c are ordinary eigenfunctions of L to λ 0 = 0. Similar, the response functions W † ξx and W † ξy are determined by L † W † i = 0, i ∈ {ξ x , ξ y }. On the other hand, the rotational Goldstone mode ∂φU c is a generalized eigenfunction of rank k = 2 to eigenvalue 0 such that L 2 ∂φU c = 0 or L∂φU c = v y 0 ∂ ξx U c − v x 0 ∂ ξy U c . In a similar fashion, the corresponding rotational response function can be obtained as a solution of In the case of solely rotating solutions with v 0 = 0, only the rotational Goldstone mode ∂φU c (ξ) is an usual eigenfunction of L to eigenvalue zero, L∂φU c (ξ) = 0. The same holds for the adjoint operator and the rotational response function, i.e., L † W † ϕ (ξ) = 0. Although the translational Goldstone modes ∂ i U c , {ξ x , ξ y } are not eigenfunctions of L, linear combinations of them, viz. Y ∓ (ξ) = ±i∂ ξx U c + ∂ ξy U c , yield eigenfunctions to eigenvalues λ ±1 = ±iω 0 , LY ∓ (ξ) = ±iω 0 Y ∓ (ξ). Since eigenvalues of L and the adjoint operator L † are complex conjugate to each other, there must be corresponding eigenfunctions Y † ∓ (ξ) to the adjoint operator L † , L † Y † ∓ (ξ) = ±iω 0 Y † ∓ (ξ).

C.3 Inverse problem -solving for the control
Folllowing [26], now we consider the inverse problem and view (C.17) as a conditional equation for the control f . The latter is assumed to be a linear superposition of Goldstone modes B f (r, t) = K 1 ∂ x U c (ξ(t)) + K 2 ∂ y U c (ξ(t)) + K 3 ∂ ϕ U c (ξ(t)) , (C. 19) where U c is evaluated at the presumed position in spaceξ(t) = A (−Φ ϕ (t)) (r − Φ (t)). Plugging (C. 19) into 20) and, hence, the expansion coefficients K j are given as Finally, we obtain the analytical expressions for Goldstone mode control f Gold based on asymptotic perturbation analysis .