A dynamical system for prioritizing and coordinating motivations

TL;DR

This paper introduces a dynamical systems framework for robots to prioritize and switch between multiple recurring tasks, ensuring balanced task execution through evolving task values and bioinspired decision models.

Contribution

It presents a novel approach combining control theory and bioinspired decision making to enable robots to dynamically prioritize tasks and achieve stable task switching behavior.

Findings

Guarantees for recurrent task switching via stable limit cycles.

Large basin of attraction for the task coordination pattern.

Numerical evidence of robustness to perturbations.

Abstract

We develop a dynamical systems approach to prioritizing and selecting multiple recurring tasks with the aim of conferring a degree of deliberative goal selection to a mobile robot confronted with competing objectives. We take navigation as our prototypical task, and use reactive (i.e., vector field) planners derived from navigation functions to encode control policies that achieve each individual task. We associate a scalar "value" with each task representing its current urgency and let that quantity evolve in time as the robot evaluates the importance of its assigned task relative to competing tasks. The robot's motion control input is generated as a convex combination of the individual task vector fields. Their weights, in turn, evolve dynamically according to a decision model adapted from the literature on bioinspired swarm decision making, driven by the values. In this paper we study a simple case with two recurring, competing navigation tasks and derive conditions under which it can be guaranteed that the robot will repeatedly serve each in turn. Specifically, we provide conditions sufficient for the emergence of a stable limit cycle along which the robot repeatedly and alternately navigates to the two goal locations. Numerical study suggests that the basin of attraction is quite large so that significant perturbations are recovered with a reliable return to the desired task coordination pattern.

Figures

Equations

$$\underset{t\to +\infty }{\lim }x(t)=x_i^{\ast }.$$

$${\phi }_i(x)=\frac{1}{2}\Vert x-x_i^{\ast }{\Vert }_2^2.$$

$$\Phi (x)=\left[\begin{array}{cccc}{F}_1(x)& \dots & F_N(x)& 0\end{array}\right]\in {\mathbb{R}}^{d\times (N+1)}.$$

$$\dot{x}$$

$$\dot{m}$$

$$\dot{v}$$

$$\dot{x}=f_x(x,m)=-\Phi (x)m.$$

$${\dot{m}}_i={\tilde{v}}_i{m}_U-m_i\left(1/{\tilde{v}}_i-{\tilde{v}}_i{m}_U+{\sigma }_i(1-m_i-m_U)\right).$$

$${\dot{v}}_i$$

$$\dot{\xi }=f_{\xi }(\xi ),$$

$$\dot{x}$$

$${\dot{m}}_1$$

$${\dot{m}}_2$$

$${\dot{v}}_1$$

$${\dot{v}}_2$$

$$\Delta m=m_1-m_2,\bar{m}=\frac{m_1+m_2}{2};\Delta v=v_1-v_2,\bar{v}=\frac{v_1+v_2}{2}.$$

$$z=(x_1,x_2,\Delta m,\bar{m},\Delta v,\bar{v})$$

$$\dot{z}=f_z(z),$$

$${\dot{x}}_1=f_{x_1}(z)=$$

$${\dot{x}}_2=f_{x_2}(z)=$$

$$\dot{\Delta m}=f_{\Delta m}(z)=-$$

$$+\bar{v}\Delta m(1-2\bar{m})/{ϵ}_v+\Delta v(1-2\bar{m})(1+\bar{m})/{ϵ}_v,$$

$$\dot{\bar{m}}=f_{\bar{m}}(z)=\frac{1}{2}$$

$$+\frac{2\bar{v}+\Delta v}{2{ϵ}_v}(1-2\bar{m})(1+\frac{2\bar{m}+\Delta m}{2})$$

$$+\frac{2\bar{v}-\Delta v}{2{ϵ}_v}(1-2\bar{m})(1+\frac{2\bar{m}-\Delta m}{2})$$

$$\displaystyle-\frac{\sigma}{2}(2\bar{m}+\Delta m)(2\bar{m}-\Delta m)\biggr{)},$$

$${ϵ}_{\lambda }\dot{\Delta v}=f_{\Delta v}(z)=-(\Delta v-\Delta \phi (x)),$$

$${ϵ}_{\lambda }\dot{\bar{v}}=f_{\bar{v}}(z)=-(\bar{v}-\bar{\phi }(x)),$$

$$h(z_{r,1},z_{r,2},z_{r,3},z_{r,4}):=(z_{r,1},z_{r,2},z_{r,3},z_{r,4},z_{r,1},z_{r,2});h^{†}(z):=(z_1,z_2,z_3,z_4).$$

$${\dot{z}}_r=f_r(z_r,{ϵ}_v):=D{h}^{†}\cdot f_z\circ h(z_r),$$

$$-2(1+\sigma {ϵ}_v){\bar{m}}_d^2-(4{ϵ}_v^2+1){\bar{m}}_d+1=0,$$

$$\frac{-(4{ϵ}_v^2+1)+\sqrt{(4{ϵ}_v^2+1{)}^2+8(1+\sigma {ϵ}_v)}}{4(1+\sigma {ϵ}_v)}$$

$$=$$

$$z_{rd}:=h^{†}(z_d)=(x_{1rd},x_{2rd},\Delta m_{rd},{\bar{m}}_{rd}),$$

$$\dot{x}=f(t,x)$$

$$E=\{x;W(x)=0,x\in \bar{G}\}\text{ and }{E}_{\infty }=E\cup \{\infty \}.$$

Full text

{AMS}

37G25, 37D10, 37N35

1 Introduction

This paper develops a dynamical systems framework for planning and executing recurrent coverage tasks such as might be assigned a robot night watchman, sentry, or other patrol agent. Such problems have historically occasioned the introduction of temporal logic [26] or probabilistic [23] specification with subsequent hybrid implementation [20]. We seek to explore a more “embodied” approach, i.e., one that is capable of direct integration with a robot’s real-time sensorimotor models and controllers yet still offering a degree of deliberative judgment. A prototypical example of a real-time deliberative system is a foraging animal that achieves its basic needs for food and shelter by periodically revisiting different locations in its environment at different times apparently governed by some internal sense of relative urgency or satiety. In the vocabulary of psychology, the animal can be said to have drives which motivate it to perform actions that reduce those drives [37, 17, 30]. Inspired by the flexibility and robustness of such natural “reactively deliberate” systems, we seek a simple model of their drive-based decision-making mechanisms that might be robustly embodied within the dynamical sensorimotor layers of autonomous physical systems — a motivational dynamics for robots.

Dynamical systems approaches have been successful in understanding mechanisms for decision making in biological systems such as human choice behavior in two-alternative forced choice tasks [3], migration behavior in animal groups [24], and nest site selection behavior [35] in honeybee swarms. This literature has further motivated the exploration of such bio-inspired models for application to engineered systems [11]. Often, these decision mechanisms are value based in the sense that the organism can be interpreted as associating a numerical value with each available alternative and selecting the alternative with the highest value. Decision making in biological systems tends to be embodied in the sense that animals implement their decisions by moving their bodies in some way [25]. In the standard two-alternative forced choice task, an animal registers a decision by pushing a button or by looking at a particular point on a screen. In the context of migration or nest site selection, the animal moves its entire body to a new location. We take navigation, interpreted broadly as the task of steering a system to a desired goal state while avoiding obstacles, as the prototypical task for a mobile robot.

Vector field methods provide a natural way to encode the sensorimotor activity required to perform navigation tasks in dynamical systems language. When the vector field arises as the gradient of a well-chosen function, such as a navigation function [32], the system dynamics readily admit performance guarantees, such as proofs of convergence to the desired state while avoiding obstructions along the way. Furthermore, such vector field methods naturally map to control inputs for mechanical systems described by Lagrangian dynamics [19] and can be composed via linear combination or more intricate sequential [4] and parallel [6] operations.

These features motivate the consideration of vector field methods as an interesting alternative to the logical approaches to deliberation mentioned above, e.g., [20], whose hybrid (event-based) transitions require separately derived logical representation of the underlying dynamics. Instead, we seek an intrinsically dynamical systems approach to the composition and prioritization of potentially competing tasks that interprets the coefficients of their representative fields’ linear combinations as a kind of motivational state to be continuously adjusted in real time in a way that is flexible and robust to perturbations. In this initial work we seek to merely encode an activity composed of cyclicly-repeating base tasks, targeting, e.g., the recurrent patrol missions mentioned at the outset. The natural dynamical systems object by which to encode such an activity is a limit cycle. In subsequent work [31] we are studying empirical implementation of these ideas on physical robots and seek to enhance our limit cycle framework to allow the system to respond to external stimuli while keeping the limit cycle as the base behavior.

The main result of this paper is captured in Figure 1 which summarizes a numerical study illustrating two central analytical insights stated as Theorem 2.2 and Theorem 2.3. The motivational feedback path has a gain parametrized by $\epsilon_{v}>0$ and a time scale parametrized by $\epsilon_{\lambda}>0$. Numerical studies summarized by the four subsequent plots referenced by the numbered points of the figure indicate the presence of a stable limit cycle for a wide range of these parameter values. Analysis reveals that $\epsilon_{v}$ plays the role of an $\epsilon_{\lambda}$-dependent bifurcation parameter. Specifically, in the fast timescale limit $\epsilon_{\lambda}\to 0$, Theorem 2.2 establishes the existence of a Hopf bifurcation at a critical value of the feedback gain parameter $\epsilon^{*}_{v}(0)$. Further numerical study confirms the value of that formally-determined parameter, and suggests that the Hopf bifurcation persists along a curve of critical values, $\epsilon^{*}_{v}(\epsilon_{\lambda})$ for positive $\epsilon_{\lambda}$.

Seeking formal confirmation of the limit cycles suggested by those simulations at the physically interesting parameter values where $\epsilon_{\lambda}>0$, we next take recourse to a singular perturbation analysis. Specifically, we consider the joint limit $\epsilon_{v}\to 0,\epsilon_{\lambda}\to 0$ and carry out a dimension reduction of the system in this limit yielding planar dynamics exhibiting a limit cycle established by application of the Poincaré-Bendixson theorem. Arguments from geometric singular perturbation theory together with its conjectured (numerically corroborated) hyperbolicity then imply that this limit cycle persists for finite $\epsilon_{v},\epsilon_{\lambda}>0$.

This work is related to prior literature on dynamical decision-making in biological systems. Seeley et al. [35] studied nest site selection behavior in honeybee swarms and discovered a mechanism called a stop signal, by which bees who were committed to one nest site physically wrestled bees committed to other sites in order to get them to abandon their commitment. Seeley et al. constructed a dynamical systems model of this behavior and showed that the introduction of a stop signal allowed the system to avoid the deadlock state where no clear majority emerges in favor of any given option. We use the dynamical system from Seeley et al. [35] which models value-based nest site selection in honeybee swarms to modulate the motivation state. We let the value associated with each task be modulated by how far the agent is from the goal state associated with that task. This introduces feedback into the motivation dynamics by making the current system state influence the task values and thereby the motivation state.

Pais et al. [29] studied Seeley et al.’s model [35] using singular perturbation theory and showed that the stop signal also makes the model sensitive to the absolute value of the alternatives, allowing the system to remain in deadlock if all alternatives are equally poor. Pais et al. suggested that this sensitivity is useful to avoid prematurely committing to a suboptimal alternative, and show that it results in hysteresis as a function of the difference in the value of the alternatives. More recently, Gray et al. [11] proposed a model of decision-making dynamics for networked agents inspired by the above-referenced works on nest site selection and showed that singularity theory [9] can provide a powerful tool to engineer desired outcomes in decision-making models. These convincing accounts of the utility and potential analytical tractability of such bioinspired decision models provide a direct point of departure for our work. Specifically, [29, 35] studied one-off decisions where the value of each option (i.e., task) is static. In contrast, we allow the values of the tasks to change dynamically as they are completed by feeding back the system state, which allows the agent to determine the status of each task.

