Stochastically Perturbed Sliding Motion in Piecewise-Smooth Systems

Sliding motion is evolution on a switching manifold of a discontinuous, piecewise-smooth system of ordinary differential equations. In this paper we quantitatively study the effects of small-amplitude, additive, white Gaussian noise on stable sliding motion. For equations that are static in directions parallel to the switching manifold, the distance of orbits from the switching manifold approaches a quasi-steady-state density. From this density we calculate the mean and variance for the near sliding solution. Numerical results of a relay control system reveal that the noise may significantly affect the period and amplitude of periodic solutions with sliding segments.


Introduction
Nonsmoothness and noise are two features of a dynamical system that may be the cause of important qualitative behaviour. Hybrid and piecewise-smooth systems are utilized in a wide variety of fields to model phenomena that involve switching, impacts or other nonsmooth elements [1,2,3,4,5,6]. Recent studies have successively explained novel behaviour that may occur in such systems. For instance, so-called discontinuity maps have been developed and analyzed in order to understand border-collision scenarios in piecewise-smooth systems of ordinary differential equations [7,8,9]. Moreover, parameter uncertainty, background vibrations and other sources of noise are ubiquitous in real systems. Studies of stochastic differential equations have, for example, led to an understanding of noise-induced dynamics such as stochastic resonance and coherence resonance in excitable systems [10,11,12,13]. However, investigations into systems that are both piecewise-smooth and involve noise are relatively uncommon.
One-dimensional, piecewise-linear maps with noise have been the subject of some isolated investigations [14,15,16,17,18]. In [19], noise-induced transitions for a two-dimensional, piecewise-smooth system of ordinary differential equations are explained through an analysis of a one-dimensional, piecewise-linear return map with additive noise. Simple vibro-impacting systems have been analyzed with stochastic averaging [20]. Exact results are attainable for a classical, unforced, linear oscillator with elastic impacts [21,22], whereas more complex scenarios have been investigated asymptotically and numerically [23,24]. Stochasticity in switched control systems is particularly important in regards to robustness of an output signal to noise [25,26,27]. Here Lyapunov functions are an invaluable mathematical tool for determining stability because they do not necessitate smoothness in the vector field. Noise-induced oscillatory motion has been studied in piecewise-linear systems for which local linearity makes some key calculations tractable [28,29].
In this paper we consider piecewise-smooth, stochastic differential equations for which the underlying deterministic dynamics are described by the ODĖ where each Ω i ⊂ D ⊂ R N is open, nonempty and pairwise-disjoint, ∪ i Ω i = D, and each F i : Ω i → R N is a smooth function. Equation (1) is a Filippov system [30] and well-suited to model phenomena that alternate between different dynamical regimes, such as vibrating systems experiencing impacts or friction [31,32,33,34,35], and switching in electrical circuits [4,5,36]. Boundaries between the neighbouring subdomains, Ω i , are codimension-one surfaces termed switching manifolds. Often, a section of a switching manifold has the property that on either side of the manifold the vector field points towards the manifold. In this case any orbit that reaches the switching manifold becomes trapped on the manifold for some time. The resulting motion on the switching manifold is known as sliding motion, Fig. 1-A. Formally this is achieved by Filippov's method [37,30,1] which defines a vector field on the switching manifold by the unique convex combination of the two limiting vector fields on either side that is tangent to the switching manifold. Sliding motion corresponds to the sticking phase of stick-slip oscillators [38,39] and the coalesced regime of a piecewise-linear relay control system [40,41]. The addition of small noise pushes orbits off the switching manifold, but large excursions are curbed by the deterministic component of the system, Fig. 1-B. Thus the motion is balanced by the competing actions of noise and drift. In this paper we explore these dynamics more carefully.
We here summarize the basic effect of adding small noise to a smooth system. Below we compare this to our results for the piecewise-smooth system (1). Consider an N-dimensional where Φ is a smooth function. The addition of small amplitude, time-independent, white, Gaussian noise gives the stochastic differential equation where 0 < ε ≪ 1, B(x) is an N × N matrix with a smooth dependency on x, and W (t) is a standard Brownian motion [42,43]. Let p ε (x, t|x 0 ) denote the transitional probability density function (PDF) for the point x(t), given x(0) = x 0 . A straight-forward expansion in powers of √ ε reveals that the mean of p ε (x, t|x 0 ) differs from the deterministic solution (the solution to (2) with x(0) = x 0 ) by O(ε), whereas deviations are O( √ ε) [43,44].
A B Figure 1: Schematics of a Filippov system near a switching manifold that attracts orbits from both sides in the absence of noise, panel A, and with small amplitude additive noise, panel B.
In this paper we derive an analogous result for sliding motion for which the above method of expansion does not work because the vector field is discontinuous. Instead we analyze a one-dimensional, discontinuous stochastic differential equation for a quantity representing the distance from the switching manifold. We find that the mean solution differs from Filippov's deterministic sliding solution by O(ε), and deviations are O( √ ε), matching the smooth case. Moreover, the calculations suggest conditions necessary for the noise to induce a change in the dynamics that is not dominated by randomness.
Here we outline the remainder of the paper. We motivate our work in §2 by illustrating that noise may significantly reduce the amplitude and period of a solution to a prototypical relay control model involving segments of sliding motion. In §3 we introduce a system of twodimensional stochastic differential equations, (8)- (11), that describes stochastically perturbed sliding motion relating to a linear switching manifold in the case that system is the same in directions tangent to the switching manifold. In this case, the equation for motion in the direction orthogonal to the switching manifold, x, is independent of the variable representing displacement tangent to the switching manifold, y, and for this reason is amenable to an exact analysis. We leave a description of more general scenarios for subsequent work. In §4 we analyze the stochastic differential equation for x(t) and derive its quasi-steady-state distribution. In §5 we analyze the equation for y(t) and obtain expressions for the mean and variance of y(t). Derivations for this section are given in §6 and Appendix A. Conclusions are presented in §7.
are commonly modelled byẋ = Ax + Bu , where x ∈ R N , ϕ is the signal measurement and u is the control response, [40,5,1]. The system (4) is a Filippov system with a single switching manifold, {x | C T x = 0}, on which sliding may occur. In this system sliding corresponds to the idealized scenario of discrete switching events occurring continuously in time. Periodic orbits of (4) that involve sliding are described in [41,49,50,51].
As an example we consider the following canonical form, taken from [41] (also given in [1]), where ζ ∈ R is a parameter. Fig. 2 illustrates a stable periodic orbit of (4) with (5) involving sliding motion on the switching manifold. There are 12 separate sliding segments per period. These correspond to time intervals for which x 1 (the first component of the vector, x) is zero. The stability of periodic orbits with sliding may be determined by analyzing the Jacobian of a return map [41,51,52]. The robustness of periodic orbits with sliding has been briefly investigated by studying the size of the basin of attraction of the periodic orbit [52], and imposing a short time between consecutive switching events [53]. Fig. 3-A shows a typical orbit for the system when noise is added to the control signal as By comparing with Fig. 2-B, this shows that the noise may dampen the oscillations and induce an increase in frequency. Fig. 3-B shows that the average length of time per oscillation of x 3 decreases as the noise amplitude increases. We have observed similar behaviour for different values of ζ. In §7 we use the results below to speculate on the reason for this behaviour.

