1 Introduction↩︎

A recurrent problem across the sciences is decision–making, which can be summarized as a rational agent choosing among alternatives so as to maximise some expected profit [1], [2], on time scales ranging from operational to strategic and political [3], [4], with applications from biology [5], [6] to economics [7] and psychology [8]. A standard way to capture the resulting bounded rational behaviour is through exponential weighting [9], which in reinforcement learning [10], [11] appear as softmax, Gibbs, or Boltzmann exploration policies, in which the probability of selecting a strategy is proportional to an exponential function of its expected reward and which balance exploitation against exploration [12].

For our research, we study a resource’s stock dynamics as being consumed by two distinguished groups following different exploitation strategies, one consuming an unlimited but highly costly common resource, such as wind energy, and the second employing a comparably cheaper but restricted and slowly renewable resource, for instance biomass. Let \(y \geq 0\) represent the limited resource stock and \(x \in[0,1]\) the share of agents exploiting it, while \(1-x\) is the portion of agents consuming the unlimited resource instead.

We model the joint evolution of \(x\) and \(y\) as the fast–slow system \[\renewcommand{\theequation}{\theparentequation.\arabic{equation}} \begin{align} \dot{x} &= \gamma_{1}(1-x)\left( \eta_{1} + \frac{1-\eta_{1}}{1+e^{-\beta_{1}(\alpha_{1} + \delta(x,y))}} \right) - \gamma_{2} x \left( \eta_{2} + \frac{1-\eta_{2}}{1+e^{-\beta_{2}(\alpha_{2} - \delta(x,y))}} \right), \tag{1}\\ \dot{y} &= \varepsilon y (1-rx), \tag{2} \end{align} \tag{3}\]

where \(0<\varepsilon\ll1\) indicates the timescale separation, involving the slow recovery speed of the limited resource stock \(y\) and the comparatively fast change of \(x\) due to the agents’ adaptation of their exploitation strategies.

The natural and metabolic component of 3 is a simple equation for \(y\) governed by the growth rate \(\varepsilon\) and each agents’ relative harvesting rate \(r\), resulting in an effective total harvesting rate of \(\varepsilon r x\). On the other hand, the economic element is a model of bounded rational behaviour governed by a set of parameters as follows. The terms \(\gamma_{1,2}>0\) represent the rate at which agents from one strategy consider switching to the opposite strategy. If an agent considers switching, they either change strategy independently of the possible profits, which happens with probability \(\eta_{1,2}\in[0,1]\) and can be called unconditional exploration, or they base their decision whether to switch strategies on the profit difference \[\delta(x,y) = y + \frac{c}{d-x}-b, \label{Eq:ProfitDifference}\tag{4}\] which happens with probability \(1-\eta_{1,2}\), and where \(b>1\) represents the benefits of harvesting one unit of the unlimited resource, while \(c>0\) and \(d>1\) are parameters governing the costs of harvesting one unit of the unlimited resource. Note that the function \(\delta(x,y)\) is chosen such that these costs decrease from \(c/(d-1)\) to \(c/d\) as the share of agents exploiting it, i.e. \(1-x\), grows. Particularly, agents following the latter strategy employ a Boltzmann (or softmax) policy to decide their future strategy [12], [13]. Therefore, the logits, i.e., the logarithm of the odds of the probability of a certain event occurring, of the relative probabilities for both strategies are proportional to the corresponding profits, multiplied by an inverse temperature parameter \(\beta_{1,2}>0\), which is a common way to model bounded rationality [14], [15]. In addition to the payoff comparison, agents favour switching to a certain extent governed by an offset parameter \(\alpha_{1,2}>0\), which can be interpreted as a conditional exploration. Altogether, considering soft optimization and conditional exploration lead to the logistic terms \(1 / (1 + e^{-\beta(\alpha\pm\delta)})\) of 3 . Similar approaches are found in economics, concerning the aggregate saving rate in a large population [16], as well as in ecology, where the cost–profit difference between distinct energy resources is expressed [17] or the rate of succession of grassland to forest as described in [18].

Beyond reinforcement learning, exponential and logit choice rules have a long history in economics and learning theory [19], [20]. In parallel, fast–slow dynamical systems exhibiting fold singularities and canard dynamics have been extensively studied from a mathematical perspective [21], [22] and in a variety of applied contexts.

A representative recent case study is the tutorial treatment of the FitzHugh–Nagumo system [23], which combines a codimension-one bifurcation analysis with a numerical determination of the canard locus in parameter space. Our fold-point diagrams (Figure 2) and singular-Hopf criticality maps (Figure 4) play an analogous role for the present non-polynomial Boltzmann-type model, identifying where folds, singular Hopf bifurcations and canard explosions occur as the cost parameters vary.

For generic initial conditions, the flow generated by 3 quickly converges to attracting equilibria of the fast dynamics 1 . Once a trajectory is sufficiently close to such equilibria, it evolves according to the slow dynamics 2 until the equilibria of the fast dynamics undergo a bifurcation, repeating the aforementioned process to either produce sustained oscillations or stabilize in an equilibrium point of 3 . For more details on fast–slow systems, see 5.1. In particular, our numerical results show that under parameter variations, the critical manifold of system 3 presents either zero, two or four fold points, see Section 2. The previous fact provides our model with the capability to produce interesting dynamics, such as those shown in Figure 1 for a case with two fold points. In the left panel, open–loop responses of system 3 are depicted, including the stabilization of trajectories (red, orange) to equilibrium points near the fold points \(F_{1,2}\), and a relaxation oscillation (blue), a response formed by alternating fast and slow segments [24]. On the other hand, the right panel exhibits actuated responses (blue, purple) of 3 for the stabilization of canard orbits (red, orange), particular solutions known for their lack of robustness towards perturbations as well as their capability to travel for considerable amounts of time along repelling regions of the critical manifold, geometric object associated to 3 which describes the boundary between the fast 1 and slow 2 dynamics when \(\varepsilon=0\) [22], [25], [26].

We now state the motivation for our work and highlight its novelty.

Motivation: In the uncontrolled system, the canard cycles associated with the fold points \(F_{1,2}\) are structurally unstable, i.e., they exist only within an exponentially small window of the harvesting rate \(r\) and are destroyed by arbitrarily small perturbations, so the criticality that drives the abrupt transitions between strategies cannot be exploited reliably. This motivates a control approach that renders a chosen canard cycle a robust attractor, enabling a controlling entity to exploit the system’s criticality, and hence to produce or suppress drastic consumption changes, with minor effort. From an applied perspective, stabilising an orbit near \(F_1\) amounts to preserving a high stock of the limited resource \(y\), for instance a public good such as a forest providing ecological habitats or recreational value, by limiting the share of agents consuming it. In contrast, stabilising near \(F_2\) amounts to consuming \(y\) in much larger but still sustainable amounts for purely economic reasons.

Novelty: Existing studies of such decision dynamics typically rely on polynomial or locally truncated models, and control objectives are seldom posed near non–generic folds. Our contribution is threefold. First, we analyse a genuinely non–polynomial, Boltzmann–type decision model within a fast–slow framework, and we locate its folds together with the associated singular Hopf bifurcations and canard explosions (Section 2). Second, and central to this work, we design compatible fast–slow controllers, that is, schemes that stabilise a target canard cycle without dramatically altering the original dynamical structure (see 5.2), around a possibly non–generic fold point. This extends the generic theory of [26] to higher contact orders \(k>1\). The relevance of non-generic folds, where higher-order terms better capture the local geometry of the critical manifold, is detailed in Section 3. Third, the proposed controllers are robust to bounded parametric uncertainty and act on the agents’ decision rules rather than on the resource dynamics directly, which makes them less intrusive, and more realistic, from an applied standpoint.

The rest of this work is organized as follows. In Section 2 we present and justify some assumptions that simplify our analysis by reducing the effective number of parameters to be considered. Next, our main results are detailed in Section 3. First, we derive compatible fast–slow controllers for a generalized quadratic system, see Theorem 1. Thereafter, we give arguments for applying our results in the decision–making system 3 and demonstrate its effectiveness to control canard cycles, even in the presence of modelling uncertainties. Finally, we discuss the reach of our controllers and potential future study directions in Section 4. In 5, we present the mathematical background necessary for our studies.

2 Preliminary analysis and motivation↩︎

System 3 has thirteen parameters \((\gamma_1,\gamma_2,\eta_1,\eta_2,\alpha_1,\alpha_2,\beta_1,\beta_2,b,c,d,r,\varepsilon)\), which adds to the complexity of its analysis. To mitigate this complication, we describe and justify the main assumptions we make about 3 to reduce the number of parameters. In addition, we provide a motivation for our main results, namely, the control of ‘small’ canard cycles. The following assumptions are made for 3 .

A1:

