Controlling Canard Cycles

Canard cycles are periodic orbits that appear as special solutions of fast-slow systems (or singularly perturbed Ordinary Differential Equations). It is well known that canard cycles are difficult to detect, hard to reproduce numerically, and that they are sensible to exponentially small changes in parameters. In this paper we combine techniques from geometric singular perturbation theory, the blow-up method, and control theory, to design controllers that stabilize canard cycles of planar fast-slow systems with a folded critical manifold. As an application, we propose a controller that produces stable mixed-mode oscillations in the van der Pol oscillator.


Introduction
Fast-slow systems (also known as singularly perturbed ordinary differential equations, see more details in Section 2) are often used to model phenomena occurring in two or more time scales. Examples of these are vast and range from oscillatory patters in biochemistry and neuroscience [18,26,6,25], all the way to stability analysis and control of power networks [10,14], among many others [41,Chapter 20]. The overall idea behind the analysis of fast-slow systems is to separate the behavior that occurs at each time scale, understand such behavior, and then try to elucidate the corresponding dynamics of the full system. Many approaches have been developed, such as asymptotic methods [17,34,51,52], numeric and computational tools [24,32], and geometric techniques [20,31,33], see also [41,46,56]. In this article we take a geometric approach.
Although the time scale separation approach has been very fruitful, there are some cases in which it does not suffice to completely describe the dynamics of a fast-slow system, see the details in Section 2. The reason is that, for some systems, the fast and the slow dynamics are interrelated in such a way that some complex behavior is only discovered when they are not fully separated. An example of the aforementioned situation are the so-called canards [7,8,15], see Section 2.1 for the appropriate definition. Canards are orbits that, counter-intuitively, stay close for a considerable amount of time to a repelling set of equilibrium points of the fast dynamics. Canards are extremely important not only in the theory of fast-slow systems, but also in applied sciences, and especially in neuroscience, as they have allowed, for example, the detailed description of the very fast onset of large amplitude oscillations due to small changes of a parameter in neuronal models [18,26] and of other complex oscillatory patterns [9,12,47]. Due to their very nature, canard orbits are not robust, meaning that small perturbations may drastically change the shape of the orbit.
On the other hand the application of singular perturbation techniques in control theory is far-reaching. Perhaps, as already introduced above, one of the biggest appeals of the theory of fast-slow systems is the time scale separation, which allows the reduction of large systems into lower dimensional ones for which the control design is simpler [29,38,37]. Applications range from the control of robots [53,54,28], all the way to industrial biochemical processes, and large power networks [13,35,36,43,50,49]. However, as already mentioned, not all fast-slow systems can be analyzed by the convenient time scale separation strategy, and although some efforts from very diverse perspectives have been made [2,3,4,5,22,23,29,30], a general theory that includes not only the regulation problem but also the path following and trajectory planning problems is, to date, lacking.
The main goal of this article is to merge techniques of fast-slow dynamical systems with control theory methods to develop controllers that stabilize canard orbits. The idea of controlling canards has already been explored in [16], where an integral feedback controller is designed for the FitzHugh-Nagumo model to steer it towards the so-called "canard regime". In contrast, here we take a more general and geometric approach by considering the folded canard normal form, see Section 2.1. Moreover, we integrate control techniques with Geometric Singular Perturbation Theory (GSPT) and propose a controller design methodology in the blow-up space. Later we apply such geometric insight to the van der Pol oscillator where we provide a controller that produces any oscillatory pattern allowed by the geometric properties of the model, see Section 4.
The rest of this document is arranged as follows: in Section 2 we present definitions and preliminaries of the geometric theory of fast-slow systems and of folded canards, which are necessary for the main analysis. In Section 3 we develop a controller that stabilizes folded canard orbits, where the main strategy is to combine the blow-up method with state-feedback control techniques to achieve the goal. Afterwards in Section 4, as an extension to our previously developed controller, we develop a controller that stabilizes several canard cycles and is able to produce robust complex oscillatory patters in the van der Pol oscillator. We finish in Section 5 with some concluding remarks and an outlook.

