Sliding modes of high codimension in piecewise-smooth dynamical systems

We consider piecewise-smooth dynamical systems, i.e., systems of ordinary differential equations switching between different sets of equations on distinct domains, separated by hyper-surfaces. As is well-known, when the solution approaches a discontinuity manifold, a classical solution may cease to exist. For this reason, starting with the pioneering work of Filippov, a concept of weak solution (also known as sliding mode) has been introduced and studied. Nowadays, the solution of piecewise-smooth dynamical systems in and close to discontinuity manifolds is well understood, if the manifold consists locally of a single discontinuity hyper-surface or of the intersection of two discontinuity hyper-surfaces. The present work presents partial results on the solution in and close to discontinuity manifolds of codimension 3 and higher.


Introduction
In this article, we consider piecewise-smooth dynamical systems, i.e., systems of ODEs which are characterized by a vector field that has jump discontinuities along certain hyper-surfaces of the Euclidean space. As general references, we address the reader to [2,13,15] and-from a more computational perspective-to [1].
There are several examples of this kind of systems in the scientific literature. E.g., many biological systems are represented in terms of switching functions (see, e.g., [20,Section 6]), which describe the effect of a certain variable on another, based on the fact that a certain threshold is reached and activates the coupling. Switching functions are frequently described by sigmoid (continuous) functions, but in many cases it is convenient to consider them as step functions (i.e., discontinuous ones). Gene regulatory networks (see, e.g., [16,19]), which are sets of genes (or parts of genes) that interact with each other to control a specific cell function, are also typically modeled by step functions that describe the switching interaction of this class of complex biological networks. Moreover, a strong interest in sliding mode control (see, e.g., [18]) has been shown in recent years, where the finite-time reaching of a certain discontinuity (sliding) manifold is the goal of a control system. Finally, when considering models based on state-dependent differential-algebraic delay equations, the manifolds are induced by the breaking points which characterize the solution [7]. The presence of discontinuities in all these models may give rise to termination or bifurcation phenomena. The lack of smoothness when the solution reaches one of the discontinuity hyper-surfaces leads to difficulties in dealing with its analysis, because the classical theory does not apply.
To circumvent this problem, the concept of weak solution (sliding mode) has been introduced [6,23]. Two methods have been proposed to analyze this situation: the regularization framework based on singular perturbation techniques (see, e.g., [22]) and the Filippov-like theory of differential inclusions. The theory developed by Filippov covers successfully the codimension-1 case, which is determined when a classical solution meets a single discontinuity hyper-surface. In such a case, it provides a concept of weak solution (a codimension-1 motion called sliding mode) in those cases where a classical solution ceases to exist. However, Filippov theory is not able to treat the higher codimension case, where the analysis is increasingly more difficult when the solution approaches the intersection of 2 or more discontinuity hyper-surfaces, which-even less common-can likely happen in the dynamical systems described above.
In some recent articles, the case has been considered, where a solution approaches the intersection of a pair of such discontinuity hyper-surfaces, with the possible occurrence of codimension-2 sliding modes. Due to the lack of uniqueness in classical Filippov theory, the analysis of the dynamics is non-trivial. It is based either on the analysis of regularizations of the piecewise-smooth dynamical system by a suitable singular perturbation or on the application of some specific selection principles. In [4] and [3], several possibilities have been proposed, based on some suitable criteria decided a priori, on how to select a Filippov sliding vector field on a codimension-2 discontinuity surface. Differently, a complete taxonomy has been obtained in [8], based on the analysis of a natural regularization which replaces step functions by continuous functions. A common feature of regularization techniques (e.g., [8,12,14]) is that of obtaining two time-scales for the dynamics, the faster of which describes the so-called hidden dynamics of the system and in many cases allows-passing to the limit with respect to the regularization parameter-to select a proper sliding mode. Based on the analysis of [8], a complete algorithm for the integration of piecewise-smooth dynamical systems with sliding modes of codimension at most 2 has been proposed in [10]. This algorithm is, up to our knowledge, the only one available to compute codimension-1 and 2 sliding modes in the numerical integration of a piecewise-smooth dynamical system. This is preferrable to the direct integration of the regularized system, because the stiffness of the singularly perturbed problem and the possible presence of high oscillations make the numerical integration very challenging.
This article is devoted to the case of higher dimensional intersections and related higher codimension sliding modes. Unfortunately, the analysis of the hidden dynamics is much more complicated and a complete classification is not possible. Nevertheless, we are able to prove some interesting results and to analyze to some detail special cases which are important in applications. E.g., situations where the vector fields have specific geometric properties (like the so-called nodal attractivity [3]) and where the systems are characterized by certain properties (like cooperativity [11]).
Organization of the article. In Section 2, we introduce the main concepts (sliding modes, regularization, and hidden dynamics). The main part of the paper is devoted to the study of the hidden dynamics in codimension 3 and higher. The maximal number of stationary points (which corresponds to sliding modes of maximal codimension) is studied in Section 3 and an example with six stationary points in codimension 3 is presented. This is obtained by applying Bezout's theory. Section 4 discusses the case of centrally symmetric vector fields and Section 5 the case of nodally attractive vector fields. An interesting class of examples shows that a vector field that is centrally symmetric and nodally attractive can have more than one stationary point. The study of the limit behaviour for τ → ∞ of the solution of the hidden dynamics is very challenging for d ≥ 3, because geometric arguments (e.g., Poincaré-Bendixson theorem) are not available or difficult to apply. A positive convergence result for nodally attractive vector fields with a hidden dynamics that is strictly cooperative is given in the final Section 6.