The parameters \(\eta_i\) represent the probability that agents switch strategy independently of payoffs. In the context considered here, namely switching between exploiting a slowly renewable but cheap resource and an unlimited but costly resource, such payoff-independent switches are expected to be rare, since strategic changes are mainly driven by relative profitability. It is therefore natural to assume \(\eta_i\) is small. Hence, for purposes of analysis, we assume that \(\eta_1=\eta_2=0\).

A2:

The offset parameters \(\alpha_i\) represent conditional exploration biases, independent of the actual payoff difference. In the absence of compelling evidence for an asymmetric bias, and to reduce parameter dimensionality, we set \(\alpha_1=\alpha_2=\alpha>0\).

A3:

The inverse temperature parameters \(\beta_i\) model how strictly agents react to payoff differences. Since both groups are drawn from the same population of decision-makers and face similar informational constraints, it is reasonable to assume a common value \(\beta_1=\beta_2=\beta>0\).

We note that the qualitative features of our results remain valid under small perturbations of A1–A3. In fact, all our simulations are performed without the above considerations to verify such a claim.

Next, let \(\gamma\mathrel{\vcenter{:}}=\frac{\gamma_1}{\gamma_2}\). After a time rescaling \(t\mapsto \frac{t}{\gamma_2}\) and recycling the notation \(\varepsilon\) for \(\frac{\varepsilon}{\gamma_2}\) we obtain the simplified model \[\begin{align} \dot{x} &= \frac{\gamma(1-x)}{1+e^{-\beta(\alpha + \delta(x,y))}} - \frac{x}{1+e^{-\beta(\alpha - \delta(x,y))}}, \\ \dot{y} &= \varepsilon y (1-rx). \end{align}\] Finally, we perform the rescaling \((\tilde{\alpha},\tilde{y}, \tilde{c}, \tilde{b})=(\beta\alpha,\beta y,\beta c,\beta b)\) to eliminate the parameter \(\beta\). For notational convenience we recycle the notation \(\alpha,y,c,b\) for the rescaled terms obtaining \[\label{eq:main2} \begin{align} \dot{x} &= \frac{\gamma(1-x)}{1+e^{-(\alpha + \delta(x,y))}} - \frac{x}{1+e^{-(\alpha - \delta(x,y))}}, \\ \dot{y} &= \varepsilon y (1-rx). \end{align}\tag{5}\] With this, we have reduced the number of parameters to seven \((\gamma,\alpha,b,c,d,r,\varepsilon)\).

The critical manifold of 5 (see the necessary background on fast-slow systems in 5) is defined as \[\label{eq:C0} \mathcal{C}_0=\left\{ (x,y)\in\mathbb{R}^2\,|\frac{\gamma(1-x)}{1+e^{-(\alpha + \delta(x,y))}} - \frac{x}{1+e^{-(\alpha - \delta(x,y))}}=0\,\right\}.\tag{6}\] This already highlights the difficulty of analyzing 5 , as the critical manifold does not have a clear explicit expression, with the exponential function playing the most important role. Hence, we combine analytical and numerical methods to describe \(\mathcal{C}_0\). The following propositions provides certain sufficient conditions describing the properties of the critical manifold, and is complemented numerically afterwards.

Intuitively, normal hyperbolicity asks that the fast flow be uniformly attracting or repelling transverse to \(\mathcal{C}_0\). Equivalently, that the fast linearisation at \(\mathcal{C}_0\) have no eigenvalue on the imaginary axis. For the scalar fast variable here this is simply \(\partial_x f\neq0\) (\(\partial_x f<0\) attracting, \(\partial_x f>0\) repelling), and it fails precisely at the folds, where \(\partial_x f=0\). The formal definition is given in 5.1 (Definition 2). Contrasting Proposition 1, the next Proposition qualitatively characterizes the non-hyperbolic scenario.

We close this section by locating the bifurcations of the full system \(0<\varepsilon\ll 1\) that organize the transition between equilibria, canard cycles, and relaxation oscillations. Although the critical manifold \(\mathcal{C}_0\) does not admit a closed-form expression, the structure of the slow equation in 5 allows us to determine the location of the Hopf bifurcations exactly. Throughout, let \[f(x,y)=\gamma(1-x)\,\sigma\big(\alpha+\delta(x,y)\big) -x\,\sigma\big(\alpha-\delta(x,y)\big)\] denote the right-hand side of the fast equation in 5 , where \(\sigma(z)=\left(1+e^{-z}\right)^{-1}\).

Proposition 3 has the following consequences, which we highlight individually:

• In contrast with the generic singular Hopf scenario, in which the bifurcation occurs at an \(\mathcal{O}(\varepsilon)\)-distance from the parameter value at which the equilibrium sits exactly on the fold [21], [22], here the bifurcation parameter value carries no \(\varepsilon\)-correction at all. This is a consequence of the fact that \(g(x,y)=y(1-rx)\) and \(\partial_y g\) vanish simultaneously at the equilibrium.

• Item 3 connects this analysis to the degenerate folds studied in Section 3. At a non-generic fold (\(k> 1\)) the transversality of the eigenvalue crossing fails, and the Hopf bifurcation itself degenerates.

• Since each fold point yields one Hopf bifurcation, the fold-point diagrams of Figure 2 also account for the singular Hopf bifurcations of the full system as \(r\) varies.

• The criticality of these Hopf bifurcations is governed by the first Lyapunov coefficient, which due to the transcendental nature of \(f\) does not admit a closed form. We therefore evaluate it numerically through the intrinsic, coordinate-free criticality invariant \(\sigma\) of [27] (with \(\sigma<0\) corresponding to a supercritical bifurcation), which avoids any normal-form transformation and depends only on Lie derivatives of the fast field at the fold. The resulting criticality maps are shown in Figure 4. At the parameters used for the controllers in Section 3 the bifurcation is supercritical over most of the two-fold region, but its sign changes across an interior curve. The corresponding mechanism by which criticality changes is out of the scope of the current presentation.

3 Main Results↩︎

In this section, we detail the local analysis and the derivation of fast–slow controllers that stabilize canard orbits in a neighborhood of a fold point of the decision–making system 3 . We note that equilibrium stabilization near a fold can be viewed as a singular limit of canard stabilization, but the latter provides substantially improved robustness with respect to parameter uncertainty and structural perturbations. Indeed, stabilizing an equilibrium at a fold by tuning a single parameter (e.g.the harvesting rate \(r\)) requires precise knowledge of the fold location: arbitrarily small parameter errors can move the equilibrium onto the repelling branch of \(\mathcal{C}_0\), immediately triggering large relaxation oscillations. By contrast, the canard controller targets the transition mechanism itself and remains effective under moderate uncertainty; moreover, it acts on the agents’ decision-making rules rather than on the resource dynamics directly, making it less intrusive from an applied perspective.

First, we study a general fast–slow system and design compatible controllers capable of stabilizing canards around a not necessarily generic fold point, results summarized in Theorem 1, and extending the application of our control methods to more degenerate cases. Although degenerate folds are nongeneric, they naturally arise as effective models when higher-order terms dominate the local geometry near a fold over finite parameter ranges. Moreover, depending on the dynamical features of the system at hand, we provide two different alternatives to accomplish the stabilization of canards, when having either a semi–actuated or a fully controlled system. Thereafter, we employ our general results to system 3 in order to stabilize canard orbits in two different fold points, entitling the controlling entity with the capability to stay near the bifurcation or produce abrupt transitions by exploiting the system’s criticality. In particular, in our first example, detailed in Section 3.3, we stabilize orbits near a fold point that maximizes the stock of renewable resource \(y\) by restricting the share of agents exploiting it. On the other hand, our second example, explained in Section 3.4, denotes an optimal exploitation strategy in which the maximum number possible of agents \(x\) benefit from consuming the limited resource \(y\) while reducing its stock to the local minimum, producing a favorable scenario both for the population and the controlling entity.

3.1 Generalization of the normal form controls↩︎

We begin by introducing a general fast–slow system in the form \[\begin{align} \begin{aligned} \dot{x} &= F(x,y, \varepsilon),\\ \dot{y} &= \varepsilon G(x,y, \varepsilon), \label{Eq:NormalFormApplicationsGeneralVectorField} \end{aligned} \end{align}\tag{9}\] where \(0<\varepsilon\ll1\) is the timescale separation parameter, while \(F:\mathbb{R}^{m}\times\mathbb{R}^{n}\times\mathbb{R}\xrightarrow{}\mathbb{R}^{m}\) and \(G:\mathbb{R}^{m}\times\mathbb{R}^{n}\times\mathbb{R}\xrightarrow{}\mathbb{R}^{n}\) are sufficiently smooth functions such that the associated critical manifold \(\mathcal{C}_{0}\) of the singularly perturbed vector field 9 has at least one point \(p=\left( x_{p}, y_{p} \right)\in\mathcal{C}_{0}\) as a not necessarily generic fold point, i.e., it satisfies the following conditions:

1. \(F(x_p,y_p,0)=0\) and \(G(x_p,y_p,0)\neq 0\) (the point lies on \(\mathcal{C}_0\) and the slow flow does not vanish);

2. \(\partial_x F(x_p,y_p,0)=\partial_x^2 F(x_p,y_p,0)=\cdots =\partial_x^{\,2k-1} F(x_p,y_p,0)=0\), \(k\in\mathbb{N}\) (the first \(2k-1\) fast derivatives vanish);

3. \(\partial_x^{\,2k} F(x_p,y_p,0)\neq 0\) (leading fast derivative non-zero, contact order \(2k\));

4. \(\partial_y F(x_p,y_p,0)\neq 0\) (nondegeneracy: \(\mathcal{C}_0\) is a graph over \(y\) at the fold).

The parameter \(k\) determines the contact order of the parabola–like shape in the critical manifold of 9 near the bifurcation point \(p\), which for a generic fold point is \(k=1\). It is precisely the parameter \(k\) that enables the use of our control techniques in applications for which the associated critical manifold \(\mathcal{C}_{0}\) is better locally approximated by a higher–order quadratic–like system (\(k>1\)) around a not necessarily generic fold point \(p\in\mathcal{C}_{0}\). For instance, in Figure 6 the critical manifold \(\mathcal{S}_{0}\mathrel{\vcenter{:}}=\left\{ (w,z)\in\mathbb{R}^{2}: 0.1w^{5} + 0.25w^{4} + 0.1w^{2} = z \right\}\) (black) is approximated by \(\mathcal{S}_{2}\mathrel{\vcenter{:}}=\left\{ (w,z) \in \mathbb{R}^{2}: 0.1w^{2} = z \right\}\) (red) and \(\mathcal{S}_{4}\mathrel{\vcenter{:}}=\left\{ (w,z) \in \mathbb{R}^{2}: 0.25w^{4} = z \right\}\) (blue). The coefficients of \(\mathcal{S}_{0}\) are illustrative and carry no physical meaning; the subscript on each manifold records the leading power retained, so that \(\mathcal{S}_{0}\) is the original, untruncated manifold, while \(\mathcal{S}_{2}\) and \(\mathcal{S}_{4}\) are the generic (\(k=1\)) and quartic (\(k=2\)) fold approximations. Only even subscripts arise because, by conditions (ii)–(iii), a fold has even contact order \(2k\). At larger scales \(\mathcal{S}_{4}\) provides the better fit, even though \(\mathcal{S}_{2}\) represents the first derivative different from zero in the expansion of 9 around the fold point located at the origin. In other words, although a generic approximation provides a better fit in an arbitrarily small neighbourhood of the fold, higher–order approximations yield a more accurate representation of the critical manifold over the finite region explored by canard trajectories. We note that while arbitrary polynomial approximations of the critical manifold (including linear combinations of different–order truncations) may improve pointwise accuracy, the control design in Section 3 relies on the Hamiltonian structure of the leading part, and therefore truncations that destroy this structure fall outside the present framework.

Hence, as long as conditions i–iv are fulfilled, the dynamics in a region of 9 sufficiently close to the fold point \(p\) are approximated by a system in the form \[\renewcommand{\theequation}{\theparentequation.\arabic{equation}} \begin{align} \dot{x} &= a_{c}x^{2k} + b_{c}y + \tilde{F}\left(x,y,\varepsilon, \alpha(x,y;\lambda)\right), \tag{10}\\ \dot{y} &= -\varepsilon \left( \sigma k a_{c} x^{2k-1} - \alpha(x,y;\lambda) + \tilde{G}\left(x,y,\varepsilon, \alpha(x,y;\lambda)\right) \right), \tag{11} \end{align} \tag{12}\] with the parabola–like shape and its stability conditions determined by the non–zero constants \(a_{c} = \frac{1}{(2k)!}\frac{\partial^{2k} F}{\partial x^{2k}}\left(x_{p}, y_{p}, 0\right)\), and \(b_{c} = \frac{\partial F}{\partial y}\left(x_{p}, y_{p}, 0\right)\), obtained through the expansion of 9 around the fold point \(p\in \mathcal{C}_{0}\), and where we reuse the same variables for the coordinate system as in 9 for simplicity, with the parameter \(\sigma = \text{sign}\left(a_{c}b_{c}\right)\). Additionally, the functions \(\tilde{F}\) and \(\tilde{G}\) collect any higher order terms, while the term \(\alpha(x,y;\lambda)\) sets the equilibrium point of the slow problem of 12 . Furthermore, if we consider that for a neighborhood sufficiently close to the fold point \(p\) the higher order terms \(\tilde{F}(x,y,\varepsilon,\alpha(x,y;\lambda))\) and \(\tilde{G}(x,y,\varepsilon,\alpha(x,y;\lambda))\) in 12 are zero, then for \(\varepsilon>0\) and \(\alpha(x,y;\lambda)=0\) orbits of 12 are given by level sets of \[H\left( x, y, \varepsilon \right) = \frac{\sigma}{2}e^{2y/\sigma \varepsilon}\left( \frac{b_{c}y}{\varepsilon} + \frac{a_{c}x^{2k}}{\varepsilon} - \frac{\sigma b_{c}}{2} \right), \label{Eq:HamiltonianGeneralVectorField}\tag{13}\] and it is known that canard cycles exist for \(H(x,y,\varepsilon)=h\) with \(0<-h/b_c<1/4\) [21], [26]; the classical range \(h\in(0,1/4)\) corresponds to the normalization \((a_c,b_c)=(1,-1)\). Thus, a region of the critical manifold of 9 sufficiently close to the fold point \(p\) can be approximated by a manifold as \[\mathcal{S}_{0} = \left\{ (x,y)\in\mathbb{R}^{2}: a_{c}x^{2k} + b_{c}y + \tilde{F}(x,y,0,\alpha(x,y;\lambda)) = 0 \right\}.\] Now, let us introduce the fast–slow control entries into 12 , as \[\begin{align} \begin{aligned} \dot{x} &= a_{c}x^{2k} + b_{c}y + \tilde{F}(x, y, \varepsilon, \alpha(x,y;\lambda)) + u(x,y;\xi),\\ \dot{y} &= -\varepsilon \left( \sigma k a_{c}x^{2k-1} - \alpha(x,y;\lambda) + \tilde{G}(x, y, \varepsilon, \alpha(x,y;\lambda)) + v(x,y;\xi)\right), \end{aligned} \label{Eq:GeneralizedQuadraticSystemHOTControl} \end{align}\tag{14}\] where \(\xi\) collects all possible controller parameters. The controlled system 14 is obtained by adding controllers to each of the equations of the normal form 12 ; as such, each controller is added to actuate at different scales; \(u\) is the fast controller (or acts on the fast variable), while \(v\) is in the slow regime. Moreover, observe that without control restrictions both \(u(x,y;\xi)\) and \(v(x,y;\xi)\) are able to compensate the higher order terms \(\tilde{F}(x,y,\varepsilon,\alpha(x,y;\lambda))\) and \(\tilde{G}(x,y,\varepsilon,\alpha(x,y;\lambda))\), respectively, when enough knowledge on these higher order terms is at hand. Hence, assuming that the higher order terms are compensated by the control entries \(u(x,y;\xi)\) and \(v(x,y;\xi)\), two possible alternatives arise in order to bring the slow dynamics of 14 to the necessary form depending on the \(\alpha(x,y;\lambda)\) term, that is, setting the equilibrium of the slow dynamics in the origin so that a canard point is obtained. First, if \(\alpha(x,y;\lambda)=\alpha\) is a constant value, then a coordinate transformation \(\hat{x} = x-\alpha\) is enough to set the equilibrium of the slow dynamics of 14 at the origin. On the other hand, if \(\alpha(x,y;\lambda)\) presents a different variable or parameter dependence, then the slow control \(v(x,y;\xi)\) is able to compensate its effect. Hence, the fast controller \(u(x,y;\xi)\) acts as the main control, while the slow input \(v(x,y;\xi)\) becomes a support control when the slow dynamics require it. Thus, after compensating the effect of the higher order terms \(\tilde{F}(x,y,\varepsilon,\alpha(x,y;\lambda))\), and \(\tilde{G}(x,y,\varepsilon,\alpha(x,y;\lambda))\), as well as \(\alpha(x,y;\lambda)\), the control problem reduces to \[\begin{align} \begin{aligned} \dot{x} &= a_{c}x^{2k} + b_{c}y + u(x,y;\xi), \\ \dot{y} &= -\varepsilon \sigma k a_{c}x^{2k-1}. \end{aligned} \label{Eq:GeneralizedQuadraticSystem} \end{align}\tag{15}\] Then, the critical manifold of 15 is \(\mathcal{C}_{0} \mathrel{\vcenter{:}}= \{ (x, y) \in \mathbb{R}^{2}: y = -(a_{c}/b_{c})x^{2k} \}\), with stability conditions given by the sign of \(a_{c}\) and \(u(x,y;\xi)\) being a compatible controller. Moreover, notice that 15 defines a conservative system with the Hamiltonian 13 .