Preliminaries
A fast-slow system is a singularly perturbed ordinary differential equation (ODE) of the form where x ∈ R m is the fast variable, y ∈ R n the slow variable, 0 < ε 1 is a small parameter accounting for the time scale separation between the aforementioned variables, λ ∈ R p denotes other parameters, and f and g are assumed sufficiently smooth. In this document the over-dot is used to denote the derivative with respect to the slow time τ . It is well-known that, for ε > 0, an equivalent way of writing (1) is where now the prime denotes the derivative with respect to the fast time t := τ /ε. One of the mathematical theories that is concerned with the analysis of (1)-(2) is Geometric Singular Perturbation Theory (GSPT) [41]. The overall idea of GSPT is to study the limit equations that result from setting ε = 0 in (1)- (2). Then, one looks for invariant objects that can be shown to persist up to small perturbations. Such invariant objects give a qualitative description of the behavior of (1)- (2). Accordingly, setting ε = 0 in (1)-(2) one gets known, respectively, as the reduced slow subsystem (which is a Constrained Differential Equation [55] or a Differential Algebraic Equation [42]) and the layer equation. The aforementioned limit systems are not equivalent any more, but they are related by the following important geometric object.
Definition 1 (The critical manifold). The critical manifold is defined as We note that the critical manifold is the phase-space of the reduced slow subsystem and the set of equilibrium points of the layer equation. The properties of the critical manifold are essential to GSPT, in particular the following.
Definition 2 (Normal hyperbolicity). Let p ∈ C 0 . We say that p is hyperbolic if the matrix D x f (p, 0, λ)| C0 has all its eigenvalues away from the imaginary axis. If every point p ∈ C 0 is hyperbolic, we say that C 0 is normally hyperbolic. On the contrary, if for some p ∈ C 0 the matrix D x f (p, 0, λ)| C0 has at least one of its eigenvalues on the imaginary axis, then we say that p is a non-hyperbolic point.
It is known from Fenichel's theory [19,20] that a compact and normally hyperbolic critical manifold S 0 ⊆ C 0 of (3) persists as a locally invariant slow manifold S ε under sufficiently small perturbations. In other words, Fenichel's theory guarantees that in a neighborhood of a normally hyperbolic critical manifold the dynamics of (1)-(2) are well approximated by the limit systems (3). Remark 1. Along this paper we use the notation S a 0 and S r 0 to denote, depending on the eigenvalues of D x f (x, y, 0, λ)| S0 , the attracting an repelling parts of the (compact) critical manifold S 0 . Accordingly, the corresponding slow manifolds are denoted as S a ε and S r ε . On the other hand, critical manifolds may lose normal hyperbolicity, for example, due to singularities of the layer equation, see Figure 1. It is in fact due to loss of normal hyperbolicity that, as in this paper, some interesting and complicated dynamics may arise in seemingly simple fastslow systems. Fenichel's theory, however, does not hold in the vicinity of non-hyperbolic points. In some cases, depending on the nature of the non-hyperbolicity, the blow-up method [27] is a suitable technique to analyze the complicated dynamics that arise. In the forthcoming section we introduce the particular type of orbits that we are concerned with and that arise due to loss of normal-hyperbolicity of the critical manifold; the so-called canards.

Planar Folded Canards
In this section we briefly describe folded canards and folded canard cycles in the plane. As we mention below, the adjective "folded" is due to a fold singularity. However, we remark that canards (and canard cycles) can be related to other types of singularities. The interested reader is refereed to, e.g. [15,39,57], references therein and, in particular, [41,Chapter 8] and [27,Section 3] for more detailed information.
The critical manifold is locally (near the origin) a perturbed parabola and is given by The (slow and fast) reduced flow corresponding to (7) is as shown in Figure 1.  (7) near the origin. The gray parabola depicts the critical manifold S 0 which is partitioned in its attracting S a 0 = S 0 | {x<0} and repelling S r 0 = S 0 | {x>0} parts, while the origin (the fold point) is non-hyperbolic. If α = 0 the origin is also called canard point. In this latter case, the orbit along the critical manifold is also known as singular maximal canard.

Remark 2.
To fix ideas, consider for a moment (7) with zero higher order terms 1 , that is 1 Refer to [39] for the much more complicated case that includes the higher order terms.
Then, it is straightforward to check that, for ε > 0 and α = 0, the orbits of (10) are given by level sets of Some orbits of (10) are shown in Figure 2, and in fact it is known [39] that canard cycles exist for H ∈ (0, 1 4 ). Compare with α = 0 in Figure 1.
What is remarkable is that there are orbits that closely follow the unstable branch of the critical manifold for slow time of order O(1). Such type of orbits are known as canards. There is a particular canard, which is called maximal canard and is given by {H = 0} that connects the attracting slow manifold S a ε with the repelling one S r ε . More relevant to this paper are periodic orbits with canard portions, which called canard cycles.
In the following section we design feedback controllers for (5) that render a particular canard cycle asymptotically stable. In other words, we consider the path following control problem where a canard orbit is the reference.

Controlling Folded Canards
We propose to study two control problems, namely which we call the fast control problem and to be referred to as the slow control problem. Recall that f and g stand for the higher order terms as in (8). The objective is to stabilize a certain reference canard cycle to be denoted by γ h .