Other authors, particularly in the evolutionary dynamics literature, have studied systems with similar types of feedback. In evolutionary dynamics [16], which seeks to formalize Darwin’s ideas about natural selection, a set of populations each representing different strategies interact with each other and the interaction determines the level of fitness of each strategy. Fit populations thrive and grow, while unfit populations die off. Pais, Caicedo, and Leonard [28] studied the replicator-mutator equations from evolutionary dynamics with a particular network structure to the fitness function and showed conditions under which the dynamics exhibit Hopf bifurcations resulting in limit cycles. Mitchener and Nowak [27] studied evolutionary dynamics as a model of language transmission and showed conditions under which the dynamics of distinct grammars can exhibit limit cycles corresponding to periodic changes in the dominant grammar. The feedback model adopted by [28] and [27] captures the evolutionary process in which the fitness of a given strategy is determined by the relative fractions of the population adopting that strategy. Such a model is inappropriate for our robotic application, where the value of a task need not arise from competitive interactions between tasks. Our Hopf analysis in Section 4 is similar to that in [28], but we go on to show the existence of limit cycles in a two-dimensional region in parameter space using tools from geometric singular perturbation theory.

The contributions of this paper are twofold. First, we lay out an intrinsically dynamical systems approach to the composition and prioritization of potentially competing tasks for mobile robots. This approach provides an alternative to existing approaches using temporal logic and suggests future work investigating the connections between logic-based and dynamical-systems-based approaches to modeling decision making. Second, we show how a bistable system with a pitchfork bifurcation (the motivation dynamics) can be incorporated into a feedback system and produce a system that exhibits a Hopf bifurcation. This result is intuitive, and, in particular, the results of [9] suggest interpreting the oscillation our system exhibits as arising from the feedback modulation of the unfolding parameters of the pitchfork embedded in the motivation dynamics. Elucidating this structure, e.g., by finding the normal form of the Hopf bifurcation embedded in (9) is a nontrivial undertaking worthy of future study that may usefully inform further such designs.

The remainder of the paper is structured as follows. In Section 2 we lay out the broad class of systems under consideration before specifying the instance of the model which we study and stating our formal results. In Section 3 we show the result of several illustrative simulations, which suggest the existence of a Hopf bifurcation. In Section 4 we study the system in the limit $\epsilon_{\lambda}\to 0$ and show the existence of a Hopf bifurcation that results in stable limit cycles. In Section 5 we study the system in the joint limit $\epsilon_{\lambda}\to 0,\epsilon_{v}\to 0$ and show the existence of a stable limit cycle in the resulting two-dimensional reduced system; in Section 6 we show that this limit cycle persists for finite values of $\epsilon_{v}$ and $\epsilon_{\lambda}$. Finally, we conclude in Section 7.

2 Model, Problem Statement, and Formal Results

In this section we define our system model, state the problem we address and the formal results we obtain.

2.1 Model

Our model consists of three interconnected dynamical subsystems: states representing the navigation tasks and associated control actions (vector fields); the motivation state $m$; and the value state $v$. Implicit in the definition of the navigation tasks is the definition of the physical agent, which comprises the agent’s body and its workspace, or environment.

2.1.1 Body, Environment, and Motivational States

We model the robot as a point particle located at $x\in\mathcal{D}$, where the environment $\mathcal{D}\subseteq\mathbb{R}^{d}$ is a domain within Euclidean space. In general, $\mathcal{D}$ may be punctured by obstacles, but in this initial work we restrict ourselves to unobstructed domains.

We represent motivation by the state $m\in\Delta^{N}$, where $\Delta^{N}=\{m\in\mathbb{R}^{N+1}:m_{i}\geq 0,\sum_{i=1}^{N+1}m_{i}=1\}$ is the $N$-simplex. We index the first $N$ elements of $m$ by $i\in\{1,\ldots,N\}$: $m_{i}$ represents the motivation to perform task $i$. The last element we label as $m_{U}$: this represents undecided motivation, i.e., the decision to not perform any task.

2.1.2 Tasks

The agent has a set of $N$ tasks. Each task $i\in\{1,\ldots,N\}$ requires navigating the agent to the location $x_{i}^{*}\in\mathcal{D}$. Inspired by the navigation function framework [32], we assume the existence of a distance function $\varphi_{i}:\mathcal{D}\to\mathbb{R}_{+}$ for each task $i$. The function yields a gradient field $-\nabla\varphi_{i}$ such that $\dot{x}=-\nabla\varphi_{i}$ obeys

[TABLE]

That is, the gradient field $F_{i}=-\nabla\varphi_{i}$ is a vector field that accomplishes task $i$. In the following, where the domain is assumed to be unobstructed, we define the navigation functions by the squared Euclidean distance

[TABLE]

Finally, we define the matrix-valued function consisting of the $N$ task navigation vector fields plus the null gradient field associated with indecision

[TABLE]

By taking convex combinations of these vector fields we can assign the agent weighted combinations of the instantaneous (“greedy”) task plans they represent; the motivation state, defined below, will specify the convex combination to be taken at any given time.

The agent’s high-level mission is to repeatedly carry out each of the $N$ low-level tasks, i.e., visit each of the $N$ locations, in a specified order. In the vocabulary of the LTL hybrid systems literature, this corresponds to a recurrent patrol or coverage mission [7]. We now develop the detailed model, introducing its states and dynamics, then finally present statements of the problem we address and the formal results.

2.2 Model Dynamics

Having specified the system model and its state space, we now define its dynamics. The system has state $(x,m,v)\in\mathcal{D}\times\Delta^{N}\times\mathbb{R}_{+}^{N}$. The state variables evolve according to the dynamics

[TABLE]

The specific forms of the functions $f_{x},f_{m},$ and $f_{v}$ are given in the following paragraphs.

2.2.1 Navigation dynamics

The body’s location dynamics are the convex combination of the $N$ navigation vector fields (plus the null field associated with indecision), weighted by the motivation state:

[TABLE]

For example, when $m=[1,0,\cdots,0]^{T}$, the navigation dynamics are $\dot{x}=-\Phi(x)m=-\nabla\varphi_{1}(x)$, and when $m=[0,\cdots,0,1]^{T}$, the dynamics are $\dot{x}=0.$

2.2.2 Motivation dynamics

We take the motivation state dynamics from Pais et al.’s work [29] studying group decision making behavior in honeybee swarms:

[TABLE]

We set $\tilde{v}_{i}=v^{*}_{i}v_{i}$, where $v_{i}\in\mathbb{R}_{+}$ is the value of task $i$ and $v_{i}^{*}>0$ is a gain parameter that scales $v_{i}$. Equation (7) holds for each $i\in\{1,\ldots,N\}$, with the dynamics for $m_{U}$ following from the constraint that defines the simplex.

The dynamics (7) were derived for group decision making in [35] from a microscopic individual-level Markov process model that incorporates commitment, abandonment, recruitment, and stop signal mechanisms. The term $\tilde{v}_{i}m_{U}$ represents spontaneous commitment of an uncommitted individual to option $i$ at a rate which is proportional to the value $\tilde{v}_{i}$, $-m_{i}/\tilde{v}_{i}$ represents spontaneous abandonment, $\tilde{v}_{i}m_{i}m_{U}$ represents recruitment of an uncommitted individual by one committed to option $i$, and $-\sigma_{i}m_{i}(1-m_{i}-m_{U})$ represents a signal from individuals committed to options other than $i$ telling individuals committed to option $i$ to abandon their commitment. In our context where $m$ represents a single decision maker’s motivation state, each of these mechanisms can be interpreted as modeling specific processes between neurons in the decision maker’s brain rather than between individuals in a group.

2.2.3 Value dynamics

We define the dynamics of $v_{i}\in\mathbb{R}_{+}$, the value associated with state $i$, by

[TABLE]

Recall that the distance functions $\varphi_{i}$ take the value zero at the goal and increase with distance from the goal $x_{i}^{*}$. The intuition behind these dynamics is that the value (i.e., urgency) of a task should increase when the agent is far from that task’s goal state and decrease when the agent reaches that state; a larger $v_{i}$ means that task $i$ is more valuable, i.e., more urgent. We encode this intuition in the linear time-invariant dynamics (8), which cause $v_{i}$ to follow $\varphi_{i}(x)$ with a lag associated with the integration time scale $\lambda_{i}$. The dynamics correspond closely to the concept of drive reduction theory in social psychology, where motivation is thought to arise from the desire to carry out actions that satisfy various intrinsic drives [37, 17].

2.3 Formal Problem Statement and Analytical Results

The foregoing presentation introduces a broad class of models whose application to specific problems of reactive task planning and motivational control of multiple competing tasks we intend to explore empirically on physical robots. For the analytical purposes of this paper we find it expedient to consider a severely restricted instance from that class entailing only two, greatly simplified tasks and affording, in turn, a low-dimensional parametrization through imposition of various symmetries. In this section we first introduce the details of that restricted problem class and then state the analytical results we obtain.

2.3.1 Two Tasks, Their Essential Parameters, and New Coordinates

We have four parameters for each task $i\in\{1,\ldots,N\}$. Each task requires navigating to a goal location $x_{i}^{*}\in\mathcal{D}$. In the motivation dynamics (7), there is a positive stop signal parameter $\sigma_{i}>0$ and value gain $v_{i}^{*}>0$. Finally, in the value dynamics (8), there is a time scale $\lambda_{i}>0$. We show that the number of parameters can be greatly reduced and that the system’s behavior can be largely understood by varying the value of $v_{i}^{*}$.

For many parameter values, the system (6)–(8) exhibits a stable limit cycle in numerical simulations. To systematically study the system, we specialize to the case of a planar workspace $\mathcal{D}=\mathbb{R}^{2}$ and $N=2$ tasks. Then the state space of the system (6)–(8) is $\mathbb{R}^{2}\times\Delta^{2}\times\mathbb{R}_{+}^{2}$, for which we pick the coordinates $\xi=(x_{1},x_{2},m_{1},m_{2},v_{1},v_{2})$. Furthermore, we set the following parameter values.

We set the two goal locations to $x_{1}^{*}=(1,0)\in\mathcal{D},x_{2}^{*}=(-1,0)\in\mathcal{D}$. This choice is made without loss of generality, as it amounts to a translation, rotation, and scaling of the coordinates $x$ for $\mathcal{D}$. Scaling the coordinates $x$ by a factor $\alpha$ requires scaling the $v$ coordinates and parameters $v_{i}^{*}$ by a factor $\alpha^{2}$ to account for the fact that $\varphi_{i}(x)$ is homogenous of degree 2. For the stop signal $\sigma_{i}$, we follow Pais et al. [29] and impose the symmetry $\sigma_{1}=\sigma_{2}=\sigma$. Similarly, for ease of exposition and analysis we equate the value gain parameters $v_{1}^{*}=v_{2}^{*}=v^{*}>0$ as well as the value time scale parameters $\lambda_{1}=\lambda_{2}=\lambda>0$. With these choices the set of system parameters is reduced to $\sigma,v^{*},$ and $\lambda$, each of which must be positive. Fixing $\sigma$ at a nominal value, e.g., 8, leaves $v^{*}$ and $\lambda$ as free parameters whose values determine the behavior of the system.