Now, the aim of the control strategy is to stabilize one orbit of 15 , namely \(\gamma_{h}=\{ (x,y) \in \mathbb{R}^{2} : H=h \}\), related to the canard trajectory of interest. Hence, the control objective is to reduce the error existing between a level set of 13 and the desired trajectory \(\gamma_{h}\), determined by the constant \(h\). Therefore, we want to stabilize the error dynamics \(\tilde{H} \mathrel{\vcenter{:}}= H - h\) to zero. Specifically, we select \(h = -\tfrac{1}{4}\,\mathrm{sign}(b_c)\,e^{-c_c/\varepsilon}\) with \(c_c \in (0,\infty)\). Then \(-h/b_c = \tfrac{1}{4|b_c|}\,e^{-c_c/\varepsilon} > 0\), and \(-h/b_c < 1/4\) whenever \(e^{-c_c/\varepsilon} < |b_c|\), so for \(\varepsilon\) small enough \(h\) lies in the canard-cycle range \(0 < -h/b_c < 1/4\). Note that this selection of \(h\) is done in order to improve the precision by using exponential variations. Thus, it is clear that \[\dot{\tilde{H}} = \frac{e^{2y/\sigma \varepsilon}}{\varepsilon^{2}}\left( \dot{y}\left( a_{c}x^{2k} + b_{c}y \right) + \varepsilon \sigma k a_{c}x^{2k-1}\dot{x} \right). \label{Eq:HDerivativeRaw}\tag{16}\] Then, by substituting 15 in 16 yields \[\dot{\tilde{H}} = \frac{e^{2y/\sigma \varepsilon}}{\varepsilon}\left(\sigma k a_{c}x^{2k-1}u \right). \label{Eq:HDerivativeControl}\tag{17}\] Subsequently, an alternative to guarantee the local asymptotic stability of the desired level set \(\gamma_{h}\) is to define the error dynamics as \(\dot{\tilde{H}} = -A_{c}\tilde{H} = -A_{c}(H - h)\), for some \(A_{c} > 0\). Thus, the fast control entry is given by \[u(x,y;\xi) = -\frac{\varepsilon B_{c} x}{\sigma k a_{c}}(H-h)e^{-2y/\sigma \varepsilon}, \label{Eq:FastControlEntry}\tag{18}\] where we have desingularized the origin by setting \(A_{c} = B_{c}x^{2k}\), positive for every \(x\in\mathbb{R}\backslash\{0\}\), and \(B_{c}>0\) is the controllers’ gain.

To show stability, we define a candidate Lyapunov function as \[L(x, y) = \frac{1}{2}\tilde{H}^{2}. \label{Eq:LyapunovGeneralControl}\tag{19}\] Observe that 19 is positive for every \(\tilde{H}\neq 0\), and that \(L=0\) if and only if \(\tilde{H}=0\), if and only if \((x, y)\in \gamma_{h}\), which we recall is the control target defined earlier. Then, it is easy to check that \[\dot{L} = \tilde{H} \dot{\tilde{H}} = -B_{c}x^{2k}\tilde{H}^{2}. \label{Eq:LyapunovFastControlDerivative}\tag{20}\] Finally, to demonstrate the asymptotic stability of \(\gamma_{h}\) as 20 is only negative semidefinite, by LaSalle’s invariance principle [28], trajectories of the fast–slow vector field 15 under the control action \(u(x,y;\xi)\) as 18 reach, in finite time, the largest invariant set contained in \[\mathcal{I} = \left\{ (x, y)\in \mathbb{R}^{2}:\dot{L} = 0 \right\} = \{ x=0 \} \cup \left\{ \tilde{H} = 0 \right\}. \label{Eq:LyapunovInvariantSet}\tag{21}\] It is important to mention that, as long as \(y\neq0\), the vector field 15 does not vanish when it reaches 21 , allowing the usage of the fast control 18 in the quadratic–like system 15 . Notice, however, that \({x=0}\) is generically not invariant for the closed–loop dynamics 15 . Actually, setting \(x=0\) in 15 reduces to \(\left( \dot{x}, \dot{y} \right) = \left( b_{c}y, 0 \right)\), inducing an increment or decrement in the share of agents \(x\) depending of the values of \(b_{c}\). Hence, trajectories of 15 eventually reach \(\mathcal{I} = \{ (x,y) = (0,0)\} \cup \{ \tilde{H} = 0 \}\), and since the origin is an unstable equilibrium point of 15 with the controller as 18 , we have that every trajectory with initial condition different from zero eventually reaches the set \(\{ \tilde{H} = 0\}\) as \(t\xrightarrow{}\infty\). In Figure 7 we present an example of the effect that the fast control entry 18 has in system 15 . The target level set \(\gamma_{h}\) and the resulting trajectory are shown in solid red and blue, respectively. As it can be appreciated, the fast control scheme makes the system to behave in the desired manner, causing the orbit to follow unstable branches of the critical manifold \(\mathcal{C}_{0}\) for time of order \(\mathcal{O}(1)\). Hence, canard orbits originating from a fold point in system 15 are stable when implementing an adequate fast control law as 18 .

We summarize these results in the following theorem.

The proof of Theorem 1 directly follows the previously presented analysis and, as long as the higher order terms \(\tilde{F}(x,y,\varepsilon,\alpha(x,y;\lambda))\) and \(\tilde{G}(x,y,\varepsilon,\alpha(x,y;\lambda))\) of 12 are sufficiently small, the control scheme synthesized in Theorem 1 produces canard cycles in a neighbourhood of the fold point \(p=(x_{p}, y_{p})\in\mathcal{C}_{0}\) of 9 . This fact shows the true strength of our control technique since, regardless of the complexity of system 9 , as long as its critical manifold satisfies the previously stated conditions, it is possible to stabilize canard orbits in a neighborhood of the fold point \(p\), even if 9 is not explicitly given in normal form. Finally, it is noteworthy to mention that as the amplitude of the target level set \(\gamma_{h} = \{ (x,y)\in\mathbb{R}^{2}:H=h \}\) increases, the approximation error between the associated critical manifold \(\mathcal{C}_{0}\) and the target level set \(\gamma_{h}\) will accordingly increase, producing an undesired response in 9 caused by the disparity of both trajectories due to the effect of the higher order terms. To improve the aforementioned behaviour, a possible alternative is the use of complementary control schemes acting sufficiently far away from the bifurcation such that when the trajectory travels in a vicinity that is correctly approximated by our expansion the system is under the effect of our control scheme, but once the expansion error increases beyond a tolerance threshold, an additional controller redirect the solution in order to follow the flow along the critical manifold \(\mathcal{C}_{0}\).

Next, we discuss on the effectiveness of the fast–slow control technique developed in this section for the control of canard orbits near a fold point in the presence of parametric uncertainties.

3.2 Control of canards in the perturbed normal form↩︎