Main concepts
We closely follow the notation of [10].
We consider discontinuity hyper-surfaces where α : R n → R d (with d < n) is sufficiently differentiable and the hyper-surfaces intersect transversally. They divide the phase space into 2 d open regions where k = (k 1 , . . . , k d ) is a multi-index with k j ∈ {−1, 1}. The piecewise-smooth dynamical system is then given bẏ The smooth functions F k (y) are assumed to be defined in a neighbourhood of the closure of R k . In the discontinuity set d j =1 j , the right-hand side of (3) is thus multi-valued with values from the neighbouring domains.
We extend the concept of solution to the discontinuity set in the spirit of Filippov [5]. For an index vector k = (k 1 , . . . , k d ) with k j ∈ {−1, 0, 1}, we consider which is a discontinuity set, if at least one component k j vanishes. It reduces to (2) if all components are non-zero. For k = (k 1 , . . . , k d ), we define I k = {j | k j = 0}, and we let which collects the index vectors such that R touches R k . With this notation, we consider the differential-algebraic equation (DAE) with algebraic variables λ j , j ∈ I k . For λ j ∈ [−1, 1] the vector field in (5) is a convex combination of the vector fields F (y) (with ∈ N k ). The solution of (5) is therefore a Filippov solution of (3).
On the ε-neighbourhood of the discontinuity set (4), we consider the differential equatioṅ where the transition function π(u) is continuous, piecewise-smooth, and satisfies π(u) = −1 for u ≤ −1 and π(u) = 1 for u ≥ 1. For ease of presentation, we assume in the following that π(u) = u for |u| ≤ 1. This ordinary differential equation is on the set R ε a regularization of the discontinuous problem (3).