The symmetry $\sigma_{1}=\sigma_{2}$ makes the decision-making mechanism, i.e., the motivation dynamics (7), obey an $S_{2}$ symmetry (i.e., symmetry under the interchange of the task numbers $1\leftrightarrow 2$) when the two task values are equal, i.e., $\tilde{v}_{1}=\tilde{v}_{2}$. The concrete interpretation of this symmetry is that the decision-making mechanism has no inherent bias for either underlying task. We find this characteristic to be desirable because it means that biases, if desired, can be introduced through the task values $v$. Isolating the biases in this way simplifies the analysis and eventual programming of the resulting system.111In the case where $\sigma_{1}=\sigma_{2}$ and $v^{*}_{1}=v^{*}_{2}$, the overall feedback system is $S_{2}$ invariant, though this may not be desirable in an application where the underlying navigation tasks are not of equal importance. For such an application, our formulation makes it natural to encode the difference in importance by breaking the symmetry $v_{1}^{*}=v_{2}^{*}$. The $S_{N}$ symmetry requirement naturally suggests further work in dynamical systems theory to study the dynamics of systems on $\Delta^{N}$ with the $S_{N}$ symmetry. Previous authors, e.g, [10], have studied dynamical systems using the symmetry approach, and some recent work, e.g., [11], has begun to use this theory in engineering applications; the present work suggests another application area for this theory.

Breaking the $S_{2}$ symmetry by allowing $\sigma_{1}$ and $\sigma_{2}$ to be set independently would introduce a bias towards one or the other options in the decision mechanism. Similarly, allowing $\lambda_{1}$ and $\lambda_{2}$ to be set independently also introduces a bias, this time in the valuation. Numerical experimentation suggests that the limit cycle behavior analyzed with the assumed symmetries persists for asymmetric settings of the parameters. Pais et al. [29] extensively studied the motivation dynamics in isolation and showed that, for fixed $v$, the dynamics (7) exhibit a pitchfork bifurcation as $\sigma=\sigma_{1}=\sigma_{2}$ is increased through a threshold value. As $\sigma$ is increased, the oscillations become increasingly nonlinear (as can be expected since $\sigma$ scales a nonlinearity in the dynamics (7)). In the limit of large $\sigma$ the dynamics tend towards a relaxation oscillation where the change in the motivation state $m$ happens quickly relative to the change in location state $x$. All evidence suggests that the parameters must be such that the motivation dynamics are in the post-bifurcation regime for the limit cycle to exist, however we will make careful analysis of the connection between the pitchfork bifurcation of the motivation dynamics and the Hopf bifurcation of the closed-loop system the subject of future work.

In the case $N=2$ and assuming $\varphi_{i}$ defined by (1), the equations (6)–(8) are

[TABLE]

where the components of $f_{\xi}$ are given by

[TABLE]

We change coordinates on $\Delta^{2}\times\mathbb{R}^{2}_{+}$ by transforming into mean and difference coordinates defined by

[TABLE]

Define the coordinates

[TABLE]

on the space $\mathbb{R}^{2}\times\Delta^{2}\times\mathbb{R}\times\mathbb{R}_{+}$ and parameters $\epsilon_{v}=1/v^{*}$ and $\epsilon_{\lambda}=1/\lambda$. It is easy to see that the transformation from $\zeta$ to $z$ is a diffeomorphism. In the mean-difference coordinates, the dynamics (9) are

[TABLE]

where the components of $f_{z}$ are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\Delta\varphi(x)=\varphi_{1}(x)-\varphi_{2}(x),\bar{\varphi}(x)=(\varphi_{1}(x)+\varphi_{2}(x))/2$.

Note that the system (9) exhibits a symmetric equilibrium state which we call a (symmetric) deadlock equilibrium, adopting the nomenclature established in the literature, e.g., [29]. The word deadlock refers to the fact that in such an equilibrium, the system fails to reach a decision. Precisely, we have the following:

Definition 2.1**.**

*Let $\xi_{d}$ be an equilibrium of the dynamics (9) of the form $\xi_{d}=(x_{1,d},x_{2,d},m_{d},m_{d},v_{d},v_{d})$. Then $\xi_{d}$ is called a (symmetric) deadlock equilibrium. Equivalently, in $z$ coordinates, an equilibrium $z_{d}$ of the dynamics (11) is called a (symmetric) deadlock equilibrium if it is of the form $z_{d}=(x_{1,d},x_{2,d},0,\bar{m}_{d},0,\bar{v}_{d})$. *

It is straightforward to see that a symmetric deadlock equilibrium exists if there is a point $x_{d}\in\mathcal{D}$ such that $\varphi_{1}(x_{d})=\varphi_{2}(x_{d})=\bar{\varphi}$.

2.3.2 Formal Results

In this section we state the two theorems which constitute the formal results of the paper. As can be seen in the simulations presented in Section 3, the dynamics (11) appear to exhibit a Hopf bifurcation as the parameters $\epsilon_{v}$ and $\epsilon_{\lambda}$ approach zero, giving birth to stable limit cycles. We formalize this observation in two steps. First we consider the limit $\epsilon_{\lambda}\to 0$ which reduces the dimension of the system (11) and permits an explicit computation showing the existence of a Hopf bifurcation.

In the limit $\epsilon_{\lambda}\to 0$, the $v$ variables are directly coupled to $\varphi(x)$, so $\Delta v=\Delta\varphi(x)$ and $\bar{v}=\bar{\varphi}(x)$, which are fixed points of Equations (16) and (17), respectively. Define $z_{r}=(x_{1},x_{2},\Delta m,\bar{m})$ as the vector of the remaining state variables. Explicitly, $z$ and $z_{r}$ are related by the linear embedding $z=h(z_{r})$ with left inverse given by the linear projection $h^{\dagger}$, where

[TABLE]

Then the dynamics (11) reduce to the restriction dynamics

[TABLE]

The restriction dynamics exhibit a Hopf bifurcation, as summarized in the following theorem:

Theorem 2.2**.**

*The system $\dot{z}_{r}=f_{r}(z_{r},\epsilon_{v})$ defined by (18) has a deadlock equilibrium $z_{rd}$ given by (21). For $\sigma>6$, the dynamics undergo a Hopf bifurcation resulting in stable periodic solutions at $(z_{rd},\epsilon_{v,0})$, where $\epsilon_{v,0}$ is the unique solution $\epsilon_{v}\in[0,1/2]$ of $4(2-\sigma)\epsilon_{v}^{2}-4(2-\sigma)\epsilon_{v}+(6-\sigma)=0$. For $\sigma=8,\epsilon_{v,0}=(3-\sqrt{6})/6\approx 0.09175$. *

As we are ultimately motivated by the physically meaningful case of small but non-zero values of $\epsilon_{v}$ and $\epsilon_{\lambda}$, we study the singular perturbation limit $\epsilon_{v},\epsilon_{\lambda}\to 0$ under which the system (11) can be reduced to a planar dynamical system and show the existence of a limit cycle. We then employ tools from geometric singular perturbation theory to show the persistence of this limit cycle for sufficiently small, but finite, values of $\epsilon_{v}$ and $\epsilon_{\lambda}$:

Theorem 2.3**.**

*Accepting Conjecture 6.2, below, for $\sigma>6$, there exists a stable limit cycle of (11) for sufficiently small, but finite, values of $\epsilon_{\lambda}$ and $\epsilon_{v}$. Equivalently, fixing $\lambda$, there exists a stable limit cycle of (11) for sufficiently large, but finite, values of $v^{*}$. *

3 Illustrative Simulations

Figure 1 summarizes the behavior of the system (11) as a function of the two parameters $\epsilon_{v}$ and $\epsilon_{\lambda}$. For large values of both parameters, the system exhibits a stable deadlock equilibrium, while for sufficiently small values of both parameters the system exhibits a stable limit cycle. As $\sigma$ is increased, the limit cycle tends towards a structure composed of a slow segment followed by a fast jump, which is characteristic of relaxation oscillations [2, 12]. Section 4 studies the system in the limit $\epsilon_{\lambda}\to 0$ and analytically shows the existence of a Hopf bifurcation at $\epsilon_{v}=\epsilon_{v,0}$ ($\approx 0.09175$ for $\sigma=8$). The blue line in Figure 1 shows the numerically-computed bifurcation value $\epsilon^{*}_{v}(\epsilon_{\lambda})$ for $\epsilon_{\lambda}>0$. The numerically-computed limit $\lim_{\epsilon_{\lambda}\to 0}\epsilon^{*}_{v}(\epsilon_{\lambda})$ corresponds well to the analytical value $\epsilon_{v,0}$.

Figures 2–5 show simulations of the system (11) for four representative values of the parameters $\epsilon_{v},\epsilon_{\lambda}$. We set $\sigma=8$. In $\xi$ coordinates, the initial conditions were $x=0,m_{1}=0,m_{2}=1/2,v_{1}=v_{2}=0.1$. In the mean-difference coordinates $z$ this corresponds to $\Delta m=-1/2,\bar{m}=1/4,\Delta v=0$, and $\bar{v}=0.1$. This choice of initial conditions was made to avoid the deadlock equilibrium but was otherwise generic.

Figure 2 suggests that for large values of $\epsilon_{v},\epsilon_{\lambda}$ there is a stable deadlock equilibrium in the system where no oscillations are present. Analyzing the dynamics, we see that this equilibrium corresponds to the state $z_{d}=(x_{1,d},x_{2,d},\Delta m_{d},\bar{m}_{d},\Delta v_{d},\bar{v}_{d})$, where $x_{1,d}=x_{2,d}=0,\Delta m_{d}=\Delta\varphi_{d}=0,\bar{v}_{d}=1/2$, and $\bar{m}_{d}$ solves the following quadratic equation:

[TABLE]

which has the solution $\bar{m}_{d}(\epsilon_{v})$ given by

[TABLE]

which is clearly positive, as $\sigma$ and $\epsilon_{v}$ are both positive, which implies that the second term under the radical in (20) is positive. Figures 2–5 suggest that the system undergoes a Hopf bifurcation as the parameters $\epsilon_{v}$ and $\epsilon_{\lambda}$ are decreased. In the following sections we carry out a series of analyses to characterize the bifurcation and study the resulting limit cycle.

4 Hopf analysis in the limit $\epsilon_{\lambda}\to 0$

Motivated by the numerical evidence of a Hopf bifurcation occurring at the deadlock equilibrium, we consider the system (11) in the limit $\epsilon_{\lambda}\to 0$ and analytically show the existence of a Hopf bifurcation in this limiting case as $\epsilon_{v}$ is lowered through a critical value $\epsilon_{v,0}$. We then numerically consider the case of finite $\epsilon_{\lambda}$ and compute the bifurcation value $\epsilon_{v}^{*}(\epsilon_{\lambda})$ for a range of values of $\epsilon_{\lambda}$; the numerically-computed limit $\lim_{\epsilon_{\lambda}\to 0}\epsilon_{v}^{*}(\epsilon_{\lambda})$ matches the analytical result $\epsilon_{v,0}$, as shown in Figure 1.

4.1 Dynamics in the limit $\epsilon_{\lambda}\to 0$

The limit dynamics (18) inherits the deadlock equilibrium $z_{rd}:=h^{\dagger}(z_{d})$ from the full dynamics (11), where

[TABLE]

$x_{1rd}=x_{2rd}=\Delta m_{rd}=0$, and $\bar{m}_{rd}$ again solves Equation (19).