A simple set of equations for stochastically perturbed sliding motion
We are interested in the effects of noise on the dynamics of (1) near a switching manifold. In this section we introduce a simple set of equations that approximates (1) near a switching manifold, then formulate the inclusion of noise. If a switching manifold of (1) is locally C K we can choose a coordinate system such that in a neighbourhood of the origin, x = 0, the switching manifold is simply e T 1 x = O(K) [9]. In this paper we are not concerned with effects due to nonsmoothness in the switching manifold and for this reason suppose that the switching manifold is the coordinate plane, e T 1 x = 0. Two smooth subsystems govern the nearby flow, thus, locally, we may write the deterministic system aṡ where F (L) and F (R) are, say, C 1 . Suppose there exists a section of the switching manifold, call it Σ, for which e T 1 F (L) (x) > 0 and e T 1 F (R) (x) < 0, as in Fig. 1-A. In this scenario, forward orbits arrive at Σ from either side. We use Filippov's definition to define dynamics constrained to Σ [37,30,1]. Σ is known as an attracting sliding region and evolution on Σ is referred to as sliding motion.
However, with additive noise, the system (7) seems to be too complex for us analyze to a degree of detail that is useful. This is because both F (L) and F (R) depend on all components of the vector x, and we have been unable to analytically solve the resulting N-dimensional stochastic differential equation. Consequently, for this paper, which represents a first detailed analysis of sliding motion with noise, we ignore the dependency of F (L) and F (R) on components of x parallel to Σ as this enables us to reduce mathematical problems to one dimension but still capture what seems to be the essence of stochastically perturbed sliding motion. Moreover, this provides a useful approximation to the general case over short time-frames. To obtain each oscillation time, we computed an orbit up to t = 100 and identified the last three instances at which the value of x 3 changed sign (discounting rapid sign changes over a handful of grid points near this value due to noise), then subtracted the first time from the third time. Orbits were computed with the Euler-Maruyama method of fixed step size, ∆t = 0.0001.
Given that F (L) and F (R) are functions of only e T 1 x, the remaining N − 1 components of x may be treated identically and for this reason it suffices to study a two-dimensional system. We let x = [x y] T and add small amplitude, white, Gaussian noise independent to the state of the system. Assuming for simplicity that the noise in x is independent of the noise in y, the resulting stochastic differential equation may be written as where W 1 (t) and W 2 (t) are independent Brownian motions, 0 < ε ≪ 1 and κ > 0 are constants, and φ and ψ are piecewise-C 1 that for small |x| are given by We assume a L , a R > 0 , (11) to ensure that in the absence of noise the switching manifold (x = 0) is an attracting sliding region. Since φ and ψ are independent of ε, their coefficients are O(1). Consequently, for x(0) near zero, orbits of (8) likely remain near x = 0 for relatively long periods of time, as shown in §4.1, and for this reason we do not specify the behaviour of φ and ψ for large |x|.