• Hidden dynamics
After entering an intersection of discontinuity hyper-surfaces (say at y * ∈ R k ), there are typically more than one Filippov solutions, and it is of interest to study which solution can be interpreted as the limit ε → 0 of the regularized differential equation. In the region R k ε , it follows from (7) that u i = α i (y)/ε satisfies which is a singularly perturbed problem. For its study, we introduce the fast time τ = t/ε, we denote the derivative with respect to τ by a prime, and we substitute the constant vector y * for y. This yields which is a regular dynamical system for u i , i ∈ I k . It is called hidden dynamics (a term coined in [12]). We note that this definition and all results of the present work are of local nature, because they depend on the fixed vector y * . The initial value for (9) is determined by the incoming solution at y * ∈ R k . Let k − be an index vector corresponding to a manifold of codimension one less, i.e., there is Assuming that the solution lies in R k − before entering R k at y * , then the initial value for (9) is given by u i * (0) = k i * , and the remaining u i (0) are such that the right-hand side of (9) vanishes for i = i * .
The hidden dynamics gives us much insight into the kind of solution leaving y * (cf. [8, Section 5.2]). We have: -sliding in R k , if the solution of (9) converges, for τ → ∞, to a stationary point with all components inside (−1, 1); . . , d} with the exception of those i ∈ I k for which u i (τ ) either tends to −∞ (then k + i = −1) or to +∞ (then k + i = 1). For the remaining i ∈ I k the solution u i (τ ) has a limit in If all expressions α i (y)F (y) are constant in a neighbourhood of y * , then the differential equations (9) are exact, and the above statements are obvious. The general case is less trivial, but we expect that the system (9) credibly describes the transient behaviour of the solution of the regularized differential equation. A rigorous treatment of the limit behaviour ε → 0 of the regularized solution is challenging, and even for the codimension 2 case it is still an open problem (see [9]).
Most of the present article deals with sliding modes of codimension 3. With the transition function π(u) = u for u ∈ [−1, 1] and for d = 3 the differential equations of the hidden dynamics (9) becomes (for i = 1, 2, 3 and u i ∈ [−1, 1]) We use the notation f k i = α i (y * )F k (y * ), and for the index vector k = (k 1 , k 2 , k 3 ) we only write the sign of k j . The hidden dynamics is completely described by the eight 3-dimensional vectors Sometimes, it is convenient to write (10) as with the multilinear polynomial The coefficients of this polynomial representation are given by In the following, we also consider the surfaces on which the flow (12) is orthogonal to the u i -axis. The set of stationary points is

Maximal number of sliding modes
In codimension 2 (intersection of two discontinuity hyper-surfaces), there are at most two stationary points in the hidden dynamics-one unstable and the other stable.
Here, we are interested in the maximal number of stationary points (i.e., sliding modes) in codimension three and higher. Without loss of generality we consider the transition function π(u) = u for |u| ≤ 1. We also assume that the index set in (9) is Theorem 1 Without any assumptions on the vector fields, the hidden dynamics (9) in codimension d can have at most d! isolated stationary points in R d .
Proof This bound is a consequence of Bézout's Theorem (see [17], where a proof based on basic results in algebraic geometry [21] is indicated). Introducing additional variables u 10 , . . . , u d0 , substituting in the polynomial system (9) the variables π(u j ) = u j by u j /u j 0 , and multiplying each equation by the product u 10 · . . . · u d0 yields homogeneous polynomials of degree 1 in the variables (u j , u j 0 ). Bézout's Theorem as formulated in [17] states that there are at most d! isolated solutions of the polynomial system. For convenience of the reader, we give an elementary proof in codimension d = 3 which is also useful for the construction of the stationary points. We consider the product u 1 u 2 as an independent variable, and we write the condition for a stationary point of (13) as ⎛ Solving this linear system by Cramer's rule shows that where p(u 3 ), q(u 3 ), r(u 3 ), d(u 3 ) are polynomials of degree 3. The necessary condition u 1 · u 2 = u 1 u 2 then yields the polynomial equation of degree 6, which has at most 6 solutions. The stationary points of (13) are then (u 1 , u 2 , u 3 ), where u 3 is a root of (16) and u 1 , u 2 are given by (15).
The following example shows that there exist problems having 6 stationary points in the unit cube.
With y * at the origin, we assume that the eight vector fields (in the ordering of (11)) are given by ⎛ By checking the eigenvalues of the Jacobian matrix of (12), we find that only the last stationary point in this list is asymptotically stable (one real and a complex pair of eigenvalues in the left complex half-plane), all others are unstable.
For an illustration of the dynamics of (12), we consider the surfaces i of (14), and we note that the set of stationary points is S = 1 ∩ 2 ∩ 3 . Figure 1 shows for 20 different values of u 3 cross sections of the surfaces 1 (red), 2 (blue), 3 (black) with horizontal planes (constant u 3 ). They appear as hyperbolas with horizontal and vertical asymptotes. The projection of the stationary points to these planes is plotted as small circles. If a stationary point lies on such a horizontal plane, it is highlighted by a larger circle. Therefore, we have included the values u 3 from (17) in the figure.
The sign pattern of the vector field at the corners of the square in the (u 1 , u 2 )-space with fixed u 3 is included. Note that for some value of u 3 close to −0.7 (between the 3rd and 4th pictures if Fig. 1), the hyperbola corresponding to 3 (black curve) degenerates to the union of the two lines (its two asymptotes).
We are curious to see the behaviour of the solution of the regularized differential equation, when it enters the intersection of all three discontinuity surfaces. We use ε = 10 −3 , but due to the relative scaling any other value would give the same pictures in Fig. 2. With the initial value y(0) = (1, 1, 1) , the solution approaches 2 at (0.1875, 0, 0.5) for t = 0.625. It remains O(ε)-close to it (which, in the limit ε → 0, is a codimension 1 sliding on 2 ) until t ≈ 0.733. Then, the solution approaches 1 ∩ 2 and remains O(ε)-close to it (i.e., a codimension 2 sliding in the limit) until it approaches the intersection 1 ∩ 2 ∩ 3 slightly after t = 3, as can be seen in the left picture of Fig. 2. From there on, the solution stays O(ε)-close to 1 ∩ 2 ∩ 3 (i.e., a codimension 3 sliding in the limit ε → 0).
The right pictures show the dynamics of the regularized problem close to the codimension 3 sliding. The scaled solution u i = y i /ε does not approach the stable stationary point (as could be expected), but it remains close to a periodic solution around the 4th stationary point of (17), which is unstable (one negative real eigenvalue and a complex pair of eigenvalues with positive real part).