Remark 3.
• The choice of the above control problems is motivated by applications, especially in neuron models, see [16,26,18], where the input current appears in the fast (voltage) variable and regulates the distinct firing patterns. However, if one is interested in the fully actuated case, a combination of the techniques presented here shall also be useful.
• Throughout this document we assume that one has full knowledge of the functions f and g. This means that for the fast (resp. slow) control problem we assume f = 0 (resp. g = 0). Otherwise one considers a controller of the form Notice that in the case of the fast-control problem (12), the controller changes the fast dynamics. This means that the controller can change the type of singularities the critical manifold may present.
To be more precise, consider for a moment (12) with u = −kx, k > 0, a simple proportional feedback controller. The closed-loop system then reads as for which the origin is now normally hyperbolic. This means that the feedback controller has changed the type of singularity (at the origin) from a fold to a regular one. It is clear that these type of controllers are not compatible with our task. So, we shall design controllers that do not change the type of singularity of the open-loop system. To formalize what we mean by "not changing the type of singularity", let us first recall the following definition: Definition 3 (k-jet equivalence). Let F : R n → R n and G : R n → R n be smooth maps. We say that F and G are (k-jet) equivalent at p ∈ R n if F (p) = G(p) and An equivalence class defined by the previous notion of equivalence is called the k-jet of F at p, and shall be denoted by j k F (p) [1].
Next we have a formal definition of what we refer to as a compatible controller : Definition 4 (Compatible controller). Consider a control systeṁ where ζ ∈ R n is the state variable, λ ∈ R p denotes system parameters (possibly including 0 < ε 1) and u ∈ R m stands for the controller. Suppose that for the open-loop system, that is when u = 0, the origin ζ = 0 ∈ R n is a nilpotent equilibrium point ofζ = f (ζ, 0, 0) and that there is a k ∈ N such that k is the smallest number so that j k f (0) = 0. Let u be a state-feedback controller, that is u = u(ζ, λ, ), where ∈ R m stands for parameters of the controller such as controller gains, and denote byζ = F (ζ, λ, ) the closed-loop system. We say that u is a compatible controller if the open-loop vector field f (ζ, λ, 0) and the closed-loop vector field F (ζ, λ, ) are k-jet equivalent at the origin for λ = 0.
We emphasize that once one fixes coordinates on R n , a k-jet equivalence between two maps means that such maps coincide on their partial derivatives up to order k.
As an example of the above definition, recall that a planar fast-slow system with a generic fold at the origin is given by with the defining conditions f (0, 0, 0) = 0, ∂f ∂x (0, 0, 0) = 0, ∂ 2 f ∂x 2 (0, 0, 0) = 0, and the non-degeneracy condition ∂f ∂y (0, 0, 0) = 0. Next, let u = u(x, y, ε) be a state-feedback controller and suppose one considers the fast-slow control system Then, u is a compatible controller if the closed-loop system verifies: ∂x 2 (0, 0, 0) = 0, and ∂F ∂y (0, 0, 0) = 0, which implies that the controller does not change the class of the singularity, since the origin is still a fold point of the closed-loop system.