Properties of x(t)
Since (8) lacks dependency on y, the equation for dx is decoupled from y: Given x(0) = x 0 , let p ε (x, t|x 0 ) denote the transitional PDF for the value of x(t), as governed by (12). Despite the discontinuity at x = 0, p ε (x, t|x 0 ) is unique and continuous on R × (0, ∞) × R. For x = 0 and t > 0, the PDF satisfies the Fokker-Planck equation with the initial condition p ε (x, 0|x 0 ) = δ(x − x 0 ), [42,43,54]. In §6 we provide an explicit expression for p ε (x, t|x 0 ) in the special case that φ is piecewise-constant. For the remainder of this section we use (13) to determine the long time behaviour of x(t).
If φ(x) > 0 for all x < 0, and φ(x) < 0 for all x > 0, then (12) has a steady-state density on R centered about the origin. Otherwise, φ(x) = 0 for some x = 0, and with nonzero probability orbits may cross this value of x and undergo dynamics far from the origin not described by the expansion (9). However, regardless of the global nature of φ, since ε is small the local attraction to the origin is relatively strong. Thus we expect orbits to remain near the origin for long periods of time and be distributed by a quasi-steady-state distribution for large but finite t. In §4.1 we determine the mean escape time of orbits from an O(1) neighbourhood of the origin. In §4.2 we use (13) to derive the quasi-steady-state probability density function asymptotically.

Escape from a neighbourhood of x = 0
For the function φ(x), (9), with (11), in the case that φ(x) = 0 for some x = 0, it is necessary to identify a value x b > 0, independent of ε, such that Then in the interval [−x b , x b ], the drift of (12) is towards x = 0. With ε > 0 and any initial condition Calculating the time to escape is a standard problem in the context of a single potential well, where the potential function is given by The mean escape time, T (x 0 ), may be found exactly [55,44,43]. Via Laplace's method of asymptotic evaluation of integrals [56], it follows that whenever |x 0 | < x b there exist ε-independent constants α 1 and α 2 such that