Centrally symmetric vector fields
Sharper bounds on the maximal number of sliding modes (stationary points of the hidden dynamics) can be obtained, when the class of vector fields is restricted. A vector field is called centrally symmetric at y * ∈ R k , if (as in (10) we use the notation f i = α i (y * )F (y * )). This means that close to y * ∈ R k opposite vectors among {f | ∈ N k } differ only in its sign. For the equations Fig. 1 Cross sections with the horizontal planes u 3 = const for the problem with 6 stationary points, 1 (red), 2 (blue), 3 (black) of the hidden dynamics, this implies that a substitution ↔ − is equivalent to substituting u j ↔ −u j and multiplying the right-hand side by −1, so that only products with an odd number of factors can appear in a polynomial representation. In codimension 3, the hidden dynamics thus becomes (12) with  = (a 1 , a 2 , a 3 ) and c j = (c 1j , c 2j , c 3j ) given by the formulas after (13). For the statement of the following results, we consider the expression Theorem 2 Let D be the expression of (20).
• If D < 0, the origin is the only (real) stationary point of (19).
• If D > 0, we have three stationary points, one is at the origin, and the other two are The non-zero stationary points can be in or outside the unit cube. Proof An application of Cramer's rule gives the second and third relations of (21) and a formula for u 1 u 2 , which leads to u 2 3 = D.
For the study of stability of the dynamical system (12), we consider its Jacobian matrix

be a stationary point that is different from the origin. We then have
Proof The determinant of the rank-one perturbation (23) is given by the matrix determinant lemma as where adj(J 0 ) is the adjugate matrix of J 0 = J (0, 0, 0). At the non-zero stationary points, we have from (22) , c 1 , c 2 ) .
The formula for the adjugate matrix of J 0 yields J (0, 0, 0), which proves the statement of the theorem.
As a consequence, there cannot exist three asymptotically stable stationary points in R 3 .

Nodally attractive vector fields
A vector field is called nodally attractive at y * ∈ R k , if This means that all neighbouring vector fields point towards y * (recall that f i = α i (y * )F (y * )). It is natural to consider such vector fields, because every solution of (3) entering such a y * has to continue as a sliding in R k . (10) for the hidden dynamics in codimension 3. If the vector field is nodally attractive, the hidden dynamics generically has 1, 3, or 5 stationary points, of which 1, 2, or 3 are asymptotically stable, respectively.