Remark 4.
• The choice of the controller gain c 2 in Theorem 1 has an important impact in numerical simulations due to the fact of it appearing as an argument of the exponential function. The choice c 2 = 2 yields the better numerical results when stabilizing canard cycles, that is for h ∈ (0, 1 4 ). However, to stabilize the maximal canard (h = 0), it is necessary to choose c 2 < 2 to ensure that the controller remains bounded as y → ∞. See more detail in section 3.1.2.
• We recall that although from Theorem 1 one is able to stabilize any canard (because h ≤ 1 4 ), canard cycles exist only for h ∈ (0, 1 4 ), see Figure 2 and [39]. • The second item of Theorem 1 holds for any ε > 0.
The proof of Theorem 1 follows from the forthcoming analysis and is summarized in section 3.1.3. We show in Figure 3 a simulation of the results contained in Theorem 1.
(a) (b) (c) Figure 3: In all three columns we show, in the first row the (x, y) phase portrait of the closed-loop system (19) and in the second row the time-series of the corresponding controller. In all these simulations ε = 0.01. (a) The case for whichĝ = 0 and with parameters (α, c 1 , c 2 , h) = (−0.1, 1, 2, 1 4 e −400 ). We remark here that in order for the constant h = 1 4 e −400 to be numerically feasible one has to input h exp(c 2 yε −1 ) = 1 4 exp(c 2 yε −1 − 400) into the numerical integration algorithm. The desired canard cycle to be followed is shown in dashed-gray. (b) The maximal canard case withĝ = 0 and with parameters (α, c 1 , c 2 , h) = (0, 1, 2 − e −15 , 0). Notice that, indeed, trajectories follow the unstable branch S r ε for a large "height" and that the corresponding controller remains bounded. (c) An example of the effect of the extra term in (21) where we show two trajectories with the same initial conditions. The unstable one is obtained with the controller (20) while the stable one with (21). The desired canard cycle to be followed is shown in dashed-gray. The large spike in the controller is observed every time the trajectory crosses the y-axis long a fast fiber. For such simulation we have used (α, c 1 , c 2 , h) = (0, 5, 2, 1 4 e −400 ) and g = 100x(y − x 2 ). For more details see Sections 3.1.1 and 3.1.2.
As already anticipated, the idea is to design the controllerû in the blow-up space. Therefore, let us consider a coordinate transformation defined bŷ x =rx, y =r 2ȳ , ε =r 2ε ,û =r 2μ , α =rᾱ, (22) where (x,ȳ,ε,μ,ᾱ) ∈ S 4 , with S 4 denoting the 4-sphere, that is {x 2 +ȳ 2 +ε 2 +μ 2 +ᾱ 2 = 1}, and r ∈ [0, ∞). As is usual with the blow-up method [27], instead of working in spherical coordinates, we consider local coordinates in local charts. In our particular context, these local charts parametrize different hemi-spheres of S 4 . Analogous to the analysis of the canard point in [39] we consider the charts K 1 = {ȳ = 1} and K 2 = {ε = 1}. To distinguish the local coordinates in these charts, we let (r 1 , x 1 , ε 1 , µ 1 , α 1 ) be local coordinates in K 1 , and (r 2 , x 2 , y 2 , µ 2 , α 2 ) be local coordinates in K 2 . In this way, these local coordinates are defined by: for y ≥ 0, In particular, it is worth noting that in chart K 1 the coordinate r 1 is a rescaling of the "original coordinate" y for y ≥ 0, while in chart K 2 , the coordinate r 2 is a rescaling of ε ≥ 0. Furtheremore, in a qualitative sense, in chart K 1 one studies trajectories of (19) as they approach and leave a small neighborhood of the fold point in the positive y direction, while in chart K 2 one investigates the trajectories of (19) within a sufficiently small neighborhood of the fold point.
The coordinates in the above charts are related by the transition maps: for ε 1 > 0 and for y 2 > 0.

Analysis in the rescaling chart K 2
The blown-up (and desingularized) local vector field in this chart reads as where g 2 = g 2 (r 2 , x 2 , y 2 , α 2 ) is smooth and defined by the blow-up ofĝ. More precisely, from (8) and keeping in mind the usual desingularization step, one has that is the blown-up state-feedback controller to be designed. Observe that, analogously to what is described in Remark 2, we have that for r 2 = α 2 = µ 2 = 0 the orbits of (26) are given as level sets of the function Having this in mind, we are going to design µ 2 in such a way that for a trajectory (x 2 (t 2 ), y 2 (t 2 )) of (26) one has lim t2→∞ H 2 (x 2 (t 2 ), y 2 (t 2 )) = h, where h defines the desired canard cycle and t 2 denotes the time-parameter of (26).
We approach the design of µ 2 as follows: we start by restricting to {r 2 = 0} and define Next we define a candidate Lyapunov function given by and note that L 2 > 0 for all H 2 = 0 and that L 2 = 0 if and only if H 2 = 0, if and only if (x 2 , y 2 ) ∈ γ h , where by γ h we denote the reference canard cycle, that is It follows that where µ 0 2 = µ 2 (0, x 2 , y 2 , α 2 ). Naturally, we want to design µ 0 2 such that L 2 < 0, or at least L 2 ≤ 0. We now see that a convenient choice of µ 0 2 is where c 1 > 0 and c 2 ∈ R are the controller gains. Using (32) we have Note that, because the exponential function is positive, the previous inequality holds for every value of c 2 ∈ R, however a particular choice of c 2 may drastically change the performance of the controller, hence its inclusion in (32). This can be readily seen if we substitute H 2 in (32): Let D ⊂ R 3 be a bounded domain. We see that µ 0 2 is bounded for all (α 2 , x 2 , y 2 ) ∈ D. However, since c 2 appears inside the exponential, the upper bound of |µ 0 2 | can vary widely depending on the choice of c 2 . The relevance of c 2 shall be detailed in section 3.1.2. By Lasalle's invariance principle [44] we have that, under the controller (32) and r 2 = 0, the trajectories of (26) eventually reach the largest invariant set contained in Note, however that {x 2 = 0} is generically not invariant for the closed-loop dynamics (26). Indeed, the closed-loop system (26) (restricted to r 2 = 0) reads as where setting x 2 = 0 leads to (x 2 , y 2 ) = (−y 2 , 0). Therefore, we now have that all trajectories of (26) eventually reach I 2 = {(x 2 , y 2 ) = (0, 0)} ∪ H 2 = 0 . Since the origin is an equilibrium point of (36) 3 , we have that every trajectory with initial conditions (x 2 (0), y 2 (0)) ∈ R 2 \ {(0, 0)} eventually reaches the set H 2 = 0 as t 2 → ∞. With the previous analysis we have shown the following: Proposition 1. Consider (26). Then, for r 2 ≥ 0 sufficiently small a controller of the form where c 1 > 0 and c 2 ∈ R and with H 2 is as in (28), renders the orbit γ h locally asymptotically stable.
Proof. As described above, the stability of γ h for (37) is equivalent to the stability of the zero solution ofH Substituting (32) in (38) we get We have shown that for r 2 = 0, the origin is locally asymptotically stable for (39). An apparent issue in (39) is the term x 2 2 . However, we have also shown that {x 2 = 0} is not invariant. Therefore (39) is a particular case of the non-autonomous scalar equatioñ where a(t 2 ) ≥ 0 for all t 2 and a(t 2 ) > 0 for almost all t 2 (here t 2 is the time parameter in the chart K 2 ). The solution of the unperturbed equation (40) is H 2 (t 2 ) = k exp − t2 t0 a(s 2 )ds 2 , for some k ∈ R. So, due to the properties of a(t 2 ), the trivial solution of (40), with r 2 = 0, is asymptotically stable, which is preserved under sufficiently small perturbations O(r 2 ) [11].
We show in Figure 4 a simulation of the result postulated in Proposition 1. Let us emphasize at this point that designing the controller in the rescaling chart justifies using H 2 to define a convenient Lyapunov function, even if there are higher order terms in the original vector field (19). We also point out that the maximal canard becomes unbounded in this chart. Such a case shall be studied in chart K 1 (see section 3.1.2 below). Next we digress on how to deal with a certain class of higher order terms even if r 2 (equivalently ε) is not small. Lemma 1. Consider (26) with r 2 > 0 fixed and let Γ 2 ⊂ R 2 be a neighbourhood of γ h . Assume that the function g 2 satisfies 1. g 2 = x 2 φ 2 (r 2 , x 2 , y 2 , α 2 ), where φ 2 is smooth and vanishes at the origin.