The quasi-steady-state probability density function
In §4.1 we showed that the mean escape time from an ε-independent neighbourhood [−x b , x b ] is exponentially large in 1 ε . Consequently, we can assume that the probability an orbit escapes represent the long time scale, and look for a solution to the Fokker-Planck equation (13) as a function of x andť. By substituting (18) with the WKB-type expansion [42,43], into (13), we arrive at Therefore and so by integrating φ(x) we obtain for an ε-independent function r. As a function of x and t, the dependence of this solution on t and x 0 appears only in O(ε M ) terms which may be ignored. Consequently we treat the solution as solely a function of x and refer to it as the quasi-steady-state solution, p qss,ε (x). Specifically (19) and (22) combine to give where we must have to ensure p qss,ε is normalized. For small ε and x 0 , the transitional PDF of (12), p ε (x, t|x 0 ), quickly settles to (23). The scalingx transforms (12) to from which we infer that p ε (x, t|x 0 ) approaches (23) on an O(ε) time-scale, when x 0 = O(ε). Furthermore, for times in the range ε 1−δ ≤ t ≤ ε −M , where δ > 0, it is reasonable to suppose x ∼ p qss,ε , in which case x sgn( which are useful in the next section.

Moments of y(t)
In this section we compute the mean of y(t) and conjecture the leading order term of its variance.
We assume x ∼ p qss,ε at all times under consideration which greatly simplifies calculations. We begin by deriving y(t) when ε = 0.

Deterministic sliding motion
When ε = 0, (8) is the Filippov system: As in [37,30,1], we define a vector field for sliding motion on the switching manifold (x = 0) by the unique convex combination of the two vector fields at the manifold that is tangent to the manifold. That is, ẋ for the unique scalar, q ∈ (0, 1), for whichẋ slide = 0. Solvingẋ slide = 0 gives q = a L a L +a R and thereforeẏ slide = Consequently, if (x(0), y(0)) = (0, y 0 ), then x slide (t) ≡ 0 and

The mean of y(t)
From (8) and (10) we have Integration yields and therefore If x ∼ p qss,ε , by substituting (27)- (29) we obtain where the ε-independent terms have combined to form y slide (t), (33). Therefore as ε → 0, the mean of y(t) limits on Filippov's sliding solution, y slide (t). This is non-trivial because Filippov's method, to obtain (33), and standard stochastic dynamical systems definitions, to obtain (37), are not immediately related. Note that the perturbation of y(t) from y slide (t) is order ε, mirroring the result for smooth systems, see §1. The explicit expression for the coefficient of the O(ε)-term in (37) is particularly useful. For instance, we can see that if c L = c R and d L = d R , then we require the asymmetry a L = a R in order for the O(ε)-term to be nonzero.

The variance of y(t)
The variance of y(t) may be computed via Var(y(t)) = y(t) 2 however this requires knowledge of p ε (x, t|x 0 ), for which we have not been able to obtain a useful expression in the case of general φ. We conjecture that the leading order terms of Var(y(t)) are independent of non-constant terms in φ and ψ because x(t) = O(ε) with high probability. Indeed this is consistent with numerical simulations, Fig. 4. In view of the result for the case that φ and ψ are piecewise-constant (Theorem 2, given below), we propose the following result: (8) with (11). Suppose x(0) is random with PDF, p qss,ε , and y(0) = y 0 . Then for any δ > 0, whenever From here until the concluding section, §7, we study the case that φ and ψ are piecewiseconstant, i.e.
6 Two-valued drift and a proof of Theorem 2 The stochastic differential equation (12) with (40): has been referred to as Brownian motion with two-valued drift. The transitional PDF of this process was first derived by Karatzas and Shreve in [57]. In this section we state this PDF and use it prove Theorem 2.

The transitional probability density function for x(t)
The transitional PDF for (43) is given by  Figure 4: The six horizontal bars are 95% confidence intervals for Var(y(t)) − κ 1, 1, 0), κ = 0, ε = 0.01, t = 1, and φ and ψ are piecewise-linear using different values of c L , c R , d L and d R as indicated by the table. Each confidence interval was obtained from 10 6 orbits computed numerically with the Euler-Maruyama method of fixed step size, ∆t = 0.0001. The results are consistent with Conjecture 1 which predicts that the values of c L , c R , d L and d R affect the magnitude of Var(y(t)) by at most O(ε 2 ). where is the PDF for the first passage time to zero of Brownian motion with constant drift, is the transitional PDF for Brownian motion with constant drift and an absorbing boundary condition at zero, and is the convolution relating to Laplace transforms. Fig. 5 shows (44) at different times.
In [57], Karatzas and Shreve derive (44) for ε = 1 by using Girsanov's theorem [58,59,60] and the trivariate PDF of Brownian motion, its positive occupation time, and its local time about zero. The result for ε = 1 follows simply from the scaling (25). The Laplace transform of (44) can be written as an integral-free expression and this was achieved in the earlier paper [61]. We do this below for x 0 = 0 and x > 0. In [62,63], the authors derive p ε (0, t|x 0 ) and use this to bound PDFs for a large class of scalar stochastic differential equations. In [64], p ε (x, t|0) is studied in the case a L , a R < 0. In [65], Zhang derived an expression for the transitional PDF of Brownian motion with a general bounded piecewise-continuous drift function. This could be used to analyze the PDF of (12) with general φ asymptotically, but such a calculation is beyond of scope of this paper.