In Figure 6 we show the numerically-computed bifurcation diagram for the system (18) with bifurcation parameter $\epsilon_{v}$. For large values of $\epsilon_{v}$, the deadlock equilibrium is stable. As $\epsilon_{v}$ is lowered below the critical value $\epsilon_{v,0}$, the system undergoes a Hopf bifurcation that results in a limit cycle. In Figure 6, we plot the amplitude of the oscillations of $\Delta\varphi$ for the limit cycle. As can be seen from Equation (10), $\Delta\varphi$ is constrained to take values in $[-1,1]$, so the limit cycle’s amplitude is bounded above by $1$.

4.2 Analysis of the $\epsilon_{\lambda}\to 0$ dynamics

Inspired by the bifurcation diagram, we now seek to show the existence of the Hopf bifurcation suggested by Figure 6. The following theorem from [13] summarizes the conditions under which a system undergoes Hopf bifurcation.

Theorem 4.1** (Hopf bifurcation, [13, Theorem 3.4.2]).**

Suppose that the system $\dot{z}=f(z,\mu),z\in\mathbb{R}^{n},\mu\in\mathbb{R}$, has an equilibrium $(z_{0},\mu_{0})$ and the following properties are satisfied:

The Jacobian $\left.D_{z}f\right|_{(z_{0},\mu_{0})}$ has a simple pair of pure imaginary eigenvalues $\lambda(\mu_{0})$ and $\bar{\lambda}(\mu_{0})$ and no other eigenvalues with zero real parts, 2. 2.

$\left.\mathrm{d}(\mathrm{Re}\ \lambda(\mu))/\mathrm{d}\mu\right|_{\mu=\mu_{0}}=d\neq 0.$**

Property 1) implies that there is a smooth curve of equilibria $(z(\mu),\mu)$ with $z(\mu_{0})=z_{0}$. The eigenvalues $\lambda(\mu),\bar{\lambda}(\mu)$ of $\left.D_{z}f\right|_{(z(\mu_{0}),\mu_{0})}$ which are imaginary at $\mu=\mu_{0}$ vary smoothly with $\mu$.

*If Property 2) is satisfied, then there is a unique three-dimensional center manifold passing through $(z_{0},\mu_{0})$ in $\mathbb{R}^{n}\times\mathbb{R}$ and a smooth system of coordinates (preserving the planes $\mu$=const.) for which the Taylor expansion of degree 3 on the center manifold is given by [13, (3.4.8)]. If $\left.\ell_{1}\right|_{(z_{0},\mu_{0})}\neq 0$, there is a surface of periodic solutions in the center manifold which has quadratic tangency with the eigenspace of $\lambda(\mu_{0}),\bar{\lambda}(\mu_{0})$ agreeing to second order with the parabaloid $\mu=-(\left.\ell_{1}\right|_{(z_{0},\mu_{0})}/d)(x^{2}+y^{2})$. If $\left.\ell_{1}\right|_{(z_{0},\mu_{0})}<0$, then these periodic solutions are stable limit cycles, while if $\left.\ell_{1}\right|_{(z_{0},\mu_{0})}>0$, the periodic solutions are repelling. *

The formulae for $\left.\ell_{1}\right|_{(z_{0},\mu_{0})}$, the first Lyapunov coefficient, are given in Appendix A.

Remark 4.2**.**

*In the statement of Theorem 4.1 we used $\mu$ as the bifurcation parameter for consistency with the notation of [13]. In the analysis in this paper, $\epsilon_{v}$ plays the role of bifurcation parameter. *

The Hopf bifurcation theorem applies to our system, as summarized in Theorem 2.2 stated in Section 2.3.2 and repeated below.

Theorem 2.2. The system $\dot{z}_{r}=f_{r}(z_{r},\epsilon_{v})$ defined by (18) has a deadlock equilibrium $z_{rd}$ given by (21). For $\sigma>6$, the dynamics undergo a Hopf bifurcation resulting in stable periodic solutions at $(z_{rd},\epsilon_{v,0})$, where $\epsilon_{v,0}$ is the unique solution $\epsilon_{v}\in[0,1/2]$ of $4(2-\sigma)\epsilon_{v}^{2}-4(2-\sigma)\epsilon_{v}+(6-\sigma)=0$. For $\sigma=8$, $\epsilon_{v,0}=(3-\sqrt{6})/6\approx 0.09175$.

Proof 4.3** (Proof of Theorem 2.2).**

Let $\epsilon_{v,0}$ be the solution $\epsilon_{v}\in[0,1/2]$ of $4(2-\sigma)\epsilon_{v}^{2}-4(2-\sigma)\epsilon_{v}+(6-\sigma)=0$, which exists for $\sigma>6$. By Lemma 4.4, the Jacobian $J_{0}$ of the system $\dot{z}_{r}=f(z_{r},\epsilon_{v})$ evaluated at the deadlock equilibrium $z_{rd}(\epsilon_{v})$ has a simple pair of pure imaginary eigenvalues when $\epsilon_{v}=\epsilon_{v,0}$. Therefore, the first condition of the Hopf bifurcation theorem is satisfied.

Lemma 4.6 establishes that $\left.\mathrm{d}(\mathrm{Re}\ \lambda(\epsilon_{v}))/\mathrm{d}\epsilon_{v}\right|_{\epsilon_{v,0}}\neq 0$ for the two pure imaginary eigenvalues $\lambda$ so the second condition of the Hopf bifurcation theorem is satisfied. The result then follows: the system (18) undergoes a Hopf bifurcation as the parameter $\epsilon_{v}$ is lowered through its critical value $\epsilon_{v,0}$.

*The first Lyapunov coefficient $\left.\ell_{1}\right|_{(z_{rd},\epsilon_{v,0})}$ is negative, as summarized by Lemma 4.8. This implies that the resulting limit cycles are stable. *

The following three lemmas, corresponding to the properties required by Theorem 4.1, contain the detailed arguments behind the proof of Theorem 2.2.

Lemma 4.4**.**

*Let $\sigma>6$ and let $J_{0}:=D_{z_{r}}f_{r}(z_{rd},\epsilon_{v})$ be the Jacobian of the system $\dot{z}_{r}=f_{r}(z_{r},\epsilon_{v})$ defined by (18) evaluated at the deadlock equilibrium $z_{rd}$ given by (21), considered as a function of $\epsilon_{v}$. Then, $J_{0}$ has a simple pair of two pure imaginary eigenvalues $\lambda(\epsilon_{v})$ and $\bar{\lambda}(\epsilon_{v})$ when $\epsilon_{v}=\epsilon_{v,0}$, where $\epsilon_{v,0}$ is the unique solution $\epsilon_{v}\in[0,1/2]$ of $4(2-\sigma)\epsilon_{v}^{2}-4(2-\sigma)\epsilon_{v}+(6-\sigma)=0$. *

Proof 4.5**.**

*See Appendix B.2. *

Lemma 4.6**.**

*Let $\lambda(\epsilon_{v})$ and $\bar{\lambda}(\epsilon_{v})$ be the simple pair of pure imaginary eigenvalues and let $\epsilon_{v,0}$ be as defined in Lemma 4.4. Then, $\left.\mathrm{d}(\mathrm{Re}\ \lambda(\epsilon_{v}))/\mathrm{d}\epsilon_{v}\right|_{\epsilon_{v,0}}<0$. *

Proof 4.7**.**

*See Appendix B.3. *

Lemma 4.8**.**

*Let $\ell_{1}=\left.\ell_{1}\right|_{(z_{0},\epsilon_{0})}$ be the first Lyapunov coefficient of the dynamics (18) evaluated at the deadlock equilibrium $z_{rd}$ given by (21). Then $\left.\ell_{1}\right|_{(z_{0},\epsilon_{0})}<0.$ *

Proof 4.9**.**

*See Appendix B.4. *

Theorem 2.2 then follows as a consequence of Lemmas 4.4, 4.6, and 4.8.

The implication of Theorem 2.2 is that the system (18) resulting from the limit $\epsilon_{\lambda}\to 0$ has a Hopf bifurcation at $\epsilon_{v}^{*}(0)=\epsilon_{v,0}$. However, as can be seen in Figure 4, limit cycle behavior persists for finite $\epsilon_{\lambda}$. One can numerically compute the eigenvalues of the linearization of the system (11) evaluated at the deadlock equilibrium $z_{d}$ and numerically show that a Hopf bifurcation occurs at a value $\epsilon_{v}^{*}(\epsilon_{\lambda})$. Figure 1 shows the numerically-computed values of $\epsilon_{v}^{*}(\epsilon_{\lambda})$ for a range of values of $\epsilon_{\lambda}$. It is clear that the numerical value for the limit $\lim_{\epsilon_{\lambda}\to 0}\epsilon_{v}^{*}(\epsilon_{\lambda})$ coincides with the analytical value $\epsilon_{v,0}$.

5 Reduction to a planar limit cycle in the joint limit $\epsilon_{\lambda},\epsilon_{v}\to 0$

The results from Section 4 strongly suggest the existence of a stable limit cycle for finite $\epsilon_{v},\epsilon_{\lambda}$. In this and the following section we make this conclusion rigorous by performing a series of reductions collapsing the dynamics (11) to a planar system in the joint limit $\epsilon_{\lambda}\to 0,\epsilon_{v}\to 0$. A Poincaré-Bendixson argument affords the conclusion that the planar system exhibits a stable limit cycle. Then, in the next section, we show that this limit cycle persists for small but finite $\epsilon$ by applying results from geometric singular perturbation theory.

5.1 A five dimensional attracting invariant submanifold

In our first reduction, formalized in Lemma 5.13, we note that $x_{2}$ must asymptotically converge to [math] independent of the other states’ behavior and that $x_{1}$ is attracted to the interval $[-1,1]$. Geometrically, this can be interpreted as the robot being attracted to the convex hull of its goal states $x_{1}^{*},x_{2}^{*}$. This observation reveals an attracting invariant submanifold of dimension five whose restriction dynamics we then study.

We begin by considering the dynamics of $x_{2}$ independently of the other five dynamical variables, which gives us a nonautonomous system $\dot{x}_{2}=f_{x_{2}}(t,x_{2}).$ The following results from [22] concern the asymptotic behavior of a nonautonomous system

[TABLE]

defined on $G\subseteq\mathbb{R}^{n}$. Let $G^{*}$ be an open set of $\mathbb{R}^{n}$ containing $\bar{G}$, the closure of $G$. We assume that $f:[0,\infty)\times G^{*}\rightarrow\mathbb{R}^{n}$ is a continuous (nonautonomous) vector field.

Definition 5.1**.**

Let $V:[0,\infty)\times G^{*}\rightarrow\mathbb{R}$ be a continuous, locally Lipschitz function. The function $V$ is said to be a Lyapunov function of (22) on $G$ if

1. i.

given $x$ in $\bar{G}$ there is a neighborhood $N$ of $x$ such that $V(t,x)$ is bounded from below for all $t\geq 0$ and all $x$ in $N\cap G$. 2. ii.

$\dot{V}(t,x)\leq-W(x)\leq 0$* for all $t\geq 0$ and all $x$ in $G$, where $W$ is continuous on $\bar{G}$. For $t$ where $V(t,x(t))$ is not differentiable, $\dot{V}$ is defined using the right-hand limit.*

If $V$ is a Lyapunov function for (22) on $G$, we define

[TABLE]

The statement of Theorem 5.4, below, requires the notion of absolute continuity, which is defined as follows.

Definition 5.2** ([34, Section 6.4]).**

A real-valued function $f$ on a closed, bounded interval $[a,b]$ is said to be absolutely continuous on $[a,b]$ provided for each $\epsilon>0$, there is a $\delta>0$ such that for every finite disjoint collection $\{(a_{k},b_{k})\}_{k=1}^{n}$ of open intervals in $(a,b)$,