Remark 5.
• If the third assumption of Lemma 1 does not hold, then trajectories converge to an equilibrium point contained in the set H 2 = 0 .
We show in Figure 5 a simulation regarding Lemma 1. . Due to the way the local coordinates in this chart are defined, choosing r 2 = 1, essentially amounts to considering ε = 1 in (19). On the left we show the orbits corresponding to ν 2 = 0, and on the right those for ν 2 given as in Lemma 1. Observe on the left that trajectories do not follow the desired canard while on the right they do. This means that the extra term ν 2 is suitable to render the canard asymptotically stable when the perturbations of order O(r 2 ) in (26) are not small.

Analysis in the directional chart K 1
We are now going to look at the controlled dynamics in the chart K 1 . This serves two purposes: the first is of giving a more precise meaning to the constant c 2 in the controller (37); the second is to corroborate that the controller designed previously is indeed able to also stabilize the (unbounded) maximal canard. Using the definition on K 1 as in (23), we have that the dynamics in this chart read as where, in particular, µ 1 denotes the controller written in the local coordinates of this chart. Since we have already designed a controller in the chart K 2 , see (37), we can use the transformation (25) to express µ 1 as where, analogous to what we have done in chart K 2 , we define Remark 6.
• If h = 0, then µ 1 becomes unbounded as ε 1 → 0 unless H 1 = 0 (previous observation). This is to be expected as, in the limit ε 1 → 0 the only canard orbit to stabilize is the maximal canard since lim ε1→0 H 1 = 0. Therefore, we are going to study the closed-loop dynamics (46) for the particular choice of h = 0 and for the limit ε 1 → 0. Our goal is to refine the constant c 2 so that µ 1 remains bounded whenever h = 0 and ε 1 → 0. Moreover, recalling that for this chart we have ε 1 = y ε , the limit ε 1 → 0 corresponds to the limit y → ∞ for fixed ε > 0. So from now on we let h = 0, that is We also restrict to {r 1 = 0}. In such a case we have and the closed loop system reads as It shall also be relevant to consider H 1 , namely First of all we note that lim ε1→0 H 1 = 0, and lim ε1→0 H 1 = 0 for c 2 < 4. Next, we focus on (49) where we observe that in order for the controller to be bounded as ε 1 → 0the constant c 2 should be less than 2. To be more precise: Proof. Straightforward computations.
From Lemma 2 we have that, to follow the maximal canard (h = 0) one must choose c 2 < 2 to ensure that the controller is bounded. Although analytically any choice of c 2 < 2 suffices, a particular choice may influence drastically numerical simulations since c 2 appears in the exponential. For instance, we see from the first line of (51) that c 2 < 2 but arbitrarily close to 2 reduces the contribution of the exponential term, which may induce issues in numerical simulations. For all other canards, c 2 ∈ R is sufficient. However, again from the computational perspective, c 2 = 2 is the appropriate choice as it eliminates the exponential term in (49) and in (51), which is rather convenient for simulations. We remark that a completely analogous analysis, which we omit for brevity, follows for the chart K 3 = {x = 1} where canards corresponding to h < 0 can be considered. The arguments and the conclusion are the same, namely, for h < 0 one should set c 2 < 2 so that the controller remains bounded along the unbounded canards.