Consider the control problem 15 under the effect of parametric perturbations in the form \[\begin{align} \begin{aligned} \dot{x} &= \left(a_{c} + \delta_{a}\right)x^{2k} + \left( b_{c} + \delta_{b} \right)y + u(x,y;\xi),\\ \dot{y} &= -\varepsilon \sigma k\left(a_{c} + \delta_{a}\right)x^{2k-1}, \label{Eq:PerturbedNormalForm} \end{aligned} \end{align}\tag{22}\] where \(\delta_{a}\) and \(\delta_{b}\) represent sufficiently small modelling uncertainties in 15 , with the same parameter interpretation as before. Hence, orbits of 22 are given by level sets of \[H_{p}\left( x, y, \varepsilon \right) = \frac{\sigma}{2}e^{2y/\sigma \varepsilon}\left( \frac{\left(b_{c}+\delta_{b}\right)}{\varepsilon}y + \frac{\left(a_{c}+\delta_{a}\right)}{\varepsilon}x^{2k} - \frac{\sigma \left(b_{c}+\delta_{b}\right)}{2} \right), \label{Eq:HamiltonianPerturbedNormalForm}\tag{23}\] with the subscript \(p\) standing for perturbation in \(H_{p}(x,y,\varepsilon)\). Therefore, our control objective is to stabilize canard orbits given by level sets of 23 in a vicinity of the fold point at the origin of 22 , by implementing the compatible fast–slow controller 18 , designed in the unperturbed normal form 15 and without explicit information regarding the perturbations \(\delta_{a}\) and \(\delta_{b}\). We begin by defining the error \(\tilde{H}_{p} = H_{p} - h\), with the associated perturbed error dynamics given by \[\dot{\tilde{H}}_{p} = \frac{e^{2y/\sigma \varepsilon}}{\varepsilon} \sigma k \left( a_{c} + \delta_{a} \right)x^{2k-1}u, \label{Eq:HDerivativePerturbed}\tag{24}\] where we follow the same procedure as for the unperturbed case, aiming to stabilize the error dynamics \(\dot{\tilde{H}}_{p}\) to zero. By substituting the fast–slow controller 18 in 24 yields \[\dot{\tilde{H}}_{p} = -\frac{a_{c}+\delta_{a}}{a_{c}}B_{c}x^{2k}(H-h), \label{Eq:HDerivativePerturbedControl}\tag{25}\] with \(B_{c}>0\) the controller gain, while \(H(x,y,\varepsilon)\) is 13 . Moreover, since the perturbations \(\delta_{a}\) and \(\delta_{b}\) are assumed to be sufficiently small, we rewrite 25 as \[\dot{\tilde{H}}_{p} = -\tilde{B}_{c}x^{2k}\left( \tilde{H}_{p} - H_{\delta}\right), \label{Eq:HDErivativePerturbedControl2}\tag{26}\] where \(\tilde{B}_{c}>0\) is the new controller gain, and \[H_{\delta}(x,y,\varepsilon) = \frac{\sigma}{2}e^{2y/\sigma \varepsilon}\left( \frac{\delta_{b}}{\varepsilon}y + \frac{\delta_{a}}{\varepsilon}x^{2k} - \frac{\sigma \delta_{b}}{2} \right), \label{Eq:HamiltonianPerturbation}\tag{27}\] since the unperturbed Hamiltonian 13 is equal to the difference between the perturbed Hamiltonian 23 and 27 , i.e. \(H(x,y,\varepsilon) = H_{p}(x,y,\varepsilon) - H_{\delta}(x,y,\varepsilon)\). Observe from 26 that the error dynamics \(\dot{\tilde{H}}_{p} = -\tilde{B}_{c}x^{2k}\tilde{H}_{p} + \tilde{B}_{c}x^{2k}H_{\delta}\) are bounded from above by \(-\tilde{B}_{c}x^{2k}\tilde{H}_{p}\). Therefore, in order for the perturbed error to converge to zero in finite time, it is sufficient to show that 27 is negative, and restricted to the perturbed critical manifold \(y = -\sigma|(a_{c}+\delta_{a})/(b_{c}+\delta_{b})|x^{2k}\), the aforementioned condition is satisfied as long as \[\sigma \delta_{a} < \delta_{b}\left| \frac{a_{c}}{b_{c}} \right|, \label{Eq:ConditionPerturbationStability}\tag{28}\] for \(\delta_{b}>0\). In Figure 8, we show examples for the stabilization of canard cycles in the perturbed normal form 22 by implementing the fast–slow controller 18 obtained from the unperturbed normal form 15 . Considering \(\delta_{b}\) positive, in the left and right columns we present the cases for which \(\delta_{a}\) is negative and positive, respectively, such that condition 28 is satisfied. In the upper row the phase portraits are depicted, where the critical manifold and response associated to the unperturbed problem 15 are represented in dotted black and blue, while the critical manifold and solution associated to the perturbed case 22 are shown in solid purple and green, correspondingly. Moreover, the target level set \(\gamma_{h}=\{ (x,y)\in \mathbb{R}^{2}:H=h \}\) is again presented in solid red, and observe that the stabilization of the canard trajectory is achieved for both cases. Notice that, under condition 28 , the perturbations \(\delta_{a}\) and \(\delta_{b}\) cause a phase shift between the responses of 15 and 22 , effect that can be appreciated in the comparison of the slow dynamics for both systems presented in the lower row of Figure 8. However, such difference does not modify the stabilization of the desired periodic pattern in the perturbed normal form 22 .

In what follows, we employ the results summarized in Theorem 1 to stabilize canard cycles in a neighborhood of a fold point of the main model 3 , even in the presence of modelling uncertainties when condition 28 is satisfied.

3.3 Control of canards in the decision–making system↩︎

Let us recall that under assumption A1-A3, the original system 3 is reduced to the simpler one 5 . For convenience let us recall that 5 reads as \[\renewcommand{\theequation}{\theparentequation.\arabic{equation}} \begin{align} \dot{x} &= \gamma(1-x)\left( \frac{1}{1+e^{-(\alpha + \delta(x,y))}} \right) - x\left( \frac{1}{1+e^{-(\alpha-\delta(x,y))}} \right),\tag{29}\\ \dot{y} &= \varepsilon y (1-rx),\tag{30} \end{align} \tag{31}\]

In particular, the choice \(\eta_{1}=\eta_{2}=0\) (Assumption A1) allows for the interesting dynamics of 3 to occur in the whole domain \(x\in[0,1]\). Nevertheless, we emphasize that every pattern observed for \(\eta_{1} = \eta_{2} = 0\) is also produced for any \(\eta_{1}, \eta_{2} \in (0,1)\) but in a reduced region of the domain of \(x\). By performing an asymptotic analysis of system 3 , we find the location of the left and right asymptotes of the associated critical manifold \(\mathcal{C}_{0}\), being \(x_{L} = (\gamma_1 \eta_1)/(\gamma_{1}\eta_{1} + \gamma_2)\), and \(x_{R} = \gamma_1/(\gamma_{1}+\gamma_{2}\eta_{2})\), when \(y\rightarrow -\infty\) and \(y\rightarrow \infty\), respectively. Hence, when \(\eta_{1} = \eta_{2} = 1\), the dynamics of 3 occur along a vertical line centered at \(x = \gamma_{1}/(\gamma_{1} + \gamma_{2})\), dramatically reducing the possible behaviours to be observed. In Figure 9 we show the effect of varying the unconditional exploration rates \(\eta_{1,2}\) when fixing the rest of the parameters in system 3 . On the left panel, we consider \(\eta_1 = \eta_2 = 0\), which represents a scenario where, once the agents have decided to switch from one strategy to another with rates \(\gamma_{1,2}>0\), they change in a completely interested manner with respect to the cost–benefit determined by 4 . In contrast, in the right panel we show the case when the agents modify their strategy following a more informed conviction, i.e., \(\eta_{1,2}\in(0,1)\). Notice that, by increasing the unconditional exploration rates, the presence of abrupt changes in the share of agents following an exploitation policy is considerably reduced without the need of an external controller. However, although desirable, this scenario is rather unreliable as it completely depends on the agents’ goodwill, motivating further the need of an outer controlling entity. Lastly, when every agent switch in an entirely informed way a complete balance between the two groups of agents exploiting each strategy is reached, allowing for the renewable resource only to increase, which would be the case when \(\eta_{1,2}=1\), and the share of agents exploiting the limited resource \(y\) would be constant with value \(x=\gamma_{1}/(\gamma_{1}+\gamma_{2})\).

In addition, we have identified in Section 2 that the critical manifold of 5 (and hence of 3 ) can have \(0\), \(2\), or \(4\) fold points. This motivates us to use the control approach of Theorem 1 to stabilize canards close to the fold points. For practical purposes, we will focus on cases where the critical manifold has exactly two fold points.

Therefore, by introducing control inputs in system 31 yields \[\begin{align} \begin{aligned} \dot{x} &= \gamma(1-x)\left( \frac{1}{1+e^{-(\alpha + \delta(x,y))}} \right) - x\left( \frac{1}{1+e^{-(\alpha-\delta(x,y))}} \right) + u(x,y;\xi),\\ \dot{y} &= \varepsilon\left( y (1-rx) + v(x,y;\xi) \right), \label{Eq:SystemControl} \end{aligned} \end{align}\tag{32}\] where \(u(x, y;\xi)\) and \(v(x,y;\xi)\) (with \(v(x,0;\xi)=0\), see Remark 7 below) represent the fast and slow control components, respectively, while \(\xi\) collects the possible control parameters involved.

Note that the selection of the controllers \(u(x,y;\xi)\) and \(v(x,y;\xi)\) in 32 is inspired from different physical applications [29]–[32]. Particularly, for our decision–making model 32 , the fast control \(u(x,y;\xi)\) represents a direct action modifying the rate at which agents change their strategy in order to adjust the stock of renewable resource \(y\) in a desired manner, for instance some regulation aimed to increase or decrease such concentration by limiting the number of agents that are allowed to exploit either resource, forcing some to rapidly switch to the other consumption strategy. By contrast, the slow control \(v(x,y;\xi)\) corresponds to actions directly favouring the recovery of the resource’s stock, for instance, in the form of infrastructure expansions both for the storage and production of the resource, as such tasks usually require extended time intervals in order to be effectively introduced or executed.