[TABLE]

Lipschitz continuity implies absolute continuity, as follows:

Proposition 5.3** ([34, Section 6.4, Proposition 7]).**

*If the function $f$ is Lipschitz on a closed, bounded interval $[a,b]$, then it is absolutely continuous on $[a,b]$. *

Theorem 5.4** ([22, Theorem 1]).**

Let $V$ be a Lyapunov function for (22) on $G$, and let $x(t)$ be a solution of (22) that remains in $G$ for $t\geq t_{0}\geq 0$ with $[t_{0},\omega)$ the maximal future interval of definition of $x(t)$.

1. a.

If for each $p\in\bar{G}$ there is a neighborhood $N$ of $p$ such that $|f(t,x)|$ is bounded for all $t\geq 0$ and all $x$ in $N\cap G$, then either $x(t)\to\infty$ as $t\to\omega^{-}$, or $\omega=\infty$ and $x(t)\to E_{\infty}$ as $t\to\infty$. 2. b.

If $W(x(t))$ is absolutely continuous and its derivative is bounded from above (or from below) almost everywhere on $[t_{0},\omega)$ and if $\omega=\infty$, then $x(t)\to E_{\infty}$ as $t\to\infty$.

We now show that $x_{2}$ converges. For clarity of exposition, we write the argument as a series of lemmas.

Lemma 5.5**.**

Let $\varepsilon>0$ and $\mathcal{M}=\mathcal{D}\times\Delta^{2}\times\mathbb{R}\times\mathbb{R}_{+}$, and let the set $G$ be the subset of $\mathcal{M}$ defined by

[TABLE]

*The set $G$ is positive invariant under the dynamics $\dot{z}=f_{z}(z)$ defined by (11). *

Proof 5.6**.**

Let $z=(x_{1},x_{2},\Delta m,\bar{m},\Delta v,\bar{v})$ be coordinates for $\mathcal{M}=\mathcal{D}\times\Delta^{2}\times\mathbb{R}\times\mathbb{R}_{+}$ and consider the dynamics $\dot{z}=f_{z}(z)$ defined by (11).

Recall that $\varphi_{1}(x)=\|x-x_{1}^{*}\|^{2}/2,\varphi_{2}(x)=\|x-x_{2}^{*}\|^{2}/2$ with $x_{i}^{*}=[\pm 1,0]^{T}\in\mathcal{D}$. Direct calculation shows that $\bar{\varphi}(x)=(\varphi_{1}(x)+\varphi_{2}(x))/2=(1+\|x\|^{2})/2\geq 1/2$.

Furthermore, recall from (17) that $\dot{}\bar{v}=-\lambda(\bar{v}-\bar{\varphi})$, so $\dot{}\bar{v}(\bar{v}=1/2)=-\lambda(1/2-\bar{\varphi})\geq 0$ by the lower bound on $\bar{\varphi}$. Therefore the set $\{z\in\mathcal{M}|\bar{v}\geq 1/2\}$ is positive invariant.

*Similarly, note that $\bar{m}\geq 0$ by definition and that $-2\bar{m}\leq\Delta m\leq 2\bar{m}$, so $\bar{m}=0$ implies that $\Delta m=0$. Therefore, from (15), $\dot{}\bar{m}(\bar{m}=0)=\bar{v}/\epsilon_{v}$, so $\bar{v}\geq 1/2$ implies that $\dot{}\bar{m}(\bar{m}=0)\geq c/2\epsilon_{v}$ and therefore that $\bar{m}>0$. Therefore, the continuity of the $\dot{}\bar{m}$ dynamics implies that for $\bar{v}>1/2$, there exists an $\varepsilon>0$ such that $\bar{m}<\varepsilon$ implies that $\dot{}\bar{m}(\bar{m})>0$. This implies that $G$ is a positive-invariant set. *

We now study the dynamics of $x_{2}$, which can be written as

[TABLE]

Lemma 5.7**.**

*Let $G$ be the set defined by (23) in Lemma 5.5. Let $\varepsilon>0$ and $W(z)=2\varepsilon z_{2}^{2}=2\varepsilon x_{2}^{2}$. The function $V(z)=\frac{1}{2}z_{2}^{2}=\frac{1}{2}x_{2}^{2}$ is a Lyapunov function of (24) on $G$. *

Proof 5.8**.**

Note that $V(z),W(z)\geq 0$, so they satisfy requirement i of Definition 5.1. Computing $\dot{V}$ on $G$, we have

[TABLE]

*The inequality derives from the definition of $G$ and shows that $V$ satisfies the second requirement of Definition 5.1, thereby establishing the result. *

Lemma 5.9**.**

*The function $W(x_{2}(t))$ is absolutely continuous in $t$. *

Proof 5.10**.**

Proposition 5.3 shows that Lipschitz continuity implies absolute continuity. It is well known that a function is Lipschitz continuous if it has bounded first derivative, so we proceed by bounding the derivative $\mathrm{d}W(x_{2}(t))/\mathrm{d}t$.

The time derivative is

[TABLE]

This implies that

[TABLE]

*where $V$ is the Lyapunov function defined in Lemma 5.7. Then $|\mathrm{d}W/\mathrm{d}t|$ is clearly lower bounded by zero and upper bounded by $8\varepsilon V(x_{2}(0))$ since $V$ is a Lyapunov function for (24). *

Lemmas 5.5–5.9 imply that $x_{2}\to 0$. Furthermore, $x_{1}$ is attracted to the interval $[-1,1]$, as formalized in the following lemma.

Lemma 5.11**.**

*Let $x_{1}(t)=z_{1}(t)$ be the first component of a solution of (11) with initial condition $z(0)\in G$, where $G$ is as defined in (23). Then, $x_{1}(t)\to[-1,1]$ as $t\to\infty$. *

Proof 5.12**.**

Let $\varepsilon>0$ and note that $G$ is positive-invariant under (11). The dynamics of $x_{1}$ are $\dot{x}_{1}=\Delta m-2\bar{m}x_{1}$. From the definitions of $\Delta m$ and $\bar{m}$ as coordinates for the simplex $\Delta^{2}$ it is easily shown that $-2\bar{m}\leq\Delta m\leq 2\bar{m}$. Then, on $G$, $\dot{x}_{1}$ can be lower bounded as

[TABLE]

*which is strictly positive for $x_{1}<-1$. Similarly, $\dot{x}_{1}$ can be upper bounded such that $\dot{x}_{1}>0$ if $x_{1}>1$. Therefore, $x_{1}(t)\to[-1,1]$ as $t\to\infty$. *

Putting the preceding pieces together, we find that there is an attracting invariant submanifold of dimension five, as formalized in the following lemma.

Lemma 5.13**.**

*The variables $x_{1}$ and $x_{2}$ asymptotically converge to the interval $[-1,1]$ and the point 0, respectively. Therefore, the five-dimensional submanifold $\mathcal{M}_{5}$ defined by $\mathcal{M}_{5}=\{z=(x_{1},x_{2},\Delta m,\bar{m},\Delta v,\bar{v})\in\mathcal{M}|x_{1}\in[-1,1],x_{2}=0,\bar{m}>0,\bar{v}\geq 1/2\}$ is an attracting invariant submanifold under the (autonomous) dynamics (11). *

Proof 5.14**.**

Lemmas 5.7 and 5.9 show that $V$ is a Lyapunov function for the dynamical system $\dot{x}_{2}=f_{x_{2}}(x_{2})$ whose total derivative $\dot{V}$ is upper bounded by $-W(x_{2})$. Therefore the conditions of Theorem 5.4 hold.

*Applying Theorem 5.4, we find that $x_{2}$ converges to the unique root of $W(x_{2})$, which is located at $x_{2}=0$. Applying Lemma 5.11, we conclude that $x_{1}\to[-1,1]$. Therefore, $\mathcal{M}_{5}$ is an attracting subset of the set $G$ defined in (23). *

5.2 A two-dimensional slow manifold in the limit $\epsilon_{\lambda}\to 0,\epsilon_{v}\to 0$

We now proceed to eliminate $\Delta v$ and $\bar{m}$ by taking the singular perturbation limit $\epsilon_{\lambda}\to 0,\epsilon_{v}\to 0$. The first limit couples $\Delta v$ to $\Delta\varphi(x)$, as in Section 4, while the second limit pushes $\bar{m}$ to a slow manifold $\bar{m}(\Delta\varphi(x),\Delta m)$. This results in a planar “singular” system.

We study the dynamics (11) restricted to $\mathcal{M}_{5}$ defined in Lemma 5.13 by singular perturbation in the small parameter $\epsilon_{v}$. Further, let $\ell>0$ be a positive number such that $\epsilon_{\lambda}=1/\lambda=\ell\epsilon_{v}$. For a given value of $\ell$, this links $\epsilon_{\lambda}$ to $\epsilon_{v}$. Defining222This definition, which we make for consistency with the literature on singular perturbation problems, overloads the definition of $x\in\mathcal{D}$ from Section 2. $x=(x_{1},\Delta m)\in\mathcal{X}:=[-1,1]^{2}$ and $y=((1-2\bar{m})/\epsilon_{v},\Delta v,\bar{v})\in\mathbb{R}_{+}\times\mathbb{R}\times\mathbb{R}_{+}$, the restricted dynamics can be written in the standard form for a singular perturbation problem:

[TABLE]

where the slow dynamics is given by

[TABLE]

and the fast dynamics is given by

[TABLE]

It can be easily verified that the vector field $f_{x}$ points inward on the boundary of $\mathcal{X}=[-1,1]^{2}$, so $\mathcal{X}$ is positive invariant.

The slow manifold $\mathcal{M}_{s}$ is given by the implicit function solution $y=h_{y}(x)$ of the limiting fast dynamics $g(x,h_{y}(x),0)=0$:

[TABLE]

Figure 8 compares the analytically-computed slow manifold $h_{y,1}(x(t))$ to the value of $y_{1}(t)$ computed based on a numerical simulation of the full six-dimensional system (11). The small value of the error between the two values shows that the low-dimensional slow manifold is a good approximation to the trajectory of the high-dimensional system. Figure 9 shows the analytically-computed slow manifold surface $h_{y,1}(x)$ along with the trajectories $h_{y,1}(x(t))$ and $y_{1}(t)$.

The planar reduced dynamics on the slow manifold are given by the restriction of (26) to the slow manifold $\mathcal{M}_{s}$, now expressed in the coordinates of $\mathcal{X}=[-1,1]^{2}$,

[TABLE]

where the components are

[TABLE]

As seen in Appendix C, the slow manifold $\mathcal{M}_{s}$ is hyperbolic if the initial layer equation, (44), $\dot{\delta y}=D_{y}g_{y}(x,0,0)\delta y$ (where $\delta y=y-h(x)$) has a hyperbolic equilibrium at the origin. Away from the boundary of the set $\mathcal{X}=[-1,1]^{2}$, the eigenvalues of the linearization $\left.\frac{\partial g_{y}}{\partial y}\right|_{(x,0,0)}$ are equal to $-1/\ell$ (with multiplicity two) and $-(3-2x_{1}x_{2}+3x_{1}^{2})/2$. By definition $\ell>0$, so all three eigenvalues are strictly negative for $x$ in the interior of $\mathcal{X}$. Therefore they never intersect the imaginary axis and the slow manifold $\mathcal{M}_{s}$ is hyperbolic on the interior of $\mathcal{X}$. For points $x$ on the boundary of $\mathcal{X}$, the eigenvalues are not well defined, so we cannot conclude that $\mathcal{M}_{s}$ is hyperbolic at those points.