Proof of Theorem 1
To prove Theorem 1 we first blow-down the controller µ 2 . To keep it simple we shall blow-down (37), but of course the same holds for (41). So, recall from (37) that the blown-up controller is Next, from (23) we havê whereĤ =Ĥ(x, y, ε) = 1 2 exp − 2y ε y ε −x 2 ε + 1 2 as stated in the first item of Theorem 1. Under (53) the closed-loop system corresponding to (19) reads aŝ Next, it is important to observe that lim ε→0Ĥ = 0. This means that for ε = 0 the only reference canard that is reachable is the maximal canard 4 . The maximal canard corresponds to h = 0. So, setting h = 0, and since one chooses c 2 < 2 (recall Section 3.1.2), it follows that lim ε→0 c 1x ε 1/2 exp c 2 yε −1 Ĥ = 0, meaning that the layer equation for (54) iŝ which indeed has the same type of singularity at the origin as the open-loop system, a fold. This shows that (53) is a compatible controller in the sense of Definition 4.

The slow control problem
In this section we consider the slow-control problem where the objective is, as in Section 3.1, to stabilize a prescribed canard γ h . Due to space constraints, and because the analysis is similar to the one performed in Section 3.1, we only state the relevant result.
Theorem 2. Consider (56) and letĤ = H(x, y, ε) be defined by (11). Then, the compatible controller where c 1 > 0, c 2 ∈ R and h ≤ 1 4 renders the canard orbit γ h = (x, y) ∈ R 2 | H = h locally asymptotically stable for ε > 0 sufficiently small. A convenient choice of controller gain c 2 for the maximal canard is c 2 < 2. By convenient we mean that such a choice ensures that the controller remains bounded as y → ∞.
In Figure 6 we illustrate the statement of Theorem 2.

Controlling Canard Cycles for the van der Pol oscillator
In this section we are going to extend the ideas developed previously to control canard cycles in the van der Pol oscillator. The main idea is to adapt and extend the controller proposed in Theorem 1, and to use it to control canard cycles of the van der Pol oscillator. In this context we distinguish two types of canard cycles: a) canards with head and b) canards without head. Canards with head refer to canard cycles with two fast segments, while canards without head have only one fast segment, see Figure 8. Furthermore, due to its relationship with some neuron models, like the Fitzhugh-Nagumo model [21,48], we shall consider that the controller acts on the fast variable only. The idea is that the controller represents input current. Thus, let us study Remark 7. For simplicity, we have chosen to present the case α = 0. However, the case α = 0 follows straightforwardly from considering the arguments at the beginning of Section 3.1.
The corresponding critical manifold reads as, The repelling and attracting parts of S 0 are denoted respectively by S r 0 and S a 0 , and are given by (60) Furthermore, system (58) has two fold points, one at the origin and one at (x, y) = (2, 4 3 ). In fact, the origin is a canard point and the singular limit of (58) is as shown in Figure 7.
To state our main result, let N 1 ⊂ R 2 be a region containing a subset of the repelling critical manifold S r 0 and N 2 ⊂ R 2 a small region containing a subset of S 0 around the origin. Although it is not necessary to be precise on such regions, since several choices are possible, an example of N 1 and N 2 is as follows where the defining positive constants are such that N 1 and N 2 have a non-empty intersection in the first quadrant, and 0 < y min ∈ O(ε) and y min < y h < 4 3 . The precise meaning of these bounds is given in sections 4.1 and 4.2, and is already sketched in Figure 7.
Proposition 2. Consider (58), let ψ i be a bump function with support N i , and let the repelling slow manifold S r ε be given by the graph of x = φ(y, ε). Then, one can choose N i , positive constants c 1 and k 1 , and a small constant x * , |x * | 1, such that the controller where and with stabilizes a canard cycle with height y h . Moreover, if x * < 0 then the canard is without head, while if x * > 0 then the canard is with head.
Sketch of proof: As before, all the analysis is carried-out in the blow-up space. The overall idea is as follows: the controller to be designed acts only within a small neighbourhood of {0} ∪ S r ε , mainly because the rest of the slow manifold is already stable, so there is no need of stabilizing it. The desired height of the canard is regulated by the constant y h . The controller u 2 controls the trajectories near the canard point and therefore is given by Theorem 1, where we have made the choice h = 0 and c 2 = 2. So, the new analysis is performed in the chart K 1 = {ȳ = 1}, where the objective is to stabilize the (normally hyperbolic) repelling branch of the slow manifold S r ε resulting in the controller u 1 . Later, in section 4.2, we combine the two controllers and justify the form of the controller given in the Proposition. The most important feature of u 1 is to control the location of the orbits relative to S r ε as it is precisely such location that determines the direction of the jump once the orbits reach the desired height. To avoid smoothness issues, the regions where the controllers are active are defined via bump functions. A schematic representation of this idea is provided in Figure 7, while the details of the proof follows from Sections 4.1 and 4.2.
S r 0 Figure 7: Strategy for the control design: first, within a small neighborhood of the canard point (red-shaded region), we use the controller designed in Section 3. Afterwards, a second controller is designed in chart K 1 and whose task is to stabilize the (normally hyperbolic) repelling branch S r ε . This second controller is active on a neighborhood of S r 0 (green-shaded region). Furthermore, it is via such controller that we steer the orbits towards either side of S r ε . This induces that the trajectories jump towards the desired direction once the second controller is inactive. The two orbits illustrate the aforementioned strategy.
In Figure 8 we show some simulations using the proposed controller. Before proceeding with the proof of Proposition 2, let us point out that it is straightforward to use the proposed controller to produce robust mixed mode oscillations (MMOs) [12]. One way to do this is as follows: first of all, we assume that we are able to count the number of small amplitude oscillations (SAOs) and of large amplitude oscillations (LAOs). Next, let us say that we start by following a canard without head, so we set the controller constant x * < 0 and y h to the desired height. After the number of desired SAOs has been reached, we change the controller constant x * to x * > 0 and, if desired, y h to a new height value. So, the controller will now steer the system to follow a canard with head. This process can be repeated to produce any other pattern allowed by the geometry of the van der Pol oscillator. We show in Figure 9 an example of stable MMOs that are obtained using the controller of Proposition 2.