Therefore, now we present the consequent purely fast and combined fast–slow controllers for system 32 , following the results stated in Theorem 1. First, we consider the fast control scheme, i.e. \(v(x,y;\xi)=0\), and

\[u(x,y;\xi) = - \frac{\varepsilon B_{c}}{\sigma k a_{c}}\left( x-x^{*} \right)(H-h)e^{-2\left(y - y^{*}\right)/\sigma \varepsilon}, \label{Eq:FastControlMainProblem}\tag{33}\] with \[H\left( x, y,\varepsilon \right) = \frac{\sigma}{2}e^{2\left(y - y^{*}\right)/\sigma \varepsilon}\left( \frac{b_{c}}{\varepsilon}\left( y-y^{*} \right) + \frac{a_{c}}{\varepsilon}\left( x-x^{*} \right)^{2k} - \frac{\sigma b_{c}}{2} \right), \label{Eq:HamiltonianFunction}\tag{34}\] where, as explained at the beginning of this section, \(a_{c}\) and \(b_{c}\) are constant values obtained through the expansion of 32 near the fold point \(p=(x_{p}, y_{p})\in\mathcal{C}_{0}\) of interest and are responsible of setting the parabola–like shape as well as its stability properties, while \(\sigma = \text{sign}(a_{c} b_{c})\), \(k\in\mathbb{N}\), \(B_{c} > 0\) is the controller gain, and \((x_{p}, y_{p}) = \left( x^{*}, y^{*} \right)\), numerically obtained in our studies, are the coordinates of the fold point \(p\in\mathcal{C}_{0}\) of interest in the original coordinated system of 32 , required to displace the point \(p = (x_p,y_p)\) to the origin in the coordinated system of the normal form 12 . Moreover, notice that since the slow dynamics are not actuated, we are going to stabilize canards centered at \((x, y) = \left(1/r - x^{*}, 0 \right)\). Hence, it is convenient to define the coordinate transformation \(\Hat{x} = x - \left(1/r -x^{*} \right)\), which brings 33 and 34 to \[u(x,y;\xi) = -a_{c}\left( x-x^{*} \right)^{2k} + a_{c}\left( x-1/r \right)^{2k} - \frac{\varepsilon B_{c}}{\sigma k a_{c}}\left( x-1/r \right)(H-h)e^{-2\left(y - y^{*}\right)/\sigma \varepsilon}, \label{Eq:FastControlMainProblemFinal}\tag{35}\] and \[H\left( x, y,\varepsilon \right) = \frac{\sigma}{2}e^{2\left(y - y^{*}\right)/\sigma \varepsilon}\left( \frac{b_{c}}{\varepsilon}\left( y-y^{*} \right) + \frac{a_{c}}{\varepsilon}\left( x-1/r \right)^{2k} - \frac{\sigma b_{c}}{2} \right). \label{Eq:HamiltonianFunctionFinal}\tag{36}\]

As a comparison to the resulting purely fast controller 35 and Hamiltonian 36 , now we consider a fully actuated scenario and obtain a joint fast–slow controller for the decision–making system 31 , in the form \[\begin{align} \begin{aligned} u(x,y;\xi) &= -\frac{\varepsilon B_{c}}{\sigma k a_{c}}\left( x-x^{*} \right)\left( H-h \right)e^{2\left( y-y^{*} \right)/\sigma \varepsilon},\\ v(x,y;\xi) &= -\left( 1-rx^{*} \right), \label{Eq:FastSlowControlMainProblem} \end{aligned} \end{align}\tag{37}\] with \(H(x,y,\varepsilon)\) as 34 and the same variable interpretations as before. Notice that in the joint fast–slow scenario 37 , the slow component \(v(x,y;\xi)\) effectively sets the fold point \(p\in\mathcal{C}_{0}\) of 32 at the origin in the coordinated system of the normal form 12 , and therefore the additional translation \(\hat{x} = x-(1/r - x^{*})\) necessary for the purely fast controller 35 is no longer needed.

In Figure 10, we compare the results obtained by implementing the fast control 33 (left), and the fast–slow scheme 37 (right) in system 31 for the stabilization of canards in a vicinity of the leftmost fold point \(F_1\) (see Figure 1). This particular scenario represents a situation in which the controlling entity aims to increase the slowly renewable resource \(y\) stock by limiting the share of agents consuming it. Observe that both controllers effectively cause trajectories traveling near a neighborhood of the fold point to follow the targeted level set, namely \(\gamma_{h}\), in red, thus moving sufficiently close to unstable branches of the critical manifold and producing sustained canard orbits even when the initial conditions are set near the unstable branch of the critical manifold, showing the effectiveness of our control schemes. Particularly, notice the presence of a shift between the targeted orbit and the actual response in the purely fast controller due to the translation \(\Hat{x} = x - (1/r - x^{*})\). Moreover, observe that in both cases the controllers are active only during the direction change in the limited resource stock \(y(t)\), demonstrating a remarkable energetic efficiency. Now, it is clear why stabilizing the consumption ratio in the fold point \(F_1\) is desirable, however, in order to precisely stabilize the point \(F_1\), the controlling entity would need to accurately know the value of each parameter of the system, something that is rarely, if not never seen, in real–world applications. Hence, the controlling authority, for instance a government, aims to robustly stabilize a neighboring trajectory that stays near the fold point \(F_1\), hindering any abrupt transition possible even in the presence of parametric uncertainties. To achieve the aforementioned, we extend our analysis for the stabilization of canard cycles in the perturbed normal form 22 , detailed in Section 3.2, and the results are depicted in Figure 11, where a canard cycle in a vicinity of the leftmost fold point \(F_{1}\) is stabilized even in the presence of modelling perturbations by implementing the fast–slow control 37 in 32 . In particular, for this example we consider perturbations \(\delta_{\beta}=0.05\) and \(\delta_{\gamma}=-0.001\), related to parameters \(\beta\) and \(\gamma\), respectively. In the phase portraits shown in the upper row, we represent the critical manifold and response for the unperturbed system 32 in dotted black and blue, while the critical manifold and response for the same system considering the perturbations \(\delta_{\beta}\) and \(\delta_{\gamma}\) are presented in solid purple and green, correspondingly. Additionally, the target level set \(\gamma = \{ (x,y)\in \mathbb{R}^{2}: H=h \}\) is presented in red, where \(H(x,y,\varepsilon)\) is the Hamiltonian 34 of the unperturbed model 32 . Observe that the canard orbit is effectively stabilized for both the unperturbed and perturbed scenarios. Specifically, the parameters of the parabola–like shape obtained through the expansion of 32 around the leftmost fold point \(F_{1}\) for the unperturbed case are \(a_{c}=1.6423\) and \(b_{c}=0.100927\), while for the perturbed scenario are \(a_{c}=1.75871\) and \(b_{c}=0.108357\), with the fold point \(F_{1}\) of the perturbed case centered at \((x^{*}, y^{*}) = (0.6025, 28.596378)\), once again numerically identified. Therefore, the resultant perturbations, given by the parameter difference, are \(\delta_{a}=0.116412\) and \(\delta_{b}=0.00743007\), satisfying condition 28 . Additionally, in the lower row the same phase–shift observed for the perturbed normal form 22 is appreciated in the decision–making system 32 under perturbations. Nevertheless, the qualitative behaviour of 32 with and without perturbations are similar. Finally, observe that the origin of the target level set \(\gamma_{h}\) is centered at the leftmost fold point \(F_{1}\) of the unperturbed scenario as the controller is designed only with information of the unperturbed case, granting the controlling entity with the capability to stabilize orbits in a vicinity of the desired fold point even with a bounded level of parametric uncertainty.

In the following section, we further explore the capabilities of our fast–slow controllers 35 and 37 to stabilize canard cycles in a vicinity of the rightmost fold point \(F_2\) of the decision–making model 31 and discuss on its utility in our system.

3.4 Control of the optimal consumption strategy↩︎