5.3 Existence of a stable limit cycle in the planar system

Now, we study the planar dynamics (30) and show by a straightforward application of the Poincaré-Bendixson Theorem [13, Theorem 1.8.1] that they exhibit a periodic orbit — a limit cycle — attracting an open annular neighborhood of the origin. We conjecture (and all numerical evidence corroborates) that this is an asymptotically stable limit cycle comprising the forward limit set of the entire origin-punctured state space. For present purposes it suffices to observe formally that an open neighborhood of initial conditions around the origin must take their forward limit set on this limit cycle.

Lemma 5.15**.**

*For $\sigma>6,$ there exists a periodic orbit of the reduced system (30) comprising the forward limit set of an open annular neighborhood excluding the (unstable) origin. *

Proof 5.16**.**

Note that the set $\mathcal{X}=[-1,1]^{2}$ is invariant set under the reduced dynamics $\dot{x}=f_{x}(x,h_{y}(x),0)$. Furthermore, note that the only equilibria in $\mathcal{X}$ are the origin and the two corners $(x_{1},x_{2})=(1,1),(-1,-1)$. It is easy to see that the interior of $\mathcal{X}$ is an invariant set; the edges $x_{2}=0$ and $x_{2}=1$ are part of the stable manifolds of the equilibria at the corners and the vector field is directed inwards on the other two edges. Let $\mathring{\mathcal{X}}$ be the interior of $\mathcal{X}$. Then $\mathring{\mathcal{X}}$ is pre-compact, connected, and contains a single fixed point at the origin.

The linearization of the reduced dynamics at the origin is given by

[TABLE]

The determinant $\det J_{r,0}=5\sigma/6>0$ for all $\sigma>0$ and the trace $\operatorname{tr}J_{r,0}=\sigma/6-1>0$ for all $\sigma>6$. The trace and determinant are, respectively, the sum and product of the two eigenvalues, so the fact that they are both positive implies that the eigenvalues are themselves positive. Therefore, for $\sigma>6$, the origin is an unstable focus.

*Since the entire annular region, $\mathring{\mathcal{X}}\setminus\{0\}$, is a pre-compact, positive invariant set possessing no fixed points it follows from the Poincaré-Bendixson Theorem that its forward limit set consists of proper (non-zero period) periodic orbits. In particular, in the neighborhood of the excluded repeller at the origin, there must exist an isolated periodic orbit comprising the forward limit set for that entire (punctured) neighborhood.333Although some texts, e.g., [36, Ch. 7.0], reserve the term “limit cycle” for isolated periodic orbits, the preponderant usage seems to favor that introduced in, e.g., [15, Ch. 10.6], [14, Ch. VII.1], [33, Ch. V.9], requiring merely that a periodic orbit comprise the limit set for some disjoint set of initial conditions, as we show here. Conjecture 6.2, below, will advance numerical evidence to support our assumption hereafter that this periodic orbit is hyperbolic, in which case it must be isolated [5, Ch. I.4.4], and hence qualifies unreservedly for the usage “limit cycle.” *

6 Persistence of the limit cycle for finite $\epsilon_{v},\epsilon_{\lambda}$

We now give conditions (which were previously stated as Theorem 2.3 in Section 2.3.2) under which the limit cycle whose existence was proven in the limit $\epsilon_{\lambda}\to 0,\epsilon_{v}\to 0$ in Lemma 5.15 persists for finite values of $\epsilon_{\lambda},\epsilon_{v}$. The result depends upon the conjectured hyperbolicity of that cycle, for which we establish numerical evidence below.

Theorem 2.3. * Accepting Conjecture 6.2, below, for $\sigma>6$, there exists a stable limit cycle of (11) for sufficiently small, but finite, values of $\epsilon_{\lambda}$ and $\epsilon_{v}$. Equivalently, fixing $\lambda$, there exists a stable limit cycle of (11) for sufficiently large, but finite, values of $v^{*}$.*

Lemma 5.15 shows that the singularly perturbed system with $\epsilon\to 0$, i.e., $\epsilon_{v},\epsilon_{\lambda}\to 0$ exhibits a limit cycle $\gamma_{0}$. The following result due to Fenichel [8] then allows us to show that this limit cycle persists for sufficiently small $\epsilon_{v},\epsilon_{\lambda}>0$.

Two pieces of notation are required to state the result. The symbol $\mathcal{E}_{H}$ represents the open set on which the linearization of the dynamics normal to the slow manifold has hyperbolic fixed points. In $\mathcal{E}_{H}$ the reduced vector field $X_{R}$ is defined by

[TABLE]

where $\pi^{\mathcal{E}}$ is the projection onto $\mathcal{E}$ defined in Equation (46) of Appendix C. We can now state the result.

Theorem 6.1** ([8, Theorem 13.1]).**

*Let $M$ be a $C^{r+1}$ manifold, $2\leq r\leq\infty$. Let $X^{\epsilon},\epsilon\in(-\epsilon_{0},\epsilon_{0})$ be a $C^{r}$ family of vector fields, and let $\mathcal{E}$ be a $C^{r}$ submanifold of $M$ consisting entirely of equilibrium points of $X^{0}$. Let $\gamma\in\mathcal{E}_{H}$ be a periodic orbit of the reduced vector field $X_{R}$, and suppose that $\gamma_{0}$, as a periodic orbit of $X_{R}$, has 1 as a Floquet multiplier of multiplicity precisely one. Then there exists $\epsilon_{1}>0$ and there exists a $C^{r-1}$ family of closed curves $\gamma_{\epsilon},\epsilon\in(-\epsilon_{1},\epsilon_{1})$, such that $\gamma_{0}=\gamma$ and $\gamma_{\epsilon}$ is a periodic orbit of $\epsilon^{-1}X^{\epsilon}$. The period of $\gamma_{\epsilon}$ is a $C^{r-1}$ function of $\epsilon$. *

Theorem 13.2 of [8] states that, when $\gamma_{0}$ is hyperbolic, the stability of the family $\gamma_{\epsilon}$ of periodic orbits for $\epsilon>0$ can be deduced from the stability of $\gamma_{0}$ and the stability of the linearization of $f_{x},g_{y}$ at $\epsilon=0$.

Conjecture 6.2**.**

*Let $\sigma=8$. The periodic orbit, $\gamma_{0}$, whose existence is guaranteed by Lemma 5.15 for the reduced dynamics $\dot{x}=f_{x}(x,h_{y}(x),0)$ defined by (30) is asymptotically stable with Floquet multipliers $\rho_{1}=1,\rho_{2}<1$. *

The Floquet multiplier associated with perturbations along the vector field is always equal to 1, and the remaining $n-1$ multipliers are the eigenvalues of the linearized Poincaré map $DP$ of the periodic orbit [13, Section 1.5]. We proceed by computing the numerical approximation to $P$ and $DP$ for the case $\sigma=8$ using the Poincaré section $\Sigma=\{(x_{1},0)|0.05<x_{1}<1\}$. The results are shown in Figure 10. From studying numerical solutions of the reduced dynamics, we know that the periodic orbit crosses $\Sigma$ at a point $p\approx 0.34934$, at which point the linearized map $DP(p)$ has a value of approximately 0.48. This implies that the two Floquet multipliers are $\rho_{1}=1,\rho_{2}\approx 0.48<1$.

In Conjecture 6.2, we only make claims about the Floquet multipliers in the case where $\sigma$ takes its nominal value $\sigma=8$. However, numerical investigation suggests that this result holds for generic $\sigma>6$. Theorem 2.3 then follows from applying Theorem 6.1 to the (now conjectured to be asymptotically stable) limit cycle found in Lemma 5.15.

Proof 6.3** (Proof of Theorem 2.3).**

Let $\gamma$ be a periodic orbit of the reduced system whose existence is shown in Lemma 5.15. By Conjecture 6.2, $\gamma$ has 1 as a Floquet multiplier of multiplicity precisely 1. On $\mathring{\mathcal{X}}$ the eigenvalues $\mu_{1},\mu_{2},\mu_{3}$ of the linearization $\partial g_{y}/\partial y$ can be computed in closed form and take the values

[TABLE]

On the interior of the slow manifold $(x_{1},x_{2})\in\mathring{\mathcal{X}}=(-1,1)^{2}$ these are bounded away from the imaginary axis, so $\gamma\in\mathcal{E}_{H}$. Then, by Theorem 6.1, there exists $\epsilon_{1}>0$ and a family of periodic orbits $\gamma_{\epsilon},\epsilon\in(-\epsilon_{1},\epsilon_{1})$ such that $\gamma_{0}=\gamma$.

*Specifically, for each $\epsilon\in(0,\epsilon_{1})$, there exists a stable limit cycle $\gamma_{\epsilon}$ with $\epsilon_{v}=\epsilon$ and $\epsilon_{\lambda}=\ell\epsilon$. Equivalently, fix $\lambda>0$ and define $v_{1}^{*}=1/\epsilon_{1}<+\infty$. Then for $v^{*}>v_{1}^{*}$, there exists a stable limit cycle $\gamma_{\epsilon}$ for $\epsilon=\epsilon_{v}=1/v^{*}$. *

Corollary 6.4**.**

*Theorem 2.3 establishes the existence of a stable periodic orbit $\gamma$ of the reduced dynamics for $(\epsilon_{\lambda},\epsilon_{v})\in\mathbb{R}_{+}\times[0,\epsilon_{1})$, i.e., the neighborhood of the $\epsilon_{\lambda}$ axis for sufficiently small $\epsilon_{v}>0$ and generic $\epsilon_{\lambda}$. *

Proof 6.5**.**

*Note that the fast dynamics (27) defines $\epsilon_{\lambda}=\ell\epsilon_{v}$, specifying only that $\ell>0$. Therefore, the result of Theorem 2.3 applies for parameter values $(\epsilon_{\lambda},\epsilon_{v})\in\{(\ell\epsilon_{v},\epsilon_{v})|\epsilon_{v}\in[0,\epsilon_{1}),\ell\in\mathbb{R}_{+}\}$. *

7 Conclusion

In summary, we have developed a dynamical systems method for managing motivations in physically-embodied agents, i.e., robots. This method provides a novel way for a system to autonomously and continuously switch between a set of vector fields, each of which defines a possible dynamics for the physical state corresponding to its performing a navigation task.

We specialize to the case where the system has two vector fields defined on a simply-connected subset of $\mathbb{R}^{2}$. By imposing several symmetries on our system, we are able to analyze the system in the limit where first one, and then both, of two parameters approaches zero. In the joint limit we reduce the system to a planar dynamical system by means of a singular perturbation analysis; a Poincaré-Bendixson argument shows that this planar system exhibits an isolated periodic orbit corresponding to the physical system state oscillating between two fixed points, one for each of the two vector fields. By appealing to geometric singular perturbation theory, we show that this periodic orbit persists for finite values of the two parameters.

A natural extension of this work is to consider cases where the system has more than two navigation tasks and where the domain $\mathcal{D}$ is punctured by obstacles, i.e., not simply connected. One natural way to extend this work to the case of multiple tasks is to decompose tasks into hierarchies encoded in binary trees; then, a variant of the system studied in this paper can run in each node to manage the motivations represented by each of its child nodes. By designing an appropriate method to feed the information from the child nodes back to their parent, it will be possible to maintain the limit cycle behavior for the larger number of tasks. Extending the analysis in this paper to the case of non-simply connected domains may prove more complex, as the analysis relies on several coordinate transformations that will be difficult to extend the more general case.

The other natural extension of this work is to apply it by implementing the motivational system on a physical robot. This implementation work is already in progress and will be the subject of a future report.