Proof of Theorem 2
In the case of two-valued drift (40), the quasi-steady-state density (23) is a true steady-state defined for all x ∈ R: With the notation we have the following expression for the probability that x(t) > 0, given x 0 = 0.
Proof. For x > 0, in the case of two-valued drift (40), the Fokker-Planck equation (13) is Integration over for any ∆ > 0. Then integrating with respect to t and taking ∆ → 0 produces (50) where we also use lim which may be demonstrated by noting that as t → 0 + , p ε (x, t|0) is well-approximated by a zero-mean Gaussian. Proofs of the following two lemmas are given in Appendix A. Theorem 2 is an immediate consequence of Lemma 5 combined with (27) and (38). .
Lemma 5. Consider (43) and suppose a L , a R > 0 and x(0) is random with PDF, p ss,ε . Then for In the special case, a L = a R = a, we can write V ar(y(t)) exactly. In this case, by symmetry, Consequently the integral on the right-hand side of (78) vanishes, and (81) leads to which is consistent with (42).

Conclusions
When small noise is added to a Filippov system, orbits no longer slide along an attracting sliding section of a switching manifold. Instead, with high probability, orbits follow a random path near the switching manifold, Fig. 1-B. The average size of deviations from the switching manifold is governed by the strength of the noise relative to the magnitude of the vector field in a direction orthogonal to the switching manifold. The general N-dimensional stochastic differential equation formed by adding noise to the Filippov system (7) is particularly difficult to analyze due to the discontinuity and multiple dimensions. For this reason we made the supposition that the system is invariant along the switching manifold. This prevents an exploration of the effects of noise on orbits that reach the end of an attracting sliding region, but enables calculations to be reduced to one dimension. Moreover, with this reduction directions parallel to the switching manifold may be treated identically and hence it is sufficient to study the two-dimensional system (8).
For the system (8), x(t) denotes the displacement from the switching manifold and is governed by (12) with (9). Sample paths of (12) settle to the quasi-steady-state PDF, p qss,ε (23), on an O(ε) time scale. This PDF is not a true steady-state because orbits escape a neighbourhood of zero. However, since this occurs on an exponentially long time scale, see §4.1, it is suitable to assume x(t) is distributed by p qss,ε at times in a given range ε 1−δ ≤ t ≤ ε −M , for any δ, M > 0, where we may take small δ and large M, such that this is a long time interval. This assumption has the benefit of significantly simplifying our calculations.
Equation (34) is the stochastic differential equation for y(t), which represents displacement along the switching manifold. In the limit, ε → 0, the mean of y(t), (37), limits on Filippov's sliding solution, y slide (t). For ε = 0, the perturbation of y(t) from y slide (t) is O(ε), as for a generic smooth system. The perturbation depends on linear terms in φ and ψ. In order to gauge the effect of this perturbation on overall dynamics, it is necessary to compare it to the standard deviation of y(t). In the case that φ and ψ are piecewise-constant, Var(y(t)) is O(ε), see Theorem 2. For general φ and ψ, we conjectured that the leading order terms of Var(y(t)) are unchanged. Consequently deviations of y(t) from y(t) are O( √ ε). Therefore, assuming ε ≪ 1, we expect these deviations to dominate the difference between y(t) and y slide (t).
Although our above calculations are for differential equations that are static along the switching manifold, we can apply the basic principles gained to more general systems such as the relay control system, (4) with (5). First note that since the deterministic equations are independent of ε, when we make statements asymptotically in ε, we implicitly assume that the magnitude of any parameter in the deterministic equations is much less than, in particular, 1 √ ε . But, for (4) with (5), ∂ẋ 2 ∂x 1 = −ζ − 100, which is an extremely large value (we used ζ = 0.06). Since x 1 represents displacement from the switching manifold and x 2 is a direction parallel to the switching manifold, in the context of (8)-(10), ∂ẋ 2 ∂x 1 corresponds to the values of d L and d R . These values influence, to lowest order, the perturbation of the mean from the deterministic solution, (37), but not the deviation of sample paths from the mean, (39). This is a possible explanation for the noise-induced effect identified in §2. Specifically we found that when the parameters of (4) with (5) are tuned such that in the absence of noise orbits settle to oscillatory motion with sliding, the addition of noise may significantly decrease the average oscillation time, Fig. 3. A further analysis is required in order to make more definitive statements. We believe that the decrease in oscillation time with ε is a consequence of particular geometrical orientations and that an increase in oscillation time with ε is equally possible for this type of system. This paper leaves many avenues for future investigations. Perhaps foremost, we would like to understand the perturbation of y(t) from y slide in the case that the system explicitly depends on y. In particular we would like to understand how noise effects dynamics near the end of an attracting sliding region. Also it remains to consider more general forms for the noise, such as coloured noise or noise that is correlated in x and y, and study the effects of noise on other scenarios such as sliding bifurcations [66].
A.2 Proof of Lemma 5 Since (12) has no explicit time-dependence, Furthermore, To evaluate (67) using (68), we reorder the integrals of y, x and u, but to do this we must first transfer the u-dependence from the integrand to the limits of integration. We have where in the second line we conditioned over the first passage time, s, of (12) from x < 0 to 0. By Lemma 3, and since ∞ 0 h ε (s, x, −a L ) ds = 1, for all x < 0, we obtain A similar calculation with x > 0 produces By combining (67), (68), (70) and (71) we arrive at We now simplify the four terms of (72). We define so that the first term of (72) may be written as Similarly the third term is To the second term of (72) we reorder the integration such that dv is the outer-most integral instead of the inner-most integral. This step is straight-forward but requires some care with the limits of integration. We obtain and similarly for the last term of (72): We are now able to write (72) as t 0 t 0 sgn(x(s)) sgn(x(u)) ds du = t 2 − 2 Q ε (t, a L ) + Q ε (t, a R ) Via (44), it may be demonstrated that a R p ε (0, v|0) + ε 2 ∂p + ε ∂x (0, v|0) decays exponentially to zero as v → ∞ on an O(ε) time scale. For this reason it is helpful to apply the substitutionṽ = v ε to obtain t 0 t 0 sgn(x(s)) sgn(x(u)) ds du = t 2 − 2 Q ε (t, a L ) + Q ε (t, a R ) and expand Q ε (t − εṽ, a L ) − Q ε (t − εṽ, a R ) in ε such that the integral on the right-hand side of (79) has the form i j for some coefficients α ij . To obtain the coefficients, we first evaluate (73) via multiple applications of integration by parts: In view of (79), we use (81) to obtain Q ε (t, a L ) + Q ε (t, a R ) = t 2 2 − (a 3 L + a 3 R )εt 2a 2 L a 2 R (a L + a R ) Using Laplace's method [56] to asymptotically evaluate (80), the substitution of (82) and (83) into (79) yields t 0 t 0 sgn(x(s)) sgn(x(u)) ds du = (a 3 L + a 3 R )εt a 2 L a 2 R (a L + a R ) − 2(a L − a R )t 2 (a L + a R ) + 2(a 3 L − a 3 R )εt a 2 L a 2 R (a L + a R ) By applying Lemma 4 and simplifying we finally arrive at (56).