Analysis in the directional chart K 1
Similar to the analysis in section 3.1.2 we use a directional blow-up defined by  Therefore, the local vector field associated to (58) reads as To have a better idea of what we are going to achieve with the controller, it is worth to first look at the uncontrolled dynamics. Let us define a domain Lemma 3. Consider (66) with µ 1 = 0. Then, one can choose constants ρ 1 > 0 and δ 1 > 0 such that the following properties hold within the domain D 1 .
2. Let M 1 be defined by The set M 1 corresponds to the set of equilibrium points of (66) restricted to {ε 1 = 0}. Moreover, let us denote the subsets The subset M 1,− is attracting and the subset M 1,+ can be partitioned as M 1, are the repelling and attracting branches of M 1,+ , respectively.
3. Restricted to {r 1 = 0} there exist 1-dimensional local center manifolds E 1,− and E 1,+ located, respectively, at the points p 1,− and p 1,+ . Such manifolds are given by where The flow along E 1,− is directed away from the point p 1,− and the flow along E 1,+ is directed towards the point p 1,+ . Furthermore, the center manifolds E 1,± are unique.
3. Design a variational controller that renders W cl 1 locally exponentially stable: For this it is enough to take the x 1 -component of the variational equation. So, let z 1 = x 1 −φ 1 −x * 1 . The corresponding variational equation along W cl 1 is Recall from (75) that φ 1 > 0 for r 1 ≥ 0 sufficiently small. Then, we propose to introduce in (79) a variational controller w 1 (ε 1 , z 1 ) of the form where k 1 ≥ 0. With w 1 as above, the closed-loop variational equation becomes and we readily see that, for r 1 ≥ 0 sufficiently small, z 1 → 0 exponentially as t 1 → ∞ (where by t 1 we denote the time-parameter in this chart). We also notice that the constant k 1 helps to improve the contraction rate towards W cl 1 . Moreover, since w 1 vanishes along W cl 1 , the variational controller does not change the closed-loop centre manifold W cl 1 . Finally, observe that the role of the small constant x * 1 is to shift the position of W cl 1 relative to its open-loop counterpart W 1,+ . This is important in order to tune the direction along which the trajectory jumps once the controller is inactive. 4. Restrict next to {ε 1 = 0}: Note that ν 1 (r 1 , 0, x 1 ) = 0, then we have Similar to the previous step, the new line of equilibrium points is slightly shifted according to x * 1 . In fact, the relevant set of stable equilibrium points of (82) is given as 3. If x * 1 = 0, the centre manifolds W cl 1 and W 1,+ coincide. On the other hand, if x * 1 < 0 (resp. if x * 1 > 0) then W cl 1 is located "to the left" (resp. "to the right") of W 1,+ in the x 1 -direction. 4. The image of R 1 under the flow of (66) is a wedge-like region at Σ ex 1 ∩ M cl 1 .
Proof. The proof follows directly from our previous analysis. In particular, the second item is implied by the stability properties of W cl 1 | {r1=0} , W cl 1 | {ε1=0} , and the fact that r 2 1 ε 1 = ε. The closed-loop dynamics corresponding to (66) under the controller (87) are as sketched in Figure 10. Figure 10: On the left we show the qualitative behavior of the open-loop (that is with µ 1 = 0) system (66), while on the right we show the closed-loop system obtained with the controller of Proposition 3. In both cases, the 2-dimensional surface illustrates the centre manifolds W 1,+ on the left and W cl 1 . The relative position of W cl 1 with respect to W 1,+ is determined by x * 1 . In the sketch on the right we show that W cl 1 is to the left of W 1,+ , which is indicated by the dashed curves.
To finalize this section, we blow-down the controller of Proposition 3, as it will be used in the forthcoming section. where and where φ = φ(y, ε) is defined by S r ε , that is by S r ε = {x = φ(y, ε)}.
Proof. The expression of u 1 follows from straightforward computations using (65) in (87). To check φ is as stated, note that the blow-down induces the relation {x 1 = φ 1 } ↔ x = √ yΦ(φ 1 ) = φ , 6. The sumū is, by virtue of the partition of unity, well defined as a global controller on B. Therefore, the global closed-loop vector fieldX cl :=X +ū is also well defined.
7. Let us now blow-downX cl . To be more precise, we now define the closed-loop vector field X cl on R N by Φ * (X cl ) = X cl . So, we have X cl = Φ * (X cl ) = Φ * (X +ū) = Φ * (X) + Φ * (ū) = X + Φ * (ū), where we have used the fact that the push-forward is linear [45]. Next we define u := Φ * (ū) and compute With the previous methodology we define the controller that stabilizes canard cycles of the van der Pol oscillator as u = 1 2 where u 1 is as given by Lemma 4 and u 2 as in Theorem 1, and where ψ 1 is a bump function with support N 1 containing the repelling branch S r 0 and N 2 the parabola y = x 2 around the origin. Although several choices for these neighbourhoods are possible, we recall an example given at the beginning of section 4: N 1 = (x, y) ∈ R 2 : | − y + x 2 − 1 3 x 3 | < β 1 , 0 < x < 2, y min < y < y h N 2 = (x, y) ∈ R 2 : | − y + x 2 | < β 2 , −x min < x < x max , with suitably chosen positive constants β 1 , β 2 , x min , x max , y min , y h < 4 3 . We note that one must choose 0 < y min ∈ O(ε) in order to ensure that the slow manifold S r ε is within distance O(ε) of the critical manifold S r 0 . Here y h controls the height of the desired canard cycle, therefore y h < 4 3 . The neighborhood N 1 and N 2 are sketched in Figure 7.
With the controller as in (95), and given the analysis in Section 4.1, one has that orbits of (58) passing close to the origin follow closely the repelling branch of the slow manifold S r ε up to a height determined by y h . Once orbits leave the neighborhood N 1 ∪ N 2 , they are governed by the open-loop dynamics. Finally the controller of Proposition 2 is indeed (95). We have just dropped the subscript of the constant x * 1 .