Appendix A First Lyapunov coefficient calculation

Kuznetsov [21, Section 5.4] provides the following formulae for computing $\left.\ell_{1}\right|_{(z_{0},\epsilon_{0})}$, the first Lyapunov coefficient of the dynamics $\dot{z}=f(z,\epsilon)$. Let $J_{0}=\left.D_{z}f\right|_{(z_{0},\epsilon_{0})}$. Property 1 of Theorem 4.1 implies that $J_{0}$ has two pure imaginary eigenvalues $\lambda(\epsilon_{0}),\bar{\lambda}(\epsilon_{0})=\pm i\omega_{0},\omega_{0}>0$. Let $q\in\mathbb{C}^{n}$ be a complex eigenvector corresponding to $\lambda(\epsilon_{0})$:

[TABLE]

Introduce the adjoint eigenvector $p\in\mathbb{C}^{n}$ satisfying

[TABLE]

and satisfying the normalization condition $\langle p,q\rangle=1$, where $\langle\cdot,\cdot\rangle$ is the standard inner product on $\mathbb{C}^{n}$.

Then, Taylor expand $f(z)=f(z,\epsilon_{0})$ to third order in $z$:

[TABLE]

where $B$ and $C$ are multilinear functions given by

[TABLE]

[TABLE]

which define the coefficients $B_{ijk}$ and $C_{ijkl}$. The coefficient $\left.\ell_{1}\right|_{(z_{0},\epsilon_{0})}$ is then given by Equation (5.62) of [21]:

[TABLE]

where $I$ is the identity matrix.

Appendix B Analysis of the Hopf bifurcation of (18)

In this section, we report computations relevant to the results in Theorem 2.2.

B.1 Jacobian computation

The following claim, which can be verified by direct computation, establishes the value of the Jacobian $J_{0}$ of the dynamics (18) evaluated at the deadlock equilibrium $z_{rd}$ given by (21).

Claim 1**.**

Let $J_{0}=\left.D_{z_{r}}f_{r}(z_{r},\epsilon_{v})\right|_{z_{r}=z_{rd}}$ be the Jacobian of (18) evaluated at the deadlock equilibrium defined by (21). Then

[TABLE]

where the non-zero components are given by

[TABLE]

The characteristic polynomial $p_{4}(\lambda)$ of $J_{0}$ is given by the determinant $|J_{0}-\lambda I|,$ where $I$ is the identity matrix. This determinant can be computed directly using expansion by minors:

[TABLE]

B.2 Proof of Lemma 4.4

Proof B.1** (Proof of Lemma 4.4).**

The characteristic polynomial of $J_{0}$ is computed in (35) in Appendix B.1 and can be expressed as

[TABLE]

where the final factor is given by coefficients arising directly from specific entries of the (sparse) Jacobian as

[TABLE]

and the components $j_{kl}$ of the Jacobian are as given in Equations (32)–(34) of Appendix B.1.

The roots of $p_{4}$ are given by

[TABLE]

It is clear that the first two roots are negative for all $\epsilon_{v}>0$, so the stability properties of the deadlock equilibrium are determined by the roots of $p_{2}$.

The roots $\lambda_{\pm}$ of $p_{2}$ are purely imaginary if $2\bar{m}_{rd}-j_{33}=0$. The condition $2\bar{m}_{rd}-j_{33}=0$ implies

[TABLE]

which is a valid value of $\bar{m}$ for $\epsilon_{v}\in[0,1/2]$. Inserting this expression into the expression (19) for $\bar{m}_{d}(\epsilon_{v})$, one finds that $2\bar{m}_{rd}-j_{33}=0$ implies that

[TABLE]

*For $\sigma>6$, this equation has two real-valued solutions, of which only one is in the relevant interval $[0,1/2]$ for which $\bar{m}_{rd}(\epsilon_{v})$ is defined. Therefore this solution is the relevant one defining the bifurcation value. *

B.3 Proof of Lemma 4.6

Proof B.2** (Proof of Lemma 4.6).**

Let $\epsilon_{v,0}$ and $p_{2}(\lambda)=\lambda^{2}+(2\bar{m}_{rd}-j_{33})\lambda-(j_{31}+2\bar{m}_{rd}j_{33})$ be defined as in the proof of Lemma 4.4, where $j_{kl}$ are defined in Appendix B.1. Let $\Delta_{p}=(2\bar{m}_{rd}-j_{33})^{2}+4(j_{31}+2\bar{m}_{rd}j_{33})$ be the discriminant of $p_{2}$ with respect to $\lambda$.

At the bifurcation value $\epsilon_{v,0}$, $2\bar{m}_{rd}-j_{33}=0$, so the discriminant $\Delta_{p}=4(j_{31}+4\bar{m}_{rd}^{2})$, which is negative, as can be shown by substituting in the expressions for $j_{kl}$ and grouping terms. This implies that the real part of the roots $\lambda_{\pm}$ are given by $(2\bar{m}_{rd}-j_{33})/2$. Therefore, since the coefficients of $p_{2}$ are continuous functions of $\epsilon_{v}$, the sign of the derivative $\mathrm{d}(\mathrm{Re}\ \lambda(\epsilon_{v}))/\mathrm{d}\epsilon_{v}|_{\epsilon_{v}=\epsilon_{v,0}}$ is determined by that of $(2\bar{m}^{\prime}_{rd}-j^{\prime}_{33})(\epsilon_{v,0}):=\mathrm{d}(2\bar{m}_{rd}-j_{33})/\mathrm{d}\epsilon_{v}|_{\epsilon_{v}=\epsilon_{v,0}}.$ Computing the derivative, we get

[TABLE]

which is non-zero (and in fact, positive) for all $\epsilon_{v}\in[0,1/2]$.

*Thus, the derivative $\mathrm{d}(\mathrm{Re}\ \lambda(\epsilon_{v}))/\mathrm{d}\epsilon_{v}|_{\epsilon_{v}=\epsilon_{v,0}}<0$. *

B.4 Criticality of the Hopf bifurcation in Theorem 2.2

The following result concerns the the Hopf bifurcation whose existence is proven in Theorem 2.2. It implies that the Hopf bifurcation is supercritical, so the periodic solutions created by the bifurcation are stable limit cycles.

The remainder of this section constitutes the proof of Lemma 4.8. As in Appendix A, write $\left.\ell_{1}\right|_{(z_{0},\epsilon_{0})}=\frac{1}{2\omega_{0}}\mathrm{Re}[T_{1}+T_{2}+T_{3}],$ where, from (31),

[TABLE]

$I$ is the identity matrix, and $\omega_{0},p,q,B,$ and $C$ are as defined in Appendix A. We show that $\mathrm{Re}\left[T_{i}\right]<0$ for $\sigma>6$. This implies that $\ell_{1}<0$, since $\omega_{0}>0$ by definition.

Let $J_{0}=D_{z}\left.f\right|_{(z_{0},\epsilon_{0})}$ be the Jacobian of the dynamics (18) evaluated at the bifurcation point. As shown in (36), $J_{0}$ has two purely imaginary eigenvalues when $2\bar{m}_{rd}-j_{33}=(2(1+2\epsilon_{v})\bar{m}_{rd}-(1-\epsilon_{v}^{2}))/(2\epsilon_{v})=0$. This implies that $1-2\bar{m}_{rd}=2\epsilon_{v}$ at the bifurcation point.

As shown in Claim 1, the Jacobian $J_{0}$ can be computed in closed form and takes the value

[TABLE]

where we have used the relationship established in (36) $1-2\bar{m}_{rd}=2\epsilon_{v}$ that holds at the bifurcation point in the final expression.

The eigenvalue problems $J_{0}q=i\omega_{0}q,\ J_{0}^{T}p=-i\omega_{0}p$ can be solved analytically, yielding

[TABLE]

The eigenvectors satisfy the required normalization condition $\langle p,q\rangle=1$. The sparsity of $p$ and $q$ greatly simplifies the computations of $T_{1},T_{2},T_{3}$.

B.4.1 Computing $T_{1}$

First, we compute $T_{1}$. Recall from (31) that $T_{1}=\langle p,C(q,q,\bar{q})\rangle$. Direct computation shows that $C_{i}(q,q,\bar{q})=0$ for $i\in\{1,2,4\}$. The relevant components of $C_{3jkl}$ are $C_{3111}=-96\epsilon_{v}\bar{m}_{rd},C_{3113}=C_{3131}=C_{3311}=-12\epsilon_{v}+(1-2\bar{m}_{rd})/\epsilon_{v}$. Substituting in the values of $p$ and $q$ then yields

[TABLE]

B.4.2 Computing $T_{2}$

Next, we compute $T_{2}$. Recall from (31) that $T_{2}=-2\langle p,B(q,J_{0}^{-1}B(q,\bar{q}))\rangle$. Direct computation shows that $B_{i}(q,\bar{q})=0$ for $i\in\{1,2,3\}$ and that $B_{4}(q,\bar{q})$ has relevant terms $B_{411}=(1-2\bar{m}_{rd})(1+\bar{m}_{rd})/\epsilon_{v}-12\epsilon_{v}\bar{m}_{rd},B_{413}=B_{431}=-2\epsilon_{v}-(1-2\bar{m}_{rd})/(2\epsilon_{v}),B_{433}=\sigma/2$. Then $B_{4}(q,\bar{q})=((3\sigma-9)+22\epsilon_{v}-4\sigma\epsilon_{v}^{2})/(4\omega_{0}^{2})$. The matrix $J_{0}^{-1}$ can be computed from (38) in closed form, and is equal to

[TABLE]

where $j_{kl}$ are defined in (38), so the only non-zero component of $J_{0}^{-1}B(q,\bar{q})$ is the fourth one, which is equal to $B_{4}(q,\bar{q})/j_{44}=((3\sigma-9)+22\epsilon_{v}-4\sigma\epsilon_{v}^{2})/(4\omega_{0}^{2}j_{44})$. Direct computation then shows that $B_{i}(q,J_{0}^{-1}B(q,\bar{q}))=0$ for $i\in\{2,4\}$ and then that

[TABLE]

B.4.3 Computing $T_{3}$

Finally, we compute $T_{3}$. Recall from (31) that $T_{3}=\langle p,B(\bar{q},(2i\omega_{0}I-J_{0})^{-1}B(q,q))\rangle$. As in the case of $T_{2}$, direct computation shows that $B_{i}(q,q)=0$ for $i\in\{1,2,3\}$ and that $B_{4}(q,q)$ has relevant terms $B_{411},B_{413}=B_{431}$, and $B_{433}$.

Let $\Gamma=(2i\omega_{0}I-J_{1})^{-1}$, which has the structure

[TABLE]

where $\gamma_{44}=1/(-j_{44}+2i\omega_{0})$. The first three components $(\Gamma B(q,q))_{i}=0$ for $i\in\{1,2,3\}$ and the only non-zero component is the fourth one given by $\gamma_{44}B_{4}(q,q)$. Then direct computation shows that $B_{i}(\bar{q},\Gamma B(q,q))=0$ for $i\in\{2,4\}$ so that $T_{3}=\bar{p}_{1}B_{1}(\bar{q},\Gamma B(q,q))+B_{3}(\bar{q},\Gamma B(q,q)).$ Straightforward but tedious calculations then show that

[TABLE]

B.5 Combining terms

Combining the results for $T_{1},T_{2},$ and $T_{3}$, we get that

[TABLE]