Lastly, we extend the implementation of our control schemes in the decision–making model 31 to the stabilization of a canard trajectory in the vicinity of the rightmost fold point, namely \(F_{2}\) (see Figure 1). The selection of \(F_2\) is related to an optimal exploitation strategy as this point represents the limit at which the largest group of agents benefit from consuming the slowly renewable resource \(y\), while also being the point at which the minimum stock of resource \(y\) is reserved, advantageous from an economic perspective since it considerably reduces the storage costs for the controlling entity. On top of that, the controlling authority is entitled once more with the capability to produce abrupt transitions by exploiting the system’s criticality, for instance in a scenario in which now a greater concentration of resource in needed. As an example of the implementation of our fast–slow controllers, in Figure 12 the effective stabilization of canard cycles in a vicinity of the rightmost fold point \(F_2\) is shown along with the resulting resource, agents and control responses when using the fast–slow control scheme 37 . Observe that the variation of agents consuming the renewable resource is rather small, while the element \(y\) also presents periodic oscillations with small amplitude. In addition, notice that the controller only activates periodically and remains bounded between two considerably small values, which translates in minor actions done by the controlling entity in order to stay around the desired trajectory and optimally consume the renewable resource \(y\). Finally, in Figure 13 we show a scenario in which the controlling authority effectively regulates at will the dynamical behaviour of 31 by activating our fast–slow control scheme 37 . On the upper row of Figure 13, we show the stabilization of an orbit near the leftmost fold \(F_1\), which corresponds to the scenario when the share of agents following each consumption strategy is mostly balanced, leading to the largest amount of renewable resource \(y\) in stock, without considering the right asymptote. From an authority perspective, stabilizing a canard cycle in a vicinity of \(F_1\) represents an advantageous position in order to increase the amount of limited resource, to then allow its consumption by a larger share of agents after a specific desired time, corresponding to the free dynamics of 31 . On the other hand, in the middle row of Figure 13 we present the effect of activating our controller after a given time for the stabilization of a canard cycle in a vicinity of the fold point \(F_2\), which corresponds to the aforementioned optimal consumption strategy. Similarly, in the lower row of Figure 13 we show a sequential activation of our fast–slow controllers in order to generate a pattern that oscillates in a vicinity of the leftmost fold \(F_{1}\), then deactivate the controller and let the system travel freely near the critical manifold, to later activate the controller once more when the response is in a vicinity of the rightmost fold \(F_{2}\) to stabilize the canard orbit around a neighboring target level set, deactivate again the controller and finally activate it again in a vicinity of the leftmost fold point \(F_{1}\), showing the capability of the controlling entity to exploit the systems’ criticality to reach the desired state by sequentially activating and deactivating the controller 37 . Notice that the sharp transition in the stock of renewable resource \(y\), occurring while changing the target orbit from \(F_{2}\) to \(F_{1}\), is due to the own flow of the open–loop problem 31 . Finally, observe that the oscillation frequency for the canard orbits stabilized near the fold point \(F_1\) is comparably larger with respect to the one of the canard orbits stabilized in a neighborhood of the fold point \(F_2\). This fact is due to the form of the critical manifold \(\mathcal{C}_0\), as the resulting trajectory mainly travels along the slow direction near the fold point \(F_2\), while in the case for a canard cycle close the fold \(F_1\) the systems’ response constantly alternates between the fast and slow directions.

4 Conclusion and Discussion↩︎

Decision–making represents a fundamental component in several real world phenomena as it describes a common problem in which agents need to find the maximum benefit according to different external and internal factors. In this work, we have studied a fast–slow dynamical system describing the decision–making process that two groups of clearly identified agents undergo when selecting between two harvesting strategies, one representing the exploit of an unlimited but highly costly common source, and the other describing the consumption of a comparably cheaper but limited and slowly renewable resource. With this in mind and by an extensive numerical study, we determine parameter combinations such that the associated critical manifold of the fast–slow dynamical system analyzed presents either zero, two or four fold points, allowing for the generation of canard cycles due to the system’s properties. Such special trajectories are well–known for following unstable regions of the slow–flow for considerable amounts of time, convenient in our context for the controlling entity to exploit the systems’ criticality in order to manipulate the renewable resource dynamics according to a desired outcome and enabling the possibility to produce abrupt strategic transitions in the consumption schemes that the agents follow. Therefore, we have been able to stabilize canard orbits in our system by implementing fast–slow controllers developed through the analysis of a canonical form of a not necessarily generic fold point, and shown the effectiveness of our approach for two particular scenarios when considering a parameter combination with two fold points. The first one being a situation in which the controlling authority needs to guarantee a high concentration of renewable resource for a specific period, task achieved by implementing policies that restrict the number of agents consuming such source and enabling the controlling entity to cause an abrupt transition to a scenario in which most of the population consume the limited resource with a minor control effort. On the other hand, the second implementation describes a circumstance in which, lead by the controlling entity, the population follows an optimal consumption strategy that allows the vast majority of the entire community to consume the renewable resource while reserving a local minimum of it, which represents an economic saving considering the stock’s storage. Moreover, we demonstrate that our fast–slow controllers present a certain degree of robustness to parametric perturbations, which grants the controlling entity with the capability to stabilize canard orbits even in the presence of modelling uncertainties.

However, we want to emphasize that our results are able to stabilize canard orbits in any fast–slow system that presents planar folded canard orbits, regardless on the degeneracy degree such bifurcation points may present, allowing for a wide variety of real–world processes to benefit from our controllers. For instance, consider competition, which represents a key element in the study of diverse phenomena, ranging from biological systems [33], epidemics [34], ecology [17], [35], economics [16], [36], and social models [37], describing the interaction between populations of a broad nature, behaviour that considerably resembles the decision–making process that we have studied.

Our analysis provides several future research directions:

1. Coupled decision–making units. When several populations of the form 5 are coupled, through a shared resource stock or social interaction between their shares \(x_i\), the natural question is the synchronisation of their canard–mediated transitions, and whether the controllers of Section 3 can enforce or break it. The target canards need not be of the same type: coordinating units at the same fold differs from a crossed configuration (\(F_1\)–\(F_2\)), where the coupled canards have different amplitude and frequency.

2. Multiple limited resources. The slow space of 5 is one–dimensional because a single limited resource is tracked; strategies enter only through the share \(x\). With several limited, slowly renewable resources, each contributing its own stock, the slow space becomes higher–dimensional, providing the possibility of folded–node canards and mixed–mode oscillations [25], [38]. Extending Theorem 1, which depends on the scalar Hamiltonian 36 , to this setting, where no such Hamiltonian is available, is open.

3. Stochastic canard control. Bounded–rational exploration makes noise intrinsic to the decision dynamics. As the canard window has exponentially small width \(\mathcal{O}(e^{-K/\varepsilon})\), weak noise can dominate it and trigger unintended transitions; whether the feedback 18 survives such forcing or must be redesigned is open.

4. Mechanism of the change of criticality. Figure 4 shows the criticality invariant \(\sigma\) of the Hopf bifurcation changing sign, separating a supercritical from a subcritical onset of oscillations. In turn, this determines whether a policy past the bifurcation is reversible. Identifying the mechanism, a Bautin point or a fold collision, and locating the codimension–two points analytically is open, since \(f\) admits no closed–form Lyapunov coefficient.

5. Robust global control and fold identification. Our controllers are local, assuming known fold coordinates \((x^{\ast},y^{\ast})\) and an accurate parabola–like approximation. A switching scheme incorporating a manifold–following controller once the local error grows, and the online estimation of \((x^{\ast},y^{\ast})\) for adaptive stabilisation, are natural extensions.

5 Background on slow-fast systems↩︎

In this section, we present the theoretical principles that represent the basis of our results. We begin by briefly discussing the theory of fast–slow dynamics and geometric singular perturbation theory, specifically Fenichel’s theorem [22], [39]. Thereafter, we describe planar folded canard trajectories, emphasizing their normal form [21], [26].

5.1 Fast–slow systems↩︎

Consider the class of systems given by \[\begin{align} \begin{aligned} \varepsilon \dot{x} &= f(x,y,\varepsilon),\\ \dot{y} &= g(x,y,\varepsilon), \label{Eq:FastSlowSystem95Slow} \end{aligned} \end{align}\tag{38}\] typically known as a fast–slow vector field, composed by a set of singularly perturbed ordinary differential equations, where the over–dot denotes the derivative with respect to the slow–time \(\tau\), and with \(x=x(\tau)\in\mathbb{R}^{m}\), \(y=y(\tau)\in\mathbb{R}^{n}\) the fast and slow variables, respectively. Moreover, the parameter \(0<\varepsilon\ll1\) describes the timescale separation and the mappings \(f:\mathbb{R}^{m}\times \mathbb{R}^{n} \times \mathbb{R} \xrightarrow{} \mathbb{R}^{m}\) and \(g: \mathbb{R}^{m} \times \mathbb{R}^{n} \times \mathbb{R}\xrightarrow{} \mathbb{R}^{n}\) are assumed to be sufficiently smooth. Additionally, by defining the fast time variable \(t\mathrel{\vcenter{:}}= \tau/\varepsilon\), it is possible to express the equivalent fast–time form of 38 as \[\begin{align} \begin{aligned} x' &= f(x,y,\varepsilon),\\ y' &= \varepsilon g(x,y,\varepsilon), \label{Eq:FastSlowSystem95Fast} \end{aligned} \end{align}\tag{39}\] with the prime denoting the derivative with respect to \(t\). Now, one important mathematical theory usually employed in the analysis of systems 38 –39 is geometric singular perturbation theory (GSPT), with an overall idea of analyzing both problems separately and looking for invariant objects that can persist certain small perturbations [22]. Specifically, by considering the singular limit \(\varepsilon = 0\) in 38 and 39 , two different problems arise \[\begin{align} & 0=f(x,y,0), \qquad \qquad \qquad && x'=f(x,y,0),\\ & \dot{y} = g(x,y,0), \qquad \qquad \qquad && y'=0, \label{Eq:SingularLimits} \end{align}\tag{40}\] known, correspondingly, as the reduced slow subsystem, being a constrained differential equation [40] or a differential algebraic equation [41], and the layer problem. Naturally, the aforementioned problems are no longer equivalent but they are intrinsically related through the geometric object known as the critical manifold.