Conclusions and Outlook
In this paper we have presented a methodology combining the blow-up method with Lyapunovbased control techniques to design a controller that stabilizes canard cycles. The main idea is to use a first integral in the blow-up space to regulate the canard cycle that the orbits are to follow. Later on, we have extended the previously developed method to control canard cycles in the van der Pol oscillator. Roughly speaking this procedure follows two steps: first one needs a controller that stabilizes a folded maximal canard within a small neighborhood of the canard point. Next, one needs to stabilize the unstable branch of the open-loop slow manifold and to tune the position of the closed-loop orbits with respect to it. This is essential to determine whether the closed-loop canard has a head or not. Finally one combines such controllers by means of a partition of unity.
We have further shown that the proposed controller can be used to produce stable MMOs. Several new questions and possible extensions arise from our work, and we would like to finish this paper by briefly mentioning a couple of ideas. First of all, it becomes interesting to adapt the controllers designed here to neuron models such as the FitzHugh-Nagumo, Morris-Lecar, or Hodgkin-Huxley models. Another relevant extension is to develop optimal controllers to control canards. Although from a theoretical point of view one would be interested in arbitrary cost functionals, some particular choices might be more suitable for applications. For instance, one may one want to design minimal energy controllers. It is also not completely clear whether the strategy of combining the blow-up method and control techniques still applies as the optimal controllers may be time-dependent. Finally, the notion of controlling MMOs definitely requires further investigation, as here we have just given a simple sample of the possibilities. Thus, for example, extending the ideas of this paper to higher-dimensional fast-slow systems with non-hyperbolic points is a direction to be pursued in the future.