where $p_{5}(\epsilon_{v})=27(1+\sigma)+2(13-24\sigma+9\sigma^{2})\epsilon_{v}-8(40-9\sigma+8\sigma^{2})\epsilon_{v}^{2}+16(26-17\sigma+4\sigma^{2})\epsilon_{v}^{2}+304(\sigma-1)\epsilon_{v}^{4}-32(\sigma-1)^{2}\epsilon_{v}^{5}$. Let $d(\epsilon_{v})=3+2(\sigma-2)\epsilon_{v}-4(\sigma-1)\epsilon_{v}^{2})$. Then

[TABLE]

so the sign of $\ell_{1}$ is determined by the sign of the signs of $p_{5}(\epsilon_{v})$ and $d(\epsilon_{v})$ since $\epsilon_{v},\omega_{0}>0$. We proceed by showing that both polynomials are positive.

First, evaluating $d(\epsilon_{v})$ at $\epsilon_{v}=\epsilon_{v}(\sigma)$ defined by (37), we find that $d(\epsilon_{v}(\sigma))=(2\sigma(\sqrt{\sigma-2}-1)/(\sigma-2))$, which is clearly positive for $\sigma>3$. Similarly, $p_{5}(\epsilon_{v}(\sigma=6))=189>0$ and it can be shown that $\mathrm{d}p_{5}(\epsilon_{v}(\sigma))/\mathrm{d}\sigma>0$ for $\sigma>6$, so $p_{5}(\epsilon_{v}(\sigma))>0$ for $\sigma>6$. Therefore, $\ell_{1}|_{(x_{0},\mu_{0})}<0$ for $\sigma>6$, which establishes the desired result.

Appendix C A tutorial on geometric singular perturbation theory

Singular perturbation theory is a tool for studying dynamical systems characterized by two time scales, slow time $t$ and fast time $\tau$. The time scales are related by $\tau=t/\epsilon$, where $\epsilon>0$ is a small parameter. In the slow time scale, the dynamical system is governed by differential equations that are singular at $\epsilon=0$. By taking the limit $\epsilon\to 0$, i.e., assuming that the fast dynamics are much faster than the slow dynamics, one can often reduce a system to the slow dynamics.

Fenichel did fundamental work on this theory, for which [8] is a relatively comprehensive reference. Of particular interest to this paper is the theory he developed that allows one to relate the behavior of (the invariant manifolds of) a system in the limit $\epsilon\to 0$ to the behavior with finite $\epsilon>0$. In order to do this globally on a compact subset of the state space, Fenichel developed a geometric, or coordinate-free, notion of singular perturbation. The remainder of the section constitutes a summary of the relevant material in [8]. We begin by summarizing the local results, which are expressed in a given set of coordinates, before introducing the more abstract global, coordinate-free results.

C.1 Local results

Let $M$ be an open subset of $\mathbb{R}^{\mu}\times\mathbb{R}^{\nu}$, and let $\mathcal{E}=M\cap(\mathbb{R}^{\mu}\times\{0\})$ be nonempty. We consider a system of the form

[TABLE]

where $\dot{\ }$ denotes differentiation with respect to $t$, defined for $(x,y)\in M$, for small, real $\epsilon$. When $\epsilon=0$ the system (39) degenerates to the reduced system

[TABLE]

The second equation of (40) is an implicit function that defines $y$ as a function of $x$. The relation can be expressed explicitly, at least locally, as a function $y=h(x)$ [18]. The set $\{(x,y)|y=h(x)\}$ is called the slow manifold. By translating the $y$ coordinates by $-h(x)$, we can set $y=0$ on the slow manifold, which we denote by $\mathcal{E}$. Therefore, we assume that

[TABLE]

so that (40) defines a flow in $\mathcal{E}$, and we assume that this flow has a periodic orbit $\gamma_{0}:x=p(t),y=0.$ Fenichel’s aim is to describe the orbit structure of (39) for small nonzero $\epsilon$.

By rescaling time to $\tau=t/\epsilon$, we can transform (39) to

[TABLE]

where $\ {}^{\prime}$ denotes differentiation with respect to $\tau$. The set $\mathcal{E}$ consists entirely of equilibrium points of the system (42) in the limit $\epsilon\to 0$.

The plan is to relate the orbit structure of (39) near $\gamma_{0}$, for small nonzero $\epsilon$, to the orbit structure of the reduced system (40) near $\gamma_{0}$ and to the linearization of $\lim_{\epsilon\to 0}$ (42) at points of $\gamma_{0}$. The linearization of $\lim_{\epsilon\to 0}$ (42) at $(x,0)\in\mathcal{E}$ is

[TABLE]

where $D_{2}g(x,0,0)$ denotes differentiation with respect to the second argument of $g$ evaluated at $(x,0,0)$. The second component satisfies

[TABLE]

a linear equation parametrized by $(x,0)\in\mathcal{E}$. We refer to (44) (Equation (3.8) of [8]) as the initial layer equation.

The first qualitative question that Fenichel asks about (39) is whether it has a periodic orbit $\gamma_{\epsilon}$ near $\gamma_{0}$ for $\epsilon$ near zero. Fenichel [8] claimed that Anosov [1] obtained the most general result in the literature. In particular [8, Section III], Anosov proved that $\gamma_{0}$ can be continued to a family $\gamma_{\epsilon}$ if: (i) $\gamma_{0}$, regarded as a periodic orbit of the reduced system (40), has 1 as a Floquet multiplier of multiplicity precisely one, and (ii) for each $(x,0)\in\gamma_{0}$, the initial layer equation (44) has a hyperbolic equilibrium point at $\delta y=0$. The first condition can be interpreted as a nondegeneracy requirement on the periodic orbit $\gamma_{0}$ itself, while the second condition is sometimes called normal hyperbolicity of the slow manifold defined by $g(x,y,0)=0$. Theorem 13.1 of [8] makes this result precise.

C.2 Global results

The definitions up to here have been in a given set of coordinates. In order to properly account for limit cycles, Fenichel develops a global, coordinate-free notion of the singular perturbation problem. Let $M$ be a $C^{r+1}$ manifold, $1\leq r\leq\infty$. Let $X^{\epsilon}:M\to TM$ be a family of vector fields on $M$, parametrized by $\epsilon\in(-\epsilon_{0},\epsilon_{0})$, such that $X^{\epsilon}$ is a $C^{r}$ function of $(m,\epsilon)$. Let $\mathcal{E}$ be a $C^{r}$ submanifold of $M$ consisting entirely of equilibrium points of $X^{0}$, and let $z=\phi(m)$ be a $C^{r+1}$ local coordinate in $M$. In $z$-coordinates the flow of $X^{\epsilon}$ satisfies

[TABLE]

subject to the condition

[TABLE]

Let $\mu$ be the dimension of $\mathcal{E}$ and let $\nu$ be the codimension of $\mathcal{E}$ in $M$. Because $X^{0}$ vanishes identically on $\mathcal{E}$, $T_{m}\mathcal{E}$ is in the kernel of $TX^{0}(m)$ for any $m\in M$. In coordinates, then, $Tx^{0}(m)$ will have $\mu$ eigenvalues equal to zero and $\nu$ additional eigenvalues, which we call the nontrivial eigenvalues. The subspace $T_{m}\mathcal{E}$ is invariant under $TX^{0}(m)$, and so $TX^{0}(m)$ induces a linear map

[TABLE]

on the quotient space. The eigenvalues of $QX^{0}(m)$ are the nontrivial eigenvalues of the linearization of $\lim_{\epsilon\to 0}$ (45) at $z=\phi(m)$.

Let $\mathcal{E}_{R}$ be the open set where $QX^{0}$ is invertible. For each $m\in\mathcal{E}_{R}$, $T_{m}\mathcal{E}$ has a unique complement $N_{m}$ which is invariant under $TX^{0}(m)$. Let $\pi^{\mathcal{E}}$ denote the projection on $T\mathcal{E}$ defined by the splitting $TM|\mathcal{E}_{R}=T\mathcal{E}\oplus N.$ Let $\mathcal{E}_{H}\subset\mathcal{E}_{R}$ be the open subset where $QX^{0}$ has no pure imaginary eigenvalues; this is the normally-hyperbolic component of the slow manifold.

In $\mathcal{E}_{R}$ the reduced vector field $X_{R}$ is defined by

[TABLE]

Now we can state the main theorem that asserts conditions under which periodic orbits of the reduced vector field $X_{R}$ defined in the limit $\epsilon\to 0$ persist for $\epsilon>0$.

Theorem C.1** ([8, Theorem 13.1]).**

*Let $M$ be a $C^{r+1}$ manifold, $2\leq r\leq\infty$. Let $X^{\epsilon},\epsilon\in(-\epsilon_{0},\epsilon_{0})$ be a $C^{r}$ family of vector fields, and let $\mathcal{E}$ be a $C^{r}$ submanifold of $M$ consisting entirely of equilibrium points of $X^{0}$. Let $\gamma\in\mathcal{E}_{H}$ be a periodic orbit of the reduced vector field $X_{R}$, and suppose that $\gamma_{0}$, as a periodic orbit of $X_{R}$, has 1 as a Floquet multiplier of multiplicity precisely one. Then there exists $\epsilon_{1}>0$ and there exists a $C^{r-1}$ family of closed curves $\gamma_{\epsilon},\epsilon\in(-\epsilon_{1},\epsilon_{1})$, such that $\gamma_{0}=\gamma$ and $\gamma_{\epsilon}$ is a periodic orbit of $\epsilon^{-1}X^{\epsilon}$. The period of $\gamma_{\epsilon}$ is a $C^{r-1}$ function of $\epsilon$. *

For many applications, we are only interested in the case of small positive $\epsilon$. In [8, Section V], Fenichel explains how he is able to obtain results for $\epsilon\in(-\epsilon_{0},\epsilon_{0})$. Furthermore, the stability results of [8, Theorem 13.2] are stated for $\epsilon>0$. Let us now discuss how the theorem is applied. The main conditions are 1) that the periodic orbit $\gamma_{0}$ be contained in $\mathcal{E}_{H}$, the normally-hyperbolic component of the slow manifold, and 2) that $\gamma_{0}$ have 1 as a Floquet multiplier of multiplicity precisely one.

If one has a global coordinate system for $\mathcal{E}$, testing for normal hyperbolicity reduces to verifying that the eigenvalues of $\partial g/\partial y|_{\mathcal{E}}$ have non-zero real parts; if the real parts are negative, $\mathcal{E}$ is stable. The Floquet multipliers of $\gamma_{0}$ are the eigenvalues of $B$, the linearized Poincaré map of $\gamma_{0}$, so a multiplier of 1 corresponds to a fixed point of the Poincaré map, and multipliers less than (greater than) 1 correspond to stability (instability) of the orbit. There are $\mu$ Floquet multipliers $\rho_{i},i\in\{1,\ldots,\mu\}$, where $\mu$ is the dimension of $\mathcal{E}$. It can be shown [13] that $\det(B)=\prod_{i=1}^{\mu}\rho_{i}=\exp\int_{0}^{T}\operatorname{tr}(A(s))\mathrm{d}s$, where $T$ is the period of the periodic orbit and $A(s)$ is the linearization of the reduced dynamics $Df(x,0,0)|_{\gamma_{0}(s)}$. The existence of the periodic orbit means that there is one Floquet multiplier equal to 1. In general, $\rho_{i}$ have to be found by numerically computing $B$, unless one can bound the sign of $\operatorname{tr}A$ on the slow manifold. Alternatively, if one can show that the limit cycle is asymptotically stable on the slow manifold, the Floquet multiplier condition follows.

Acknowledgement

This work was supported in part by Air Force Research Laboratory grant FA865015D1845 (subcontract 669737-1).