Moreover, one highly relevant property of fast–slow systems to analyze the different types of possible bifurcations is the normal hyperbolicity, defined as follows.

A useful consequence of normal hyperbolicity is that the dynamics of the constraint equation in 40 can be reduced to the dynamics of an ODE. Namely, due to normal hyperbolicity, the critical manifold \(\mathcal{C}_{0}\) can be given locally near a regular point \(p\in\mathcal{C}_{0}\) as the graph of a function \(f(h(y),y,0)=0\), with \(h:\mathbb{R}^{n}\mapsto\mathbb{R}^{m}\). Thus, it is possible to reduce the slow subsystem \[\begin{align} 0&=f(x,y,0),\\ \dot{y}&=g(x,y,0), \end{align}\] to the simpler form \[\dot{y}=g(h(y),y,0).\] Hence, a subset \(\mathcal{S}_{0}\) is non–normally hyperbolic if the Jacobian \((\text{D}_{x}f)(p,0)\) has at least one eigenvalue with zero real part. According to the normally hyperbolic characteristic of the problem at hand, different analysis techniques can be used. For instance, non–normally hyperbolic points, which are related to dynamic features such as relaxation oscillations, that is, limit cycles of a singularly perturbed dynamical system [24], appearing in chemistry [42], geophysics [43], biology [44]–[46], or economics [47], [48], and canard trajectories [38], can be studied with the blow–up method [22], [49]–[51]. On the other hand, for normally hyperbolic cases a widely employed technique comes from the Fenichel’s theorem, stated as follows.

Once we have revisited the main concepts of fast–slow systems, in what follows we focus our interest in the unexpected canard trajectories, solutions that follow unstable branches of the critical manifold \(\mathcal{C}_{0}\) for considerable time, as the ones shown in the right panel of Figure 1, and which can arise from bifurcations of the fold type. Such trajectories have been observed in various applied sciences, for instance in neuroscience, allowing a detailed description of the very fast onset of large amplitude oscillations caused by small parameter variations in neural models [30], [31], [52], [53], as well as in chemistry [54]–[56], to mention a few.

A particularly important phenomenon near fold points of the critical manifold is the singular Hopf bifurcation [21], [22]. Consider a planar fast–slow system 38 with a non-degenerate fold point \(p \in C_0\) satisfying the standard conditions: \(\partial_x f(p,0) = 0\), \(\partial_{xx} f(p,0) \neq 0\), \(\partial_y f(p,0) \neq 0\), and \(g(p,0) \neq 0\). Suppose additionally that the system possesses an equilibrium point that, as a parameter (say \(\mu\)) varies, crosses the fold point \(p\). For \(0 < \varepsilon \ll 1\), Krupa and Szmolyan [21] showed that this crossing produces a Hopf bifurcation of the full system at a parameter value \(\mu = \mu_H(\varepsilon)\), with \(\mu_H \to \mu_0\) as \(\varepsilon \to 0\), where \(\mu_0\) corresponds to the equilibrium lying exactly at the fold. This is called a singular Hopf bifurcation because it degenerates in the singular limit.

The singular Hopf bifurcation is accompanied by a canard explosion: as \(\mu\) increases past \(\mu_H(\varepsilon)\), the small-amplitude limit cycle born at the Hopf point grows rapidly through a sequence of canard cycles until it reaches the size of a full relaxation oscillation [22]. This entire transition occurs within a parameter interval of exponentially small width \(O(e^{-c/\varepsilon})\), \(c > 0\). The canard cycles within this explosion are level sets of a Hamiltonian that arises in the blown-up analysis near the fold (see 5.2), a fact that is central to the controller design in Section 3.

The exponential sensitivity of canard cycles to parameter variations is the fundamental reason why these trajectories, despite being dynamically significant, are difficult to observe in uncontrolled systems subject to even minor disturbances. This motivates the development of feedback controllers that stabilise canard cycles as robust attractors, which is the main objective of the present work.

5.2 Planar folded canards↩︎

Let us start by recalling the canonical form of a canard point [21], given as \[\begin{align} \begin{aligned} \dot{x} &= -yh_{1}\left( x, y, \varepsilon, \alpha \right) + x^{2}h_{2}\left( x, y, \varepsilon, \alpha\right) + \varepsilon h_{3}\left( x, y, \varepsilon, \alpha \right),\\ \dot{y} &= \varepsilon \left( xh_{4}\left( x, y, \varepsilon, \alpha \right) - \alpha h_{5}\left( x, y, \varepsilon, \alpha \right) + y h_{6}\left( x, y, \varepsilon, \alpha \right) \right), \end{aligned} \label{Eq:CanardNormalForm} \end{align}\tag{42}\] with \((x, y)\in \mathbb{R}^{2}\), \(0<\varepsilon \ll 1\), and \(\alpha\) is a parameter. Furthermore, \[\begin{align} \begin{aligned} h_{3}(x, y, \varepsilon, \alpha) &= \mathcal{O}(x, y, \varepsilon, \alpha),\\ h_{i}(x, y, \varepsilon, \alpha) &= 1 + \mathcal{O}(x, y, \varepsilon, \alpha), \quad i = 1, 2, 4, 5, \end{aligned} \label{Eq:HFunctionsNormalForm} \end{align}\tag{43}\] and \(h_{6}\) is smooth. Let us simplify 42 along 43 as \[\begin{align} \begin{aligned} \dot{x} &= -y + x^{2} + \tilde{f}(x, y, \varepsilon, \alpha),\\ \dot{y} &= \varepsilon(x - \alpha + \tilde{g}(x, y, \varepsilon, \alpha)), \end{aligned} \label{Eq:SimplifiedNormalForm} \end{align}\tag{44}\] where the functions \(\tilde{f}\) and \(\tilde{g}\) collect the higher order terms in 42 . Thus, locally the critical manifold is denoted by \[\mathcal{C}_{0} = \left\{ (x, y) \in \mathbb{R}^{2}: -y + x^{2} + \tilde{f}(x, y, 0, \alpha) = 0 \right\}. \label{Eq:CriticalManifoldNormalForm}\tag{45}\] In particular, considering 44 in the absence of higher order terms yields \[\begin{align} \begin{aligned} \dot{x} &= -y + x^{2},\\ \dot{y} &= \varepsilon(x-\alpha). \end{aligned} \label{Eq:NormalFormNonHigherOrderTerms} \end{align}\tag{46}\] It is straightforward to check that the orbits of 46 , for \(\varepsilon>0\) and \(\alpha = 0\), are given by the level sets of the Hamiltonian \[H(x, y, \varepsilon) = \frac{1}{2}e^{-2y/\varepsilon}\left( \frac{y}{\varepsilon} - \frac{x^{2}}{\varepsilon} + \frac{1}{2} \right). \label{Eq:HamiltonianNormalForm}\tag{47}\] Furthermore, some orbits of 47 are canard trajectories, as the ones shown in Figure 14, and in fact it is well–known that canard cycles exist for \(H\in\left( 0, 1/4 \right)\) [21]. Therefore, the representation of the normal form for a folded canard is of high importance for our study as the control schemes designed are based on such description. However, it is important to highlight that one should be careful while designing such controllers as the dynamical structure of 42 could be altered, rendering inaccurate behaviours. For this reason, we present the following definitions.

Stabilising a canard is a more delicate objective than stabilising a hyperbolic equilibrium, and it requires a corresponding restriction on the feedback. The fold, its contact order, and hence the entire canard family are local objects, determined by the structure of the vector field in a neighbourhood of the fold point. A feedback that ignores this structure could meet a naive stabilisation target by altering the system rather than controlling it: it could, for instance, regularise the fold into a normally hyperbolic point, or change its contact order, in either case replacing the canard one wants to stabilise by an object of a different type. To exclude this, we require the open– and closed–loop vector fields to share the leading jet that defines the singularity at the fold (Definition 3). This is made precise as follows.

In other words, a compatible controller is free to reshape the higher–order terms, and thereby to select which level set \(\{H=h\}\), i.e.which canard cycle, becomes attracting, but it leaves the local singularity, the fold and its contact order, unchanged. The stabilised orbit is then genuinely a canard of the original system rather than an artefact engineered by the control.

Acknoeledgements↩︎

We sincerely thank the anonymous reviewers whose comments helped us to improve the manuscript.

References↩︎