Theorem 4 Consider
Proof We consider the intersection 1 ∩ 2 . Nodal attractivity implies that horizontal cross sections of 1 are parts of a hyperbola connecting the front and back faces of the cube, and those of 2 connect the left and right faces of the cube. Consequently, for every fixed Along the curve 1 ∩ 2 , the vector field points vertically upwards or downwards according to g 3 (u 1 , u 2 , u 3 ) > 0 or g 3 (u 1 , u 2 , u 3 ) < 0, respectively. It changes the direction when crossing the surface 3 , which happens at stationary points. 1 By nodal attractivity, we have g 3 (u 1 , u 2 , u 3 ) > 0 at the bottom face u 3 = −1, so that the lowest intersection with 3 is an asymptotically stable stationary point. Continuing along the curve 1 ∩ 2 , the stability at stationary points alternates between asymptotically stable and unstable. Since (again by nodal attractivity) the last intersection with 3 in the cube is asymptotically stable, we must have an odd number of stationary points in the cube.

Corollary 1 Assume nodal attractivity and central symmetry.
• If the origin is the only stationary point in the cube, then it is asymptotically stable. • If there are 3 stationary points in the cube, then the origin is unstable, and the other two stationary points are asymptotically stable.
Proof This is an immediate consequence of Theorems 2 and 4.
Intuitively, it is not obvious that nodally attractive vector fields can have more than one stationary point. The constructive examples of Section 5.1 show that there exist nodally attractive vector fields having 3 stationary points. We are not aware of the existence of nodally attractive vector fields with 5 stationary points in the unit cube.

Construction of nodally attractive vector fields with 3 equilibria
The first example of a nodally attractive vector field with 3 equilibria has been communicated in 2010 by Douglas Ulmer to L. Dieci and N. Guglielmi. Here, we present a systematic way of constructing such examples. We investigate the possibility of achieving nodal attractivity among the centrally symmetric vector fields with three equilibria. We let so that the expression D of (20) becomes D = α 3 β/α 1 α 2 . Assuming D > 0 we have 3 equilibria by Theorem 2, which are in addition to the origin. If β is small compared to the α j , the equilibria are inside the unit cube. This is achieved by choosing (c 1 , c 2 , c 3 ) close to a singular matrix. We arbitrarily put so that det(c 1 , c 2 , c 3 ) = 0 for δ = 0. Moreover, we put a j = a for all j , where a is another parameter. For this choice of coefficients, we have This yields D = −δ/a, so that the problem has 3 equilibria if a and δ have opposite sign. They are inside the unit cube if |δ| < |a|. Inspecting the values of g(u 1 , u 2 , u 3 ) (with components given by (13)) at the corners of the unit cube, one can see that the corresponding vector field is nodally attractive iff −2 + δ < a < −δ.
In summary, we see that the problem is centrally symmetric and nodally attractive with three equilibria in the unit cube iff 0 < δ < 1 and − 2 + δ < a < −δ.

Limit behaviour of the solution of the hidden dynamics
The behaviour of the solution u(τ ) of the hidden dynamics for τ → ∞, corresponding to initial values u(0) ∈ ∂C from the incoming solution, gives much insight into the discontinuous system (3). It is of high interest to study this limit. The system of Example 1 illustrates that in general it is very challenging or even impossible to predict the limit behaviour. We restrict the following analysis to simple (but important) situations in codimension 3.
When a codimension 2 sliding solution along i ∩ j enters the intersection i ∩ j ∩ k , the initial value for the hidden dynamics (12) lies on the intersection of i ∩ j with one of the faces u k = ±1, depending on the side from where the solution enters. We denote by E the set of all such initial values. In the case of a nodally attractive vector field, each face of the cube has exactly one element in E.

Nodal attractivity with cooperative hidden dynamics
We consider nodally attractive vector fields and, in addition, we assume that the system (12) of the hidden dynamics is strictly cooperative on the cube C, i.e., The dynamics of cooperative systems has been thoroughly studied by Hirsch [11]. It is shown that limit sets are invariant sets of systems in one dimension lower.
Theorem 5 Assume that the dynamical system (12) is strictly cooperative and corresponds to a nodally attractive vector field. Then, the solution corresponding to an initial value u 0 ∈ E converges to the stationary point of (12) that is closest to the face on which the initial value lies.
The proof of this theorem is postponed to Section 6.3, because it is closely related to the monotonicity of solutions that we shall study next.

Remark 1
The example of Section 5.1 is nodally attractive, if the parameters satisfy (25). It is strictly cooperative on the whole cube C, if |a| < 1. Hence, for 0 < δ < 1 and −1 < a < −δ, the solution of (12) corresponding to an initial value u 0 ∈ E converges to the stationary point that is closest to the face on which the initial value lies.

Monotone convergence
We let (u 10 , u 20 , u 30 ) ∈ E be an initial value of the hidden dynamics (12). Without loss of generality, we assume that We define the signs σ 1 , σ 2 , σ 3 ∈ {+1, −1} in such a way that so that the solution of (12) with initial value (u 10 , u 20 , u 30 ) enters the cube C in the set We further denote by D 0 the connected component of D for which (u 10 , u 20 , u 30 ) lies on the closure of D 0 . Recall that S denotes the set of stationary points in C.

Theorem 6
Assume that ∂D 0 has an empty intersection with the faces u 1 = σ 1 , u 2 = σ 2 , u 3 = σ 3 of the cube and that (for i = 1, 2, 3) Then, the solution of (12) with initial value (27) stays in D 0 and converges to a unique stationary point in the closure of D 0 .
Proof By definition of σ 1 , σ 2 , σ 3 the solution u 1 (τ ), u 2 (τ ), u 3 (τ ) enters the set D 0 . As long as it stays there, each component is a monotonic function (increasing for σ i = +1 and decreasing for σ i = −1). The border of D 0 is composed of parts of the faces u 1 = −σ 1 , u 2 = −σ 2 , u 3 = −σ 3 of the cube C, and of subsets of the surfaces 1 , 2 , 3 . We show that, with the exception of stationary points, the vector field (12) points on ∂D 0 into the interior of D 0 . At the faces of the cube this follows from the definition of D. At a point (u 1 , u 2 , u 3 ) ∈ (∂D 0 \ S) ∩ i , this is a consequence of the fact that the inner product of σ i g(u 1 , u 2 , u 3 ) with the gradient of g i is positive. Therefore, the solution cannot leave D 0 .
Moreover, the assumption that ∂D 0 has an empty intersection with the faces u 1 = σ 1 , u 2 = σ 2 , u 3 = σ 3 implies that the solution cannot leave the cube C. Due to the monotonicity of its components, we have convergence for τ → ∞ to a stationary point of (12) that lies on the border of D 0 .
We still have to prove that there is only one stationary point on ∂D 0 . Let u * be such a stationary point. We denote by B u * the subset of D 0 , which consists of those points that are mapped to u * by the flow of (12), i.e., ϕ τ (u) → u * for τ → ∞. We shall show that B u * is an open set. Take u ∈ B u * and let T > 0 be the time such that ϕ T (u) is sufficiently close to u * . By asymptotic stability of u * there is a neighbourhood U ⊂ D 0 of ϕ T (u) such that ϕ τ (w) → u * for τ → ∞ for all w ∈ U . The set ϕ −1 T (U ) ∩ D 0 is a neigbourhood of u that is entirely in B u * . Consequently, B u * is an open set.
If there were more than one stationary points in D 0 , the set D 0 could be written as the distinct union of open sets B u * , which is in contradiction with the connectedness of D 0 .

Proof of Theorem 5
We shall apply Theorem 6. Strict cooperativity of the system (12) implies the assumption on the gradient of g i , and nodal attractivity implies an empty intersection of ∂D 0 with the faces u 1 = σ 1 , u 2 = σ 2 , u 3 = σ 3 of the cube C (see a few lines below). An application of Theorem 6 thus yields convergence of the solution of (12) to the unique stationary point in D 0 . This stationary point is the one that is closest to the face, where the initial value lies.
It follows from the proof of Theorem 4 that i ∩ j is a curve that connects the faces u k = −1 and u k = 1 (k is the index different from i and j ) of the cube. By definition of D 0 the segment of the curve between the initial value u 0 ∈ E and the first stationary point is entirely in ∂D 0 . This completes the proof of Theorem 5.

Further remarks
The regularization (7) depends on the transition function π(u), and so does the differential equation (9) for the hidden dynamics which, in codimension 3 becomes u i = g i π(u 1 ), π(u 2 ), π(u 3 ) instead of (12). The results of the present article are restricted to the multilinear regularization (7), which implies that the vector field of the hidden dynamics is also multilinear.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4. 0/.