Representing and computing the B-derivative of an EC^r vector field's PC^r flow
RRepresenting and computing the B-derivative of an EC r vector field’s P C r flow ∗ George Council † , Shai Revzen ‡ , and Samuel A. Burden § Abstract.
This paper concerns first-order approximation of the piecewise-differentiable flow generated by a classof nonsmooth vector fields. Specifically, we represent and compute the Bouligand (or B-)derivativeof the piecewise- C r flow generated by an event-selected C r vector field. Our results are remark-ably efficient: although there are factorially many “pieces” of the desired derivative, we providean algorithm that evaluates its action on a given tangent vector using polynomial time and space,and verify the algorithm’s correctness by deriving a representation for the B-derivative that requires“only” exponential time and space to construct. We apply our methods in two classes of illustrativeexamples: piecewise-constant vector fields and mechanical systems subject to unilateral constraints. Key words. nonsmooth dynamical system, differential equation with discontinuous right-hand side, first-orderapproximation, Bouligand derivative, saltation matrix,
AMS subject classifications.
1. Introduction.
First-order approximations – i.e. derivatives – are a foundational toolfor analysis and synthesis in smooth dynamical and control systems. For instance, derivativesplay a crucial rˆole in: stability via spectral [17, Ch. 8.3] or Lyapunov [35, Ch. 5] methods;controllability via linearization [17, Ch. 8.7] or Frobenius/Chow [35, Ch. 8/Ch. 11] techniques;optimality via stationarity [4, Ch. 1] or Pontryagin [29, Ch. 1] principles; identifiability viaadaptation [34, Ch. 2] or Expectation-Maximization [23, Ch. 10] methods. These tools alldepend on the existence of a computationally-amenable representation for the first-order ap-proximation of smooth system dynamics – namely, the
Fr´echet (or F-)derivative of the system’ssmooth flow [28, Ch. 5.6], which derivative is a continuous linear function of tangent vectors. By definition, nonsmooth systems do not generally enjoy existence (let alone computa-tional amenability) of first-order approximations. Restricting to the class of (so-called [6,Def. 1, 2]) event-selected C r ( EC r ) vector fields that (i) are smooth except along a finite num-ber of surfaces of discontinuity and (ii) preclude sliding [20, 40] or branching [37, Def. 3.11]through a transversality condition, we obtain flows that are piecewise-differentiable [6, Thm. 4](specifically, piecewise - C r ( P C r ) [36, Ch. 4.1]). By virtue of their piecewise-differentiability,these flows admit a first-order approximation, termed the Bouligand (or B-)derivative, which ∗ Submitted to the editors DATE.
Funding:
This material is based upon work supported by the U. S. Army Research Laboratory and the U. S. ArmyResearch Office under contract/grant number W911NF-16-1-0158 and by the U. S. National Science FoundationCyber-Physical Systems Award † Department of Mechanical Engineering, Carnegie Mellon University, Pittsburgh, PA, USA ([email protected]). ‡ Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI, USA([email protected]). § Department of Electrical & Computer Engineering, University of Washington, Seattle, WA, USA ([email protected], http://faculty.uw.edu/sburden). We emphasize both properties of the Fr´echet derivative (continuity and linearity) since the generalizedderivative we consider in what follows will retain one (continuity) while relaxing the other (piecewise-linearity). This manuscript is for review purposes only. a r X i v : . [ m a t h . D S ] F e b G. COUNCIL, S. REVZEN, S. A. BURDEN derivative is a continuous piecewise -linear function of tangent vectors [36, Ch. 3, 4]. Thispaper is concerned with the efficient representation and computation of this piecewise-linearfirst-order approximation.Our contributions are twofold: (i) we construct a representation for the B-derivative ofthe
P C r flow generated by an EC r vector field; (ii) we derive an algorithm that evaluatesthe B-derivative on a given tangent vector. Although there are factorially many “pieces” ofthe derivative, we (i) represent it using exponential time and space and (ii) compute it usingpolynomial time and space. In an effort to make our results as accessible and useful as possible,we provide a concise summary of the algorithm in section 2 and apply our methods in section 3 before rehearsing the technical background in section 4 needed to derive the representationin section 5 and verify the algorithm’s correctness in section 6.We emphasize that our methods are most useful when there are more than two surfaces ofdiscontinuity, as representation and computation of first-order approximations in the 1- and2-surface cases have been investigated extensively [2, 5, 10, 11, 18, 19], and these cases do notbenefit from the complexity savings touted above. Previously, we established existence of thepiecewise-linear first-order approximation of the flow [6, Rem. 1] and provided an inefficientscheme to evaluate each of its “pieces” [6, Sec. 7] in the presence of an arbitrary numberof surfaces of discontinuity. To the best of our knowledge, the present paper contains thefirst representation for the B-derivative of the P C r flow of a general EC r vector field andpolynomial-time algorithm to compute it.
2. Algorithm.
The goal of this paper is to obtain an algorithm that efficiently computesthe derivative of a class of nonsmooth flows. This computational task and our solution are easyto describe, yet verifying the algorithm’s correctness requires significant technical overhead.Thus, the remainder of this section will be devoted to specifying the algorithm and the problemit solves using minimal notation and terminology. Subsequent sections will provide technicaldetails – which may be of interest in their own right – that prove the algorithm is correct.Given vector field F : R d → T R d and trajectory x : [0 , ∞ ) → R d satisfying (2.1) ∀ t ≥ x t = (cid:90) t F ( x τ ) dτ, our goal is to approximate how x t varies with respect to x to first order for a given t > φ : [0 , ∞ ) × R d → R d denoting the flow of F satisfying(2.2) ∀ t ≥ , x ∈ R d : φ t ( x ) = (cid:90) t F ( φ τ ( x )) dτ, our goal is to evaluate the directional derivative Dφ t ( x ; δx ) given t > δx ∈ T x R d :(2.3) ∀ t > , δx ∈ T x R d : Dφ t ( x ; δx ) = lim α → + α ( φ t ( x + α δx ) − φ t ( x )) . Specifically, we seek to evaluate this derivative for vector fields that are smooth everywhereexcept a finite collection of surfaces where they are allowed to be discontinuous. We will In this section, we will denote time dependence using subscripts rather than parentheses.
This manuscript is for review purposes only. -DERIVATIVE OF AN EC r VECTOR FIELD’S
P C r FLOW 3 first recall how to obtain the derivative in the presence of zero (subsection 2.1) and one(subsection 2.2) surfaces of discontinuity before presenting our algorithm, which is applicablein the presence of an arbitrary number of surfaces of discontinuity (subsection 2.3). If F is continuously differentiable on thetrajectory x , the derivative δx t = Dφ t ( x ; δx ) satisfies the linear time-varying variationalequation [27, Appendix B](2.4) ∀ t ≥ δx t = (cid:90) t DF ( x τ ) · δx τ dτ, whence δx t = Dφ t ( x ; δx ) can be approximated to any desired precision in polynomial timeby applying numerical simulation algorithms [27, Ch. 4] to (2.1), (2.4). If F is continuously differentiable everywhere excepta smooth codimension-1 submanifold H ⊂ R d that intersects the trajectory x transversallyat only one point x s , s ∈ (0 , t ), the continuous-time equation (2.4) is augmented by thediscrete-time update [2, Eqn. (58)],(2.5) δx + s = (cid:18) I d + ( F + − F − ) · η (cid:62) η (cid:62) · F − (cid:19) · δx − s = M · δx − s , where δx ± s = lim τ → s ± δx τ and F ± = lim τ → s ± F ( x τ ) denote the limiting values of δx τ and F ( x τ ) at s from the right (+) and left ( − ) and η ∈ R d is any vector orthogonal to surface H at x s ; M ∈ R d × d is termed the saltation matrix [10, Eqn. (2.76)], [22, Eqn. (7.65)]. Overall,the desired derivative is(2.6) Dφ t ( x ; δx ) = Dφ t − s ( x s ) · M · Dφ s ( x ) · δx , where Dφ t − s ( x s ) , Dφ s ( x ) ∈ R d × d can be approximated by simulating (2.2), (2.4) since theflow is smooth away from time s . Computing the saltation matrix M requires O (cid:0) d (cid:1) timeand space, but evaluating its action on δx − s in (2.5) requires only O ( d ) time and space. If F is continuously differentiable everywhereexcept a finite set of smooth codimension-1 submanifolds { H j } nj =1 that intersect the trajectory x transversally at only one point x s (see Figure 2.1(a) for an illustration when n = 2), s ∈ (0 , t ),we showed in [6, Eqn. (65)] that the discrete-time update (2.5) is applied once for each surface.However, the order in which the updates are applied, and the limiting values of the vector fieldused to determine each update’s saltation matrix, depend on δx . If the surfaces intersecttransversally, there are n ! different saltation matrices determined by 2 n vector field values, soconsidering all update orders requires factorial time and space. To make these observationsprecise and specify the notation employed in Figures 2.1 and 2.2, we formally define the classof nonsmooth vector fields considered in this paper [6, Defs. 1, 2]: Definition 2.1. (event-selected C r ( EC r ) vector field) A vector field F : D → T D definedon an open domain D ⊂ R d is event-selected C r with respect to h ∈ C r ( U, R n ) at ρ ∈ R d if U ⊂ D is an open neighborhood of ρ and: (event functions) there exists f > such that Dh ( x ) · F ( x ) ≥ f for all x ∈ U ;This manuscript is for review purposes only. G. COUNCIL, S. REVZEN, S. A. BURDEN H H x δx δx − s δx + s δx t ρ − ρ + δρ − = δx − s δρ + = δx + s (cid:101) H x s (a) continuous-time variational dynamics (2.4) x t ρ = x s δx + s = δρ + = B ( δρ − ) = B ( δx − s )(b) discrete-time variational dynamics (2.8) (cid:101) H F ( − , +1) ( ρ ) F + ( ρ ) F − ( ρ ) F (+1 , − ( ρ ) δ ˙ x τ = DF ( x τ ) · δx τ Figure 2.1.
Variational dynamics that determine the B-derivative of an EC r vector field’s P C r flow (2.8). (a) Vector field F : R → T R is smooth everywhere except the smooth codimension-1 submanifolds H , H ⊂ R that intersect transversally at x s ∈ R , generating a piecewise-differentiable flow φ : [0 , ∞ ) × R → R satisfying φ τ ( x ) = x τ for all τ ∈ [0 , t ] , i.e. F is EC r and φ is P C r [6]. The B-derivative Dφ t ( x ; δx ) = δx t isdetermined as in (2.10) by the continuous-time variational dynamics δ ˙ x τ = DF ( x τ ) · δx τ and the discrete-timevariational dynamics δx + s = B ( δx − s ) . The algorithms in Figure evaluate the piecewise-linear function B using the auxiliary nonsmooth system in (b) determined by the tangent planes (cid:101) H , (cid:101) H and vector field limits F b ( ρ ) in (2.9) for b ∈ { ( − , (+1 , − , ( − , +1) , + } = {− , +1 } . (smooth extension) for all b ∈ {− , +1 } n = B n , with (2.7) D b = { x ∈ U : b j ( h j ( x ) − h j ( ρ )) ≥ } ,F | Int D b admits a C r extension F b : U → T U . Our algorithms in Figure 2.2 compute(2.8) δx + s = δρ + = B ( δρ − ) = B ( δx − s )given δρ − = δx − s ∈ R d , normals { η j = Dh j ( ρ ) } nj =1 ⊂ R d at x s to surfaces (cid:110) H j = h − j ( ρ ) (cid:111) nj =1 ,and a function Γ : {− , +1 } n → R d that evaluates limits of the vector field F at ρ = x s ,(2.9) ∀ b ∈ {− , +1 } n : Γ( b ) = F b ( ρ ) , using the piecewise-constant dynamics illustrated in Figure 2.1(b), which are the discrete-timeanalog of the continuous-time variational dynamics (2.4). Overall, the desired derivative is(2.10) Dφ t ( x ; δx ) = Dφ t − s ( x s ) · B ( Dφ s ( x ) · δx ) , This manuscript is for review purposes only. -DERIVATIVE OF AN EC r VECTOR FIELD’S
P C r FLOW 5
Algorithm 2.1 δρ + ← B ( δρ − , η, Γ) δt ← ∈ R δρ + ← δρ − ∈ R d b ← − ∈ {− , +1 } n while b (cid:54) = + do for j ∈ { , . . . , n } do τ j ← − (cid:16) η (cid:62) j · δρ + (cid:17) / (cid:16) η (cid:62) j · Γ( b ) (cid:17) j ∗ ← arg min j ∈{ ,...,n } { τ j : b j < } δt ← δt + τ j ∗ δρ + ← δρ + + τ j ∗ · Γ( b ) b j ∗ ← +1 return δρ + − δt · Γ(+ ) Algorithm 2.2 def B(dx,e,G): dt = 0 dx = np.array(dx) b = -np.ones(len(e),dtype=np.int) while np.any(b < 0) : tau = -np.dot(e,dx)/np.dot(e,G(b)) tau[b > 0] = np.inf j = np.argmin(tau) dt += tau[j] dx += tau[j] * G(b) b[j] = +1 return dx - dt * G(b) Figure 2.2.
Algorithms that evaluate the B-derivative of an EC r vector field’s P C r flow written inpseudocode (Algorithm ) and Python [30] sourcecode (Algorithm ; requires import numpy as np [25]).These algorithms apply at a point ρ ∈ R d where a vector field F : R d → T R d is event-selected C r with respectto n surfaces (see Figure for an illustration when d = n = 2 ), and assume the following data is given:tangent direction,surface normals at ρ ,vector field limits (2.9) , δρ − ∈ T ρ R d , η = { η j } nj =1 ⊂ R d , Γ : {− , +1 } n → R d , dx – array , dx.shape == (d,) ; e – array , e.shape == (n,d) ; G – function , G(b).shape == (d,) . where B : T ρ R d → T ρ R d is the continuous piecewise-linear function defined by our algorithmsin Figure 2.2. Our algorithms require O (cid:0) n d (cid:1) time and O ( d ) space to evaluate the directionalderivative (2.3) .Assuming for the moment that these algorithms are correct, we emphasize that theyachieve a dramatic reduction in the computational complexity of evaluating the B-derivative– from factorial to low-order polynomial – relative to na¨ıve enumeration of all pieces of theB-derivative. However, despite the apparent simplicity of our algorithms (computationallyand conceptually), verifying their correctness requires significant technical effort; the bulk ofthe present paper is devoted to this verification task.
3. Applications.
To illustrate and validate our methods, we apply the algorithm fromthe preceding section to piecewise-constant vector fields in subsection 3.1 and mechanicalsystems subject to unilateral constraints in subsection 3.2. Sourcecode implementation ofAlgorithm 2.2 and applications from the remainder of this section are provided in SM.
Consider the vector field F : R d → T R d defined by(3.1) ˙ x = F ( x ) = + ∆ (sign( x )) These algorithms can be modified as in (6.9) to determine the order of surface crossings for the perturbedtrajectory without changing the time or space complexity, so the associated saltation matrix (6.4) can beconstructed in O (cid:0) nd (cid:1) time and O (cid:0) d (cid:1) space; this construction is discussed in more detail in section 6. This manuscript is for review purposes only.
G. COUNCIL, S. REVZEN, S. A. BURDEN F ( − , = (1 , ) F (1 , = ( , ) F ( − , − = (1 , F (1 , − = ( , H H F ( − , = (1 , ) F (1 , = (1 , ) F ( − , − = (1 , F (1 , − = ( , H H Figure 3.1.
B-derivative of vector field from subsection 3.1 in linear (left) and piecewise-linear (right) cases.
The vector field F defined in (3.1) is piecewise-constant and discontinuous across the coordinate hyperplanes H , H , generating a piecewise-differentiable flow φ with B-derivative B . (left) The B-derivative is linear inthe special case defined by (3.2) . (right) The B-derivative is continuous and piecewise-linear in general, so aball of initial conditions flows to a piecewise-ellipsoid (gold and green fill). where ∆ : B d → R d ; so long as all components of all vectors specified by ∆ are larger than −
1, i.e. min b ∈ B d [∆( b )] j > − F is event-selected C ∞ with respect to the identity function h : R d → R d defined by h ( q ) = q . We regard (3.1) as a canonical form for piecewise-constantevent-selected C ∞ vector fields that are discontinuous across d subspaces, since any such vectorfield can be obtained by applying a linear change-of-coordinates to (3.1). In what follows, wefocus on the trajectory that passes through the origin ρ = 0, which lies at the intersection of d surfaces of discontinuity for F . With ρ − = ρ − F − ( ρ ), ρ + = ρ + F + ( ρ ), we note that ρ − flows to ρ + through ρ in 1 (one) unit of time.Our goal is to compute D x φ (1 , ρ − ; δρ − ) ∈ T ρ + R d for a given δρ − ∈ T ρ − R d . In the generalcase, the desired derivative is piecewise-linear with (up to) d ! distinct pieces, providing ageneral test. In the special case where ∆( b ) = − δ · b for all b ∈ B d , | δ | <
1, the desiredderivative is linear [6, Eqn. (86)],(3.2) D x φ (1 , ρ − ; δρ − ) = 1 − δ δ · δρ − , providing a closed-form expression for comparison. Figure 3.1 illustrates results from bothcases with d = 2; a more exhaustive test suite is provided in SM. Consider a mechanical sys-tem whose configuration is subject to one-sided (i.e. unilateral ) constraints. The dynamicsof such systems have been studied extensively using the formalisms of complementarity [24,Sec. 3], measure differential inclusions [3, Sec. 3], hybrid systems [21, Sec. 2.4, 2.5], and geo-metric mechanics [13, Sec. 3]. Regardless of the chosen formalism, in a coordinate chart
This manuscript is for review purposes only. -DERIVATIVE OF AN EC r VECTOR FIELD’S
P C r FLOW 7 Q ⊂ R d the dynamics governing q take the form (3.3) M ( q )¨ q = f ( q, ˙ q ) subject to a ( q ) ≥ M ( q ) ∈ R d × d specifies the kinetic energy metric; f ( q, ˙ q ) ∈ R d specifies the internal,applied, and Coriolis forces; a ( q ) ∈ R n specifies the unilateral constraints ; and we assumein what follows that M , f , and a are smooth functions. Different formalisms enforce theconstraint a ( q ) ≥ If constraints are enforced rigidly as in [3, 21, 24], meaning that they must be satisfied exactly, then the velocity must undergoimpact (i.e. change discontinuously) whenever ˙ q ∈ T q Q is such that a j ( q ) = 0 and Da j ( q ) · ˙ q < j ∈ { , . . . , n } [24, Sec. 2] [21, Eqn. (23)] [3, Eqn. (23)]. Unfortunately for ourpurposes, these discontinuities in the state vector x = ( q, ˙ q ) cannot be modeled using an event-selected C r vector field ˙ x = F ( x ), and the flow of such systems is generally discontinuous . C flows. We nowconsider the formalism in [13] that “softens” (i.e. approximately enforces) rigid constraints a ( q ) ≥ penalty functions { v j } nj =1 that scalequadratically with the degree of constraint violation [13, Eqn. (12)],(3.4) ∀ j ∈ { , . . . , n } : v j ( q ) = (cid:40) , a j ( q ) ≥ κ j a j ( q ) , a j ( q ) < a j ( q ) ≥ κ j , leading tothe unconstrained dynamics [13, Eqn. (14)](3.5) M ( q )¨ q = f ( q, ˙ q, u ) − n (cid:88) j =1 Dv j ( q ) (cid:62) = f ( q, ˙ q, u ) − (cid:88) (cid:110) ( κ j a j ( q )) · Da j ( q ) (cid:62) : j ∈ { , . . . , n } , a j ( q ) < (cid:111) . As shown by [39, Thm. 3], trajectories of (3.5) converge to those of (3.3) in the rigid limit(i.e. as stiffnesses go to infinity). Importantly for our purposes, the dynamics in (3.5) can bemodeled using an event-selected vector field along trajectories that pass transversally throughthe constraint surfaces, whence our algorithms can compute the B-derivative of the flow.However, the vector field (3.5) in this case is (locally Lipschitz) continuous, hence the B-derivative is trivial (all non-identity terms in (6.4) are zero), whence the flow is continuously-differentiable ( C ). We interpret the inequality a ( q ) ≥ We note that the flow can be
P C r at non-impact times if the constraint surfaces intersect orthogonally [26],i.e. if the surface normals are orthogonal with respect to the inverse of the kinetic energy metric [3, Theorem 20]. This manuscript is for review purposes only.
G. COUNCIL, S. REVZEN, S. A. BURDEN ψ ( x, y ) θ δ ˙ θ + δ ˙ θ − δ ˙ θ − g Figure 3.2.
Vertical-plane biped, a mechanical system subject to unilateral constraints (subsection 3.2.4), consists of a body with two rigid massless legs falling under the influence of gravity toward a substrate. Thesystem’s flow can be C (left column) or P C r (right) depending on how forces vary as limbs contact substrate. EC r vector fields, C flows. We now augmentthe unconstrained dynamics (3.5) with dissipation as in [13]:(3.6) M ( q )¨ q = f ( q, ˙ q, u ) − (cid:88) (cid:110) ( κ j a j ( q ) + β j Da j ( q ) · ˙ q ) · Da j ( q ) (cid:62) : j ∈ { , . . . , n } , a j ( q ) < (cid:111) ;in essence, each constraint penalty is augmented by a spring-damper that is only active whenthe constraint is violated as in studies involving contact with complex geometry [12] or ter-rain [1]. The dynamics in (3.6) can be modeled using an event-selected vector field along trajec-tories that pass transversally through the constraint surfaces, and the vector field is discontinu-ous along the constraint surfaces. However, we can show that the flow of (3.6) is continuously-differentiable ( C ) along any trajectory that passes transversally through constraint surfaces.Indeed, letting x = ( q, ˙ q ) denote the state of the system so that ˙ x = ( ˙ q, ¨ q ) = F ( x ) is determinedby (3.6), the saltation matrix (2.5) associated with each constraint a j has the form(3.7) I + 1 Da j ( q ) · ˙ q (cid:20) ± ( κ j a j ( q ) + β j Da j ( q ) · ˙ q ) · Da j ( q ) (cid:62) (cid:21) (cid:2) Da j ( q ) 0 (cid:3) where the sign in the column vector is determined by whether the constraint is activating( − ) or deactivating (+). Since matrices of the form in (3.7) commute, the saltation matricesassociated with simultaneous activation and/or deactivation of multiple constraints are allequal, whence the flow of (3.6) is continuously-differentiable ( C ) along any trajectory thatpasses transversally through constraint surfaces. To ground the preceding observations, we con-sider the vertical-plane biped illustrated in Figure 3.2( left ) that falls under the influence ofgravity toward a substrate. The biped body has mass m and moment-of-inertia J ; we let( x, y ) ∈ R denote the position of its center-of-mass in the plane and θ ∈ S denote its rota-tion. Two rigid massless limbs of length (cid:96) protrude at an angle of ± ψ with respect to verticalfrom the body’s center-of-mass above a smooth substrate whose height is a quadratic functionof horizontal position, yielding unilateral constraints(3.8) a ( x, y, θ ) = − y − ( x + (cid:96) cos ( θ − ψ )) − (cid:96) sin( θ − ψ ) ,a ( x, y, θ ) = − y − ( x + (cid:96) cos ( θ + ψ )) − (cid:96) sin( θ + ψ ) . This manuscript is for review purposes only. -DERIVATIVE OF AN EC r VECTOR FIELD’S
P C r FLOW 9
We consider the smoothness of the system’s flow along a trajectory that activates both con-straints simultaneously . Direct calculation shows that adopting the formalism in (3.6) yieldscontinuously-differentiable flow for this system as illustrated in Figure 3.2( middle ).To obtain a flow that is piecewise-differentiable but not continuously-differentiable, wemodify the damping coefficients in (3.6) using the following logic : β = β = if a ( q ) < a ( q ) ≥ exclusive or); β = β = 1 if a ( q ) < a ( q ) <
0. Direct calculation showsthat the saltation matrices obtained from different sequences of constraint activations (leftfoot reaches substrate before right foot or vice-versa) are distinct:(3.9) M (left , right) − M (right , left) = − β cos( ψ ) 0 − β (sin(2 ψ ) + cos( ψ )) 0 0 00 0 0 0 0 0 . The piecewise-linear B-derivative of the system’s flow is illustrated in Figure 3.2( right ).
4. Background.
To verify correctness of the algorithms specified in section 2, we utilizethe representation of piecewise-affine functions from [15], elements of the theory of piecewise-differentiable functions from [36], and results about the class of nonsmooth flows under con-sideration from [6]. In an effort to make this paper self-contained (i.e. to save the reader fromneeding to cross-reference multiple citations to follow our derivations), we include a substan-tial amount of background details in this section. The expert reader may wish to skim or skipthis section, returning only if questions arise in subsequent sections.
We let 0 d ∈ R d denote the vector of zeros, n ∈ R n the vector ofones, and I d ∈ R d × d the identity matrix; when dimensions are clear from context, we suppresssubscripts. The vectorized signum function sign : R d → {− , +1 } d is defined by(4.1) ∀ x ∈ R d , j ∈ { , . . . , d } : [sign( x )] j = sign( x j ) = (cid:40) − , x j < , x j ≥ . e.g. initial condition (cid:16) ( x , y , θ ) , ( ˙ x , ˙ y , ˙ θ ) (cid:17) = ((0 , h, , (0 , , h is the initial body height Sourcecode that verifies this fact using a computer algebra system is provided in SM. Although we introduce this logic purely for illustrative purposes, we note that non-trivial dependence offorcing on the set of active constraints could be implemented physically using clutches [8] or actuators [38].
This manuscript is for review purposes only. If A ∈ R (cid:96) × m and B ∈ R m × n then A · B ∈ R (cid:96) × n denotes matrix multiplication. Given a subset S ⊂ R d , we define [36, Sec. 2.1.1]aff S = n (cid:88) j =1 α j v j : n ∈ N , { v j } nj =1 ⊂ S, { α j } nj =1 ⊂ R , n (cid:88) j =1 α j = 1 , (4.2a) cone S = n (cid:88) j =1 α j v j : n ∈ N , { v j } nj =1 ⊂ S, { α j } nj =1 ⊂ [0 , ∞ ) , (4.2b) conv S = n (cid:88) j =1 α j v j : n ∈ N , { v j } nj =1 ⊂ S, { α j } nj =1 ⊂ [0 , , n (cid:88) j =1 α j = 1 , (4.2c)termed the affine span , cone span , and convex hull of S , respectively. The dimension of aconvex set S is defined to be the dimension of its affine span, dim S = dim aff S . A nonemptyset S ⊂ R d is called a polyhedron [36, Sec. 2.1.2] if there exists A ∈ R m × d , b ∈ R m suchthat S = (cid:8) x ∈ R d : A · x ≤ b (cid:9) ; note that S is closed and convex. The linear subspace L = (cid:8) x ∈ R d : A · x = 0 (cid:9) is called the lineality space of S . We will represent a piecewise-affine function using a triangulation ( Z − , Z + , ∆) [15, Sec. 3.1] that consists of a combinatorial simplicial complex ∆whose vertex set is in 1-to-1 correspondence with each of the finite sets of vectors Z − ⊂ R d , Z + ⊂ R c . For our purposes, a combinatorial simplicial complex ∆ is a collection of finite sets∆ = { ∆ ω } ω ∈ Ω such that S ⊂ ∆ ω = ⇒ S ∈ ∆ for all ω ∈ Ω; we call (cid:83) ω ∈ Ω ∆ ω the vertex set of∆. We assume that, for every ω ∈ Ω, the collections of vectors Z ± ω ⊂ Z ± determined by ∆ ω are affinely independent [15, Sec. 2.1.1] so that ∆ ± ω = conv Z ± ω are ( ω ) − − ω ⊂ R d , ∆ + ω ⊂ R c . We assume further that, forevery ω, ω (cid:48) ∈ Ω, the collections of vectors Z ± ω,ω (cid:48) ⊂ Z ± determined by ∆ ω ∩ ∆ ω (cid:48) coincide with Z ± ω ∩ Z ± ω (cid:48) ⊂ Z ± so that ∆ ± = { ∆ ± ω } ω ∈ Ω are geometric simplicial complexes [15, Sec. 2.2.1].With these assumptions in place, the correspondence between Z − and Z + determined by thetriangulation ( Z − , Z + , ∆) uniquely defines a piecewise-affine function P : | ∆ − | → | ∆ + | usingthe construction from [15, Sec. 3.1] where | ∆ − | = (cid:83) ω ∈ Ω ∆ − ω ⊂ R d , | ∆ + | = (cid:83) ω ∈ Ω ∆ + ω ⊂ R c aretermed the carriers [36, Sec. 2.2.1] of the geometric simplicial complexes ∆ ± . If a piecewise-affine function P : R d → R c is positivelyhomogeneous , that is,(4.3) ∀ α ≥ , v ∈ R d : P ( α · v ) = α · P ( v ) , then P is piecewise-linear [36, Prop. 2.2.1]. In this case, P admits a conical subdivision [36,Prop. 2.2.3], that is, there exists a finite collection Σ = { Σ ω } ω ∈ Ω such that: (i) Σ ω ⊂ R d is a d - There are more general definitions of ([complete] semi-)simplicial complexes and the closely-related conceptof ∆-complexes in the literature [16, Ch. 2.1], [15, App. A.3.1]. Since we employ these concepts primarily inservice of parameterizing piecewise-affine functions as in [15, Sec. 3.1], we adopt the (relatively restrictive)definitions of combinatorial and geometric simplicial complexes from [15, Sec. 2.2.1] in what follows.
This manuscript is for review purposes only. -DERIVATIVE OF AN EC r VECTOR FIELD’S
P C r FLOW 11 dimensional polyhedral cone for each ω ∈ Ω; (ii) the Σ ω ’s cover R d ; and (iii) the intersectionΣ ω ∩ Σ ω (cid:48) is either empty or a proper face of both polyhedral cones for each ω, ω (cid:48) ∈ Ω. P C r ) functions. (This section is largely repeated from [6,Sec. 3.2].) The notion of piecewise–differentiability we employ was originally introduced in [32];since the monograph [36] provides a more recent and comprehensive exposition, we adopt thenotational conventions therein. Let r ∈ N ∪ {∞} and D ⊂ R d be open. A continuous function f : D → R c is called piecewise- C r if for every x ∈ D there exists an open set U ⊂ D containing x and a finite collection { f j : U → R c } j ∈ J of C r functions such that for all x ∈ U we have f ( x ) ∈ { f j ( x ) } j ∈ J . The functions { f j } j ∈ J are called selection functions for f | U , and f is said to be a continuous selection of { f j } j ∈ J on U . A selection function f j is said to be active at x ∈ U if f ( x ) = f j ( x ). We let P C r ( D, R c ) denote the set of piecewise- C r functions from D to R c . Note that P C r is closed under composition. The definition of piecewise- C r may at firstappear unrelated to the intuition that a function ought to be piecewise-differentiable preciselyif its “domain can be partitioned locally into a finite number of regions relative to whichsmoothness holds” [33, Section 1]. However, as shown in [33, Thm. 2], piecewise- C r functionsare always piecewise-differentiable in this intuitive sense.Piecewise-differentiable functions possess a first–order approximation Df : T D → T R c called the Bouligand derivative (or B–derivative) [36, Ch. 3]; this is the content of [36,Lemma 4.1.3]. Significantly, this B–derivative obeys generalizations of many techniques fa-miliar from calculus, including the Chain Rule [36, Thm 3.1.1], Fundamental Theorem ofCalculus [36, Prop. 3.1.1], and Implicit Function Theorem [31, Cor. 20]. We let Df ( x ; δx )denote the B–derivative of f evaluated on the tangent vector δx ∈ T x D . The B-derivativeis positively homogeneous, i.e. ∀ δx ∈ T x D, λ ≥ Df ( x ; λ δx ) = λDf ( x ; δx ), and coincideswith the directional derivative of f in the δx ∈ T x D direction. In addition, the B-derivative Df ( x ) : T x D → T f ( x ) R c of f at x ∈ D is a continuous selection of the derivatives of theselection functions active at x [36, Prop. 4.1.3],(4.4) ∀ δx ∈ T x D : Df ( x ; δx ) ∈ { Df j ( x ) · δx } j ∈ J . However, the function Df is generally not continuous at ( x, δx ) ∈ T D ; if it is, then f is C at x [36, Prop. 3.1.2]. C r ( EC r ) vector fields and their P C r flows. Vector fields withdiscontinuous right-hand-sides and their associated flows have been studied extensively [14].In Definition 2.1 [6, Defs. 1, 2], a special class of so-called event-selected C r ( EC r ) vector fieldswere defined which are allowed to be discontinuous along a finite number of codimension-1submanifolds but do not exhibit sliding [20] along these submanifolds, and are C r elsewhere.Importantly, as shown in [6, Thm. 5], an event-selected C r vector field F : R d → T R d generatesa piecewise-differentiable flow, that is, there exists a function φ : F → R d that is piecewise- C r ( φ ∈ P C r ) in the sense defined in [36, Sec. 4.1] (summarized in subsection 4.4) where i.e. Σ ω = (cid:110)(cid:80) (cid:96) ω j =1 α j v ωj : { α j } (cid:96) ω j =1 ⊂ [0 , ∞ ) (cid:111) , some { v j } (cid:96) ω j =1 ⊂ R d [36, Thm. 2.1.1], and dim Σ ω = d i.e. (cid:83) ω ∈ Ω Σ ω = R d i.e. Σ ω ∩ Σ ω (cid:48) = (cid:110)(cid:80) (cid:96) ω,ω (cid:48) j =1 α j v ω,ω (cid:48) j : { α j } (cid:96) ω,ω (cid:48) j =1 ⊂ [0 , ∞ ) (cid:111) , some (cid:110) v ω,ω (cid:48) j (cid:111) (cid:96) ω,ω (cid:48) j =1 ⊂ (cid:8) v ωj (cid:9) (cid:96) ω j =1 ∪ (cid:110) v ω (cid:48) j (cid:111) (cid:96) ω (cid:48) j =1 This manuscript is for review purposes only. F ⊂ R × R d and(4.5) ∀ ( t, x ) ∈ F : φ ( t, x ) = x + (cid:90) t F ( φ ( s, x )) ds. Since φ is P C r , it admits a first-order approximation Dφ : T F → T R d termed the Bouligand (or B- )derivative [36, Sec. 3.1], which is a continuous piecewise-linear function of tangentvectors at every ( t, x ) ∈ F , that is, the directional derivative Dφ ( t, x ) : T ( t,x ) F → T φ ( t,x ) R d iscontinuous and piecewise-linear for all ( t, x ) ∈ F . EC r vector field’s P C r flow. Suppose F : R d → T R d is an EC r vector field with P C r flow φ : F → R d . Given a tangent vector ( δt, δx ) ∈ T ( t,x ) F , it was shownin [6, Sec. 7.1.4] that the value of the B-derivative Dφ ( t, x ; δt, δx ) ∈ T φ ( t,x ) R d can be obtainedby solving a jump-linear-time-varying differential equation [6, Eqn. (70)], where the “jump”arises from a matrix Ξ ω determined by the sequence ω in which the perturbed initial state x + α δx crosses the surfaces of discontinuity of the vector field F for small α > Dφ ( t, x ) (and,to the best of our knowledge, neither has subsequent work). The key theoretical contributionof this paper, obtained in section 5, is a representation of the B-derivative with respect tostate, D x φ ( t, x ), using a triangulation of its domain and codomain as defined in [15, Sec. 3.1](and recalled in subsection 4.2).To inform the triangulation of the B-derivative D x φ ( t, x ), we recall the values it takeson. Since the flow φ : F → R d is piecewise- C r ( P C r ), it is a continuous selection of a finitecollection of C r functions (cid:8) φ ω : F ω → R d (cid:9) ω ∈ Ω near ( t, x ) ∈ F , where F ω ⊂ F is an open setcontaining ( t, x ) for each ω ∈ Ω [36, Sec. 4.1], and the B-derivative D x φ ( t, x ) is a continuousselection of the classical ( Fr´echet or F- )derivatives { D x φ ω ( t, x ) } ω ∈ Ω [36, Prop. 4.1.3], that is,(4.6) ∀ δx ∈ W ω ⊂ T x R d : D x φ ( t, x ; δx ) = D x φ ω ( t, x ) · δx, where W ω ⊂ T x R d is the subset of tangent vectors where the selection function D x φ ω is essentially active [36, Prop. 4.1.1]. If s, t ∈ R and x ∈ R d are such that 0 < s < t and thevector field F is C r on φ ([0 , t ] \ { s } , x ), i.e. the trajectory initialized at x ∈ R d encountersexactly one discontinuity of F at ρ = φ ( s, x ) on the time interval [0 , t ], then D x φ ω ( t, x ) hasthe form(4.7) D x φ ω ( t, x ) = D x φ ( t − s, ρ ) · (cid:2) F + ( ρ ) I d (cid:3) · Ξ ω · (cid:20) (cid:62) d I d (cid:21) · D x φ ( s, x )where F + is the C r extension of F | Int D + that exists by virtue of condition 2 in Def. 2.1 andΞ ω ∈ R ( d +1) × ( d +1) is the matrix from [6, Eqn. (67)] corresponding to the selection functionindex ω ∈ Ω. In what follows, we will work in circumstances where the selection functions areindexed by the symmetric permutation group over n elements, i.e. Ω = S n , and combine (4.6)and (4.7) as(4.8) ∀ δx ∈ W σ ⊂ T x R d : D x φ ( t, x ; δx ) = D x φ ( t − s, ρ ) · M σ · D x φ ( s, x ) · δx This manuscript is for review purposes only. -DERIVATIVE OF AN EC r VECTOR FIELD’S
P C r FLOW 13 where the saltation matrix M σ ∈ R d × d corresponding to index σ is defined by(4.9) M σ = (cid:2) F + ( ρ ) I d (cid:3) · Ξ σ · (cid:20) (cid:62) d I d (cid:21) . EC r vector field. Suppose vector field F : R d → T R d isevent-selected C r with respect to h ∈ C r ( U, R n ) at ρ ∈ U ⊂ R d . For b ∈ B n = {− , +1 } n let(4.10) (cid:101) D b = (cid:110) x ∈ R d : b j Dh j ( ρ )( x − ρ ) ≥ (cid:111) and consider the piecewise-constant vector field (cid:101) F : R d → T R d defined by(4.11) ∀ b ∈ B n , x ∈ (cid:101) D b : (cid:101) F ( x ) = F b ( ρ )where F b is the C r extension of F | Int D b that exists by virtue of condition 2 in Def. 2.1 Notethat (cid:101) F is event-selected C r with respect to the affine function (cid:101) h defined by(4.12) ∀ x ∈ R d : (cid:101) h ( x ) = Dh ( ρ )( x − ρ ) , whence it generates a piecewise-differentiable flow (cid:101) φ : (cid:101) F → R d where (cid:101) F = R × R d . In [6,Sec. 7.1.3], (cid:101) F was referred to as the sampled vector field since it is obtained by “sampling”the selection functions F b that define F near ρ , and it was noted that the function (cid:101) φ ispiecewise-affine and it approximates the original vector field’s flow φ near ρ . We will leveragethe algebraic properties of (cid:101) φ and its relationship to φ in what follows to obtain our results. EC r vector field and its local approximation. Supposevector field F : R d → T R d is event-selected C r with respect to h ∈ C r ( U, R n ) at ρ ∈ U ⊂ R d ,and let φ ∈ P C r ( F , R d ) be its piecewise-differentiable flow. Then [6, Thm. 7] ensures thereexists a piecewise-differentiable time-to-impact function τ ∈ P C r ( U, R n ) for which(4.13) ∀ x ∈ U, j ∈ { , . . . , n } : φ ( τ j ( x ) , x ) ∈ H j = h − j ( h j ( ρ )) , i.e. the point x flows to the surface H j in time τ j ( x ). Similarly, applying [6, Thm. 7] to thesampled vector field (cid:101) F : R d → T R d and piecewise-affine flow (cid:101) φ : (cid:101) F → R d associated with F at ρ constructed in subsection 4.7 ensures there exists a piecewise-affine time-to-impact function (cid:101) τ : R d → R n for which(4.14) ∀ x ∈ R d , j ∈ { , . . . , n } : (cid:101) φ ( (cid:101) τ j ( x ) , x ) ∈ (cid:101) H j = ρ + ker Dh j ( ρ ) , i.e. the point x flows to the affine subspace (cid:101) H j in time (cid:101) τ j ( x ). Ξ σ ∈ R ( d +1) × ( d +1) is referred to as a saltation matrix in [6, Sec. 7.1.4], but this usage is inconsistent withthe original definition of M σ ∈ R d × d as the saltation matrix in [2]. Note that (cid:101) F is well-defined as the value of F b is uniquely determined at ρ by virtue of being continuous,even though the original F is undefined at ρ . This manuscript is for review purposes only.
5. Representation.
Our main theoretical result is an explicit representation for the Bouli-gand (or B-)derivative of the piecewise-differentiable flow generated by an event-selected C r vector field. To that end, let F : R d → T R d be an event-selected C r vector field and φ : F → R d its piecewise-differentiable flow. In what follows, we will assume that s, t ∈ R and x ∈ R d are such that 0 < s < t and the vector field F is C r on φ ([0 , t ] \ { s } , x ). Although a gen-eral trajectory can encounter more than one point of discontinuity for F , such points areisolated [6, Lem. 6], so the Chain Rule for B-differentiable functions [36, Thm. 3.1.1] can beapplied to triangulate the desired flow derivative by composing the triangulated flow deriv-atives associated with each point. Thus, without loss of generality, we restrict our attentionto portions of trajectories that encounter one point of discontinuity for F , which point liesat the intersection of n surfaces of discontinuity for F . We assume n > n = 1 the desired B-derivative islinear [2], so it may be represented and employed in computations as a matrix.The B-derivative D x φ ( t, x ) : T x R d → T φ ( t,x ) R d we seek is a continuous piecewise-linearfunction, so it can be parsimoniously represented using a triangulation [15, Sec. 3.1], that is,a combinatorial simplicial complex (as defined in subsection 4.2) each of whose vertices areassociated with a pair of (tangent) vectors – one each in the domain and codomain of D x φ ( t, x ).We will obtain this triangulation via an indirect route: in subsection 5.1, we triangulate thepiecewise-affine flow (cid:101) φ introduced in subsection 4.7; in subsection 5.2, we differentiate ourrepresentation of (cid:101) φ to obtain a triangulation of the B-derivative D x (cid:101) φ ; in subsection 5.3, weshow how the B-derivative D x φ can be obtained from D x (cid:101) φ , providing a triangulation of thedesired derivative. The goal of this subsection is to triangulate the piecewise-affine flow (cid:101) φ introduced in subsection 4.7. To that end, let ρ = φ ( s, x ) and suppose rank Dh ( ρ ) = n so (cid:8) δρ ∈ T ρ R d : b = sign Dh ( ρ ) · δρ (cid:9) has nonempty interior for each b ∈ {− , +1 } n = B n . Letting K = ker Dh ( ρ ) ⊂ T ρ R d denote the kernel of Dh ( ρ ) and K ⊥ its orthogonal complement, foreach b ∈ B n there exists a unique ζ b ∈ K ⊥ + { ρ } such that(5.1) Dh b> ( ρ )( ζ b − ρ ) = 0 and Dh b< ( ρ )( ζ b + F b ( ρ ) − ρ ) = 0where h b> (respectively, h b< ) denotes the function obtained by selecting components h j of h for which b j = +1 (respectively, b j = − (cid:101) φ introduced in subsection 4.7 (see Figure 5.1(a)):(5.2) ∀ b ∈ B n : ζ b ∈ (cid:101) D − , (cid:101) φ (1 , ζ b ) = ζ b + F b ( ρ ) ∈ (cid:101) D + , that is, the point ζ b lies “before” all event surface tangent planes and flows in 1 (one) unitof time to ζ b + F b ( ρ ) which lies “after” all event surface tangent planes (neither “before” nor“after” should be interpreted strictly). We denote the collections of these vectors as follows:(5.3) Z − = { ζ b } b ∈ B n , Z + = { ζ b + F b ( ρ ) } b ∈ B n . As observed in [6, Sec. 7.1.5], first-order approximations of an EC r vector field’s P C r flow are not affectedby flow between surfaces that are tangent at ρ , so we assume such redundancy has been removed. Here and in what follows we mildly abuse notation via the natural vector space isomorphism R d (cid:39) T ρ R d . Uniqueness is ensured by rank Dh ( ρ ) = n since (i) K ⊥ is n -dimensional, (ii) the rows of Dh ( ρ ) are linearlyindependent, and hence (iii) there are n independent equations in the n unknowns needed to specify ζ b in (5.1). This manuscript is for review purposes only. -DERIVATIVE OF AN EC r VECTOR FIELD’S
P C r FLOW 15 (cid:101) φ ζ − (cid:101) H (cid:101) H ζ ( − , +1) ζ (+1 , − (cid:101) φ (cid:0) ζ (+1 , − (cid:1) ρ = ζ + (cid:101) φ ( ζ + ) (cid:101) φ (cid:0) ζ ( − , +1) (cid:1) ∆ − (2 , ∆ − (1 , ∆ +(2 , ∆ +(1 , ∆ − σ K Σ σ = ∆ − σ + K (a) flow of the sampled system (b) triangulation of | ∆ − | , | ∆ + | (c) triangulation of | Σ | ρ − ρ + Figure 5.1.
Triangulation of the time-1 flow (cid:101) φ of the sampled system associated with an EC r vector field. (a) For each b ∈ {− , +1 } , the point ζ b defined by (5.1) flows from (cid:101) D − to (cid:101) D + in 1 (one) unit of time via the sampled system illustrated in Figure (b) and defined in subsection . (b) The sets (cid:8) ζ − , ζ + , ζ (+1 , − (cid:9) , (cid:8) ζ − , ζ + , ζ ( − , +1) (cid:9) indexed by (5.5) define geometric simplices ∆ − (1 , , ∆ +(2 , that pass through subspaces (cid:101) H , (cid:101) H in the same order. (c) For each σ ∈ { (1 , , (2 , } , extending ∆ − σ by direct sum with subspace K yields Σ σ . In what follows, it will be convenient to use an element σ ∈ S n of the symmetric per-mutation group over n elements to specify n + 1 elements of b ∈ B n as follows: for each k ∈ { , . . . , n } , let σ ( { , . . . , k } ) ⊂ { , . . . , n } specify the unique b ∈ B n whose j -th compo-nent is +1 if and only if j ∈ σ ( { , . . . , k } ). Note that this identification yields, with someabuse of notation, σ ( { } ) = − , σ ( { , . . . , n } ) = + . Finally, note that: (cid:8) ζ σ ( { ,...,k } ) − ρ (cid:9) n − k =0 are linearly independent;(5.4a) (cid:8) ζ σ ( { ,...,k } ) + F σ ( { ,...,k } ) ( ρ ) − ρ (cid:9) nk =1 are linearly independent . (5.4b)The former fact (5.4a) is easily verified in coordinates where Dh ( ρ ) = (cid:2) I n n × ( d − n ) (cid:3) , whencethe latter fact (5.4b) follows from (5.4a) and (5.2) via [6, Cor. 5(c)] (the time- t flow of an EC r vector field is a homeomorphism of the state space for all t ∈ R ).Let ∆ denote the combinatorial simplicial complex over vertex set B n whose maximal n -simplices are indexed by σ ∈ S n via(5.5) ∆ σ = { σ ( { , . . . , k } ) } nk =0 ∈ ∆where we regard σ ( { , . . . , k } ) as an element of B n using the same abuse of notation employedin (5.4). By associating each vertex b ∈ B n with the vector ζ b ∈ Z − ⊂ R d , every n -simplex∆ σ determines an n -dimensional geometric simplex ∆ − σ ⊂ R d , the dimensionality of which isensured by (5.4a); similarly, (5.4b) ensures that associating each b ∈ B n with ( ζ b + F b ( ρ )) ∈ This manuscript is for review purposes only. Z + ⊂ R d determines an n -dimensional geometric simplex ∆ + σ ⊂ R d from each n -simplex ∆ σ .Refer to Figure 5.1(b) for an illustration when n = 2. The triple ( Z − , Z + , ∆) parameterizes acontinuous piecewise-affine homeomorphism P : | ∆ − | → | ∆ + | using the construction from [15,Sec. 3.1] (summarized in subsection 4.2), where | ∆ ± | = (cid:83) σ ∈ S n ∆ ± σ ⊂ R d denote the carriers of the geometric simplicial complexes ∆ ± .We now show that the piecewise-affine function P constructed above is the non-linear partof the time-1 flow of the sampled system (cid:101) φ restricted to | ∆ − | . For each σ ∈ S n we extend the n -dimensional geometric simplex ∆ − σ determined by the n -simplex ∆ σ via direct sum with the( d − n )-dimensional subspace K to obtain a d -dimensional polyhedron Σ σ (see Figure 5.1(c)),and let | Σ | = (cid:83) σ ∈ S n Σ σ . Note that K is a subset of the lineality space of Σ σ for each σ ∈ S n . Lemma 5.1. (cid:101) φ | | Σ | is piecewise-affine and (5.6) ∀ z ∈ (cid:12)(cid:12) ∆ − (cid:12)(cid:12) , ξ ∈ K : (cid:101) φ ( z + ξ ) = P ( z ) + ξ. Proof.
This proof will proceed in two steps: (i) show that (cid:101) φ ( z ) = P ( z ) for all z ∈ | ∆ − | ;(ii) show that (cid:101) φ ( z + ξ ) = (cid:101) φ ( z ) + ξ for all z ∈ | ∆ − | , ξ ∈ K .(i) Recall from (5.2) that (cid:101) φ | Z − = P | Z − where Z − is the vertex set for the geometricsimplicial complex ∆ − . For each σ ∈ S n let Z σ = { ζ b } b ∈ ∆ σ denote the vertex set of the n -dimensional geometric simplex ∆ − σ . Then we claim that each z ∈ ∆ − σ passes through the samesequence of transition surfaces as each ζ b ∈ Z σ . To verify this claim, we use the piecewise-affine time-to-impact function (cid:101) τ : R d → R n from subsection 4.8. Note that ζ b impacts affinesubspace (cid:101) H j at time 1 if b j = − b j = +1, i.e.(5.7) (cid:101) τ j ( ζ b ) = (cid:40) , b j = − , b j = +1 . A convex combination α ζ b + (1 − α ) ζ b (cid:48) , α ∈ (0 , b, b (cid:48) ∈ ∆ σ , impacts (cid:101) H j at time (cid:101) τ j ( α ζ b + (1 − α ) ζ b (cid:48) ) = , b j = − ∧ b (cid:48) j = − t ∈ (0 , , ( b j = +1 ∧ b (cid:48) j = − ∨ ( b j = − ∧ b (cid:48) j = +1);0 , b j = +1 ∧ b (cid:48) j = +1 . More generally, any point z ∈ ∆ − σ is a convex combination of the vertices Z σ , whence itimpacts surfaces in the order prescribed by σ :(5.8) ∀ z ∈ ∆ − σ : 0 ≤ (cid:101) τ σ (1) ( z ) ≤ (cid:101) τ σ (2) ( z ) ≤ · · · ≤ (cid:101) τ σ ( n ) ( z ) < . Thus, (cid:101) φ | ∆ − σ is affine and agrees with P | ∆ − σ . Since | ∆ − | = (cid:83) σ ∈ S n ∆ − σ , we have (cid:101) φ | | ∆ − | = P .(ii) We now show that the piecewise-affine map (cid:101) φ is indifferent to ξ ∈ K = ker Dh ( ρ ): ∀ ξ ∈ K , z ∈ (cid:12)(cid:12) ∆ − (cid:12)(cid:12) : (cid:101) φ ( z + ξ ) = (cid:101) φ ( ρ + ( z + ξ − ρ ))(5.9a) = (cid:101) φ ( ρ ) + D (cid:101) φ ( ρ ; z + ξ − ρ )(5.9b) = (cid:101) φ ( ρ ) + D (cid:101) φ ( ρ ; z − ρ ) + ξ (5.9c) = (cid:101) φ ( z ) + ξ. (5.9d) This manuscript is for review purposes only. -DERIVATIVE OF AN EC r VECTOR FIELD’S
P C r FLOW 17
Indeed: (5.9a) since z + ξ = ρ + ( z + ξ − ρ ); (5.9b) since (cid:101) φ is affine on the segment { ρ + α ( z + ξ − ρ ) : α ∈ [0 , } ; (5.9c) since each piece of the continuous piecewise-linear B-derivative D (cid:101) φ ( ρ ) is specified by a saltation matrix (as recalled in subsection 4.4) that is theproduct of matrices of the form ( I d + g · Dh j ( ρ )) [6, Eqn. (60)], thus ξ ∈ K = ker Dh ( ρ ) istransformed by I d ; (5.9d) for the same reason as (5.9b). (cid:101) φ . The goal of this subsection is to differentiate the representationof (cid:101) φ from subsection 5.1 to obtain a triangulation of the B-derivative D (cid:101) φ : T ρ − R d → T ρ + R d between the following two points:(5.10) ρ − = ρ − F − ( ρ ) , ρ + = (cid:101) φ (1 , ρ − ) = ρ + 12 F + ( ρ ) . Lemma 5.2.
The function B = D (cid:101) φ ( ρ − ) : T ρ − R d → T ρ + R d satisfies: B specifies how (cid:101) φ varies relative to (cid:101) φ ( ρ − ) , (5.11) ∀ x ∈ | Σ | : (cid:101) φ ( x ) = (cid:101) φ ( ρ − ) + B ( x − ρ − );2. B is continuous and piecewise-linear with conical subdivision (5.12) Σ (cid:48) = (cid:8) Σ (cid:48) σ = cone (cid:0) Σ σ − ρ − (cid:1) : σ ∈ S n (cid:9) ;3. B | Σ (cid:48) σ is linear for all σ ∈ S n and (5.13) ∀ δρ ∈ Σ (cid:48) σ : B ( δρ ) = M σ · δρ ;4. L = K + span F − ( ρ ) is a ( d − n + 1) -dimensional lineality space for Σ (cid:48) and (5.14) ∀ σ ∈ S n : Σ (cid:48) σ = L + cone (cid:110) Π ⊥ L · ( ζ σ ( { ,...,k } ) − ρ ) (cid:111) n − k =1 , where Π ⊥ L is the orthogonal projection onto L ⊥ ; B | L is linear and (5.15) ∀ δρ ∈ T ρ − R d : B ( δρ ) = B (Π L · δρ ) + B (cid:16) Π ⊥ L · δρ (cid:17) , where Π L is the orthogonal projection onto L .Proof. Each point follows from straightforward application of results in [36]: (1.), (2.),and (3.) are conclusions (4.), (3.), and (2.), respectively, of [36, Prop. 2.2.6]; (4.) followsfrom the definitions of lineality space [36, Sec. 2.1.2] and the ζ b ’s (5.1); (5.) is a restatementof [36, Lem. 2.3.2]. φ . The goal of this subsection is to show that the piecewise-linearfunction B triangulated in subsection 5.2 gives the non-linear part of the desired B-derivative D x φ ( t, x ) and (5.16) W σ = D x φ ( s, x ) − (cid:0) Σ (cid:48) σ (cid:1) ⊂ T x R d is the cone of tangent vectors where the saltation matrix M σ is active in (4.8). Here and in what follows we mildly abuse notation via the natural vector space isomorphisms R d (cid:39) T ρ − R d (cid:39) T ρ + R d (cid:39) T ρ R d . This manuscript is for review purposes only.
Theorem 5.3.
Suppose the vector field F : R d → T R d is event-selected C r with respect to h : R d → R n at ρ . Let φ : F → R d be the P C r flow of F and s, t ∈ R , x ∈ R d be such that < s < t and F is C r on φ ([0 , t ] \ { s } , x ) ⊂ R d . Then with ρ = φ ( s, x ) , the B-derivative ofthe flow φ with respect to state, D x φ ( t, x ) : T x R d → T φ ( t,x ) R d , is given by ∀ δx ∈ T x R d : D x φ ( t, x ; δx ) = D x φ ( t − s, ρ ) · B ( D x φ ( s, x ) · δx ) , (5.17a) ∀ δx ∈ W σ ⊂ T x R d : D x φ ( t, x ; δx ) = D x φ ( t − s, ρ ) · M σ · D x φ ( s, x ) · δx, (5.17b) where B is the continuous piecewise-linear function from Lemma 5.2, W σ is the cone from (5.16) ,and M σ is the saltation matrix from (4.9) .Proof. Note that (5.17a) follows from (5.17b) by (5.13), and the fact that “pieces” of theB-derivative D x φ ( t, x ) are determined by the collection of saltation matrices { M σ } σ ∈ S n wasrecalled in subsection 4.4. Thus, to establish (5.17b) what remains to be shown is that M σ is the active “piece” for all δx ∈ W σ , i.e. that { W σ } σ ∈ S n is a conical subdivision for thepiecewise-linear operator D x φ ( t, x ), with W σ as defined in (5.16).Given δx ∈ Int W σ let δρ = D x φ ( s, x ) · δx ∈ Int Σ (cid:48) σ so that(5.18) (cid:101) τ σ (1) ( ρ + δρ ) < (cid:101) τ σ (2) ( ρ + δρ ) < · · · < (cid:101) τ σ ( n ) ( ρ + δρ )where (cid:101) τ is the time-to-impact function for the sampled system as defined in (4.14). Note that D x φ ( t, x ) is linear on span F ( x ),(5.19) ∀ α ∈ R : D x φ ( t, x ; δx + αF ( x )) = D x φ ( t, x ; δx ) + αF ( φ ( t, x )) , so without loss of generality we may assume δρ ∈ Int (cid:101) D − by translating δx in the − F ( x )direction. We claim that, for all α > φ ( t, x + α δx ) passes through theevent surfaces with the same sequence as (cid:101) φ (1 , ρ + α δρ ), i.e. that(5.20) τ σ (1) ( x + α δx ) < τ σ (2) ( x + α δx ) < · · · < τ σ ( n ) ( x + α δx ) , where τ is the time-to-impact function defined in (4.13). To see this, note that ∀ k ∈ { , . . . , n } : τ σ ( k ) ( x + α δx ) − τ σ ( k ) ( x ) = Dτ σ ( k ) ( x ; α δx ) + O (cid:0) α (cid:1) (5.21a) = D (cid:101) τ σ ( k ) ( ρ ; α δρ ) + O (cid:0) α (cid:1) (5.21b) = (cid:101) τ σ ( k ) ( ρ + α δρ ) − (cid:101) τ σ ( k ) ( ρ ) + O (cid:0) α (cid:1) (5.21c)where: (5.21a) since τ is P C r ; (5.21b) since δρ = D x φ ( s, x ) · δx and Dτ ( x ; δx ), D (cid:101) τ ( ρ ; δρ )are are determined by the same data, namely, Dh σ ( k ) ( ρ ) and F − ( ρ ); (5.21c) since δρ ∈ Σ (cid:48) σ .Combining the approximation (5.21) with (5.18) yields (5.20) as desired.We conclude that { W σ } σ ∈ S n is a conical subdivision for the piecewise-linear operator D x φ ( t, x ), which verifies (5.17) and completes the proof. Remark B . Although there are n ! pieces of B in general, we explicitly repre-sent all pieces using a triangulation of 2 n sample points defined in (5.3), achieving a substantialreduction – from factorial to “merely” exponential – of the information needed to representthe first-order approximation of the flow. Note that B implicitly determines the transitionsequence σ associated with the perturbation direction δx in (5.17a), whereas this sequencemust be explicitly specified to select the appropriate saltation matrix M σ in (5.17b). This manuscript is for review purposes only. -DERIVATIVE OF AN EC r VECTOR FIELD’S
P C r FLOW 19
6. Computation.
We now attend to the complexity of the computational tasks requiredto construct or evaluate the B-derivative representation from the preceding section. To thatend, let F : R d → T R d be an event-selected C r vector field with respect to h ∈ C r ( R d , R n )and φ : F → R d its piecewise- C r flow, and assume s, t ∈ R and x ∈ R d are such that 0 < s < t , ρ = φ ( s, x ), and the vector field F is C r on φ ([0 , t ] \ { s } , x ).We seek to compute D x φ ( t, x ; δx ) given δx ∈ T x R d . Since (5.17a) from Theorem 5.3 yields(6.1) D x φ ( t, x ; δx ) = D x φ ( t − s, x ) · B ( D x φ ( s, x ) · δx )where B : T ρ R d → T ρ R d , the crux of the computation is(6.2) δρ + = B ( δρ − )where δρ − = D x φ ( s, x ) · δx . In fact, Lemma 5.2 offers further simplification via (5.15): since B = B ◦ Π L + B ◦ Π ⊥ L where B ◦ Π L is the linear function(6.3) B ◦ Π L · δρ − = (cid:18) I d + ( F + ( ρ ) − F − ( ρ )) · F − ( ρ ) (cid:62) (cid:107) F − ( ρ ) (cid:107) (cid:19) · Π L · δρ − , only the piecewise-linear function B ◦ Π ⊥ L (equivalently, the restriction B | L ⊥ ) requires specialconsideration. In what follows, we will assume the following data, needed to construct the sampled system illustrated in Figure 2.1(b), is given: linearly-independent normal vectors forthe surfaces of discontinuity, i.e. Dh ( ρ ) ∈ R n × d with rank Dh ( ρ ) = n ; limiting values of thevector field at the point of intersection, i.e. F b ( ρ ) ∈ T ρ R d for each b ∈ B n ; and F-derivativesof the continuously-differentiable parts of the flow, i.e. D x φ ( s, x ) , D x φ ( t − s, x ) ∈ R d × d . Lemma 5.2 demonstrates that there are n ! pieces ofthe piecewise-linear function B , namely, the collection of saltation matrices { M σ } σ ∈ S n in (5.13)that are active on the corresponding polyhedral cones in the conical subdivision Σ (cid:48) = { Σ (cid:48) σ } σ ∈ S n in (5.12). These polyhedral cones are generated by the 2 n − points { ζ b : b ∈ B n \ {− , + }} in (5.14). For each b ∈ B n , the point ζ b ∈ K ⊥ + { ρ } where K = ker Dh ( ρ ) can be determinedby solving the n affine equations with n unknowns in (5.1). Given σ ∈ S n , the linear piece B | L ⊥ ∩ Σ (cid:48) σ can be constructed using the saltation matrix [6, Sec. 7.1.6] since B ( δρ − ) = M σ · δρ − for all δρ − ∈ L ⊥ ∩ Σ (cid:48) σ where (6.4) M σ = n − (cid:89) k =0 (cid:32) I d + (cid:0) F σ ( { ,...,k +1 } ) ( ρ ) − F σ ( { ,...,k } ) ( ρ ) (cid:1) Dh σ ( { ,...,k } ) ( ρ ) · F σ ( { ,...,k } ) ( ρ ) · Dh σ ( { ,...,k } ) ( ρ ) (cid:33) , or using barycentric coordinates [15, Eqn. (3.1)] since B ( δρ − ) = Z + σ · ( Z − σ ) † · δρ − for all δρ − ∈ L ⊥ ∩ Σ (cid:48) σ where(6.5) Z ± σ = (cid:104) z ± σ ( { , } ) z ± σ ( { , , } ) · · · z ± σ ( { , ,...,n − } ) (cid:105) ∈ R d × ( n − , We mildly abuse notation as in subsection 5.1 by using σ ∈ S n to specify n + 1 elements of b ∈ B n : foreach k ∈ { , . . . , n } , we let σ ( { , . . . , k } ) ⊂ { , . . . , n } specify the unique b ∈ B n whose j -th component is +1if and only if j ∈ σ ( { , . . . , k } ). This manuscript is for review purposes only. (6.6) ∀ b ∈ ∆ (cid:48) σ : z − b = Π ⊥ L · ( ζ b − ρ ) , z + b = B | L ⊥ ( z − b ) , (6.7) ∆ (cid:48) σ = { σ ( { , , . . . , k } ) } n − k =1 ;note that the pseudo-inverse ( Z − σ ) † is injective on L ⊥ ∩ Σ (cid:48) σ by (5.4a) and (5.14). Althoughthe matrices M σ , Z + σ · ( Z − σ ) † ∈ R d × d define the same linear transformation on the ( n − L ⊥ ∩ Σ (cid:48) σ , they are generally not the same matrix. We conclude by noting thatconstructing the saltation matrix in (6.4) requires O (cid:0) nd (cid:1) time and O (cid:0) d (cid:1) space, whereasconstructing the Barycentric coordinates in (6.5) requires O (cid:0) n d (cid:1) time and O (cid:0) d (cid:1) space(although evaluating the expression Z + σ · ( Z − σ ) † · δρ − requires only O (cid:0) nd (cid:1) time given Z ± σ ). One obvious strategy to evaluate B on δρ − ∈ T ρ R d isto (i) determine σ ∈ S n such that δρ − ∈ Σ (cid:48) σ then (ii) apply the corresponding saltation matrixor barycentric coordinates calculation from the preceding section. The general formulation of(i), termed the point location problem in the computational geometry literature, is “essentiallyopen” [9, Sec. 6.5]. For an arrangement of m hyperplanes in R d , queries can be answered in O ( d log m ) time at the expense of O (cid:0) m d (cid:1) space [7]. In our context, the conical subdivision Σ (cid:48) in (5.14) is determined by an arrangement of m = O (cid:0) n ! (cid:1) hyperplanes, so this general-purposealgorithm has time complexity O ( d log n !) = O ( d n log n ) and space complexity O (cid:0) n ! d (cid:1) .The relationship established by (5.11) between the desired B-derivative and the flow ofthe sampled system illustrated in Figure 2.1(b) suggests a different strategy, summarizedin Figure 2.2, with slightly worse O (cid:0) n d (cid:1) time complexity but dramatically superior O ( d )space complexity. To understand the strategy, interpret the tangent vector δρ − ∈ T ρ − R d as aperturbation away from the point ρ − = ρ − F − ( ρ ) that flows through ρ to ρ + = ρ + F + ( ρ )in one unit of time and observe that δρ + = (cid:101) φ ( ρ − + δρ − ) − ρ + = B ( δρ − ) as in (5.11). Theflow of the sampled system (cid:101) φ is piecewise-affine, and can be evaluated on a given perturbationvector δρ − by performing a sequence of n affine projections (one for each of the affine subspaces (cid:110) (cid:101) H j (cid:111) nj =1 where (cid:101) F is discontinuous) specified by the permutation σ ∈ S n for which δρ − ∈ Σ (cid:48) σ .Fortuitously, the sequence σ can be determined inductively as follows. First, define(6.8) δt = 0 ,δρ = δρ − ,σ (1) = arg min (cid:26) − Dh j ( ρ ) · δρ Dh j ( ρ ) · F − ( ρ ) : j ∈ { , . . . , n } (cid:27) ,τ = − Dh σ (1) ( ρ ) · δρ Dh σ (1) ( ρ ) · F − ( ρ ) . This equation only holds when (cid:13)(cid:13) δρ − (cid:13)(cid:13) is small enough to ensure ρ − + δρ − ∈ (cid:101) D − and ρ + + δρ + ∈ (cid:101) D + ;since the B-derivative is positively-homogeneous, we impose this restriction without loss of generality. This manuscript is for review purposes only. -DERIVATIVE OF AN EC r VECTOR FIELD’S
P C r FLOW 21
Then for k ∈ { , . . . , n − } inductively define(6.9) δt k +1 = δt k + τ k ,δρ k +1 = δρ k + τ k · F σ ( { ,...,k − } ) ( ρ ) ,σ ( k + 1) = arg min (cid:26) − Dh j ( ρ ) · δρ k +1 Dh j ( ρ ) · F σ ( { ,...,k } ) ( ρ ) : j ∈ { , . . . , n } \ σ ( { , . . . , k } ) (cid:27) ,τ k +1 = − Dh σ ( k +1) ( ρ ) · δρ k +1 Dh σ ( k +1) ( ρ ) · F σ ( { ,...,k } ) ( ρ ) . Finally, set δρ + = δρ n − ( δt n + τ n ) · F + ( ρ ). By construction, δρ − ∈ Σ (cid:48) σ and δρ + = B ( δρ − ).This strategy is succinctly summarized in pseudocode and sourcecode in Figure 2.2; its timecomplexity is O (cid:0) n d (cid:1) since there are n steps in the induction and each step requires O ( n ) dotproducts between d -vectors. The space complexity is O ( d ) since each step in the inductionrequires O ( d ) storage and data from preceding steps can be forgotten or overwritten.We conclude by noting that, if a general-purpose algorithm is employed to solve the pointlocation problem in O ( d n log n ) time to obtain the sequence σ ∈ S n , then the inductiondescribed in the preceding paragraph can be simplified by skipping the steps that determine σ (1) and σ ( k + 1) from (6.8) and (6.9). This simplification reduces the time complexity of theinduction to O ( nd ), so the overall algorithm retains the O ( d n log n ) time complexity of thegeneral-purpose point-location algorithm (at the expense of the superexponential O (cid:0) n ! d (cid:1) spacecomplexity of the point location algorithm). We are pessimistic these asymptotic complexitiescan be improved in general.
7. Conclusion.
We constructed a representation for the
Bouligand (or B- )derivative of the piecewise- C r ( P C r ) flow generated by an event-selected C r ( EC r ) vector field and appliedthe representation to derive a polynomial-time algorithm to evaluate the B-derivative ona given tangent vector. Our results provide a foundation that may support future workgeneralizing classical analysis and synthesis techniques for smooth control systems to theclass of nonsmooth systems considered here. In particular, we envision applying our resultsto design and control the class of mechanical systems subject to unilateral constraints thatarise in models of robot locomotion and manipulation. REFERENCES [1]
J. Aguilar and D. I. Goldman , Robophysical study of jumping dynamics on granular media , Naturephysics, 12 (2015), p. nphys3568, https://doi.org/10.1038/nphys3568.[2]
M. A. Aizerman and F. R. Gantmacher , Determination of stability by linear approximation of a peri-odic solution of a system of differential equations with discontinuous right-hand sides , The QuarterlyJournal of Mechanics and Applied Mathematics, 11 (1958), pp. 385–398, https://doi.org/10.1093/qjmam/11.4.385.[3]
P. Ballard , The dynamics of discrete mechanical systems with perfect unilateral constraints , Archive forRational Mechanics and Analysis, 154 (2000), pp. 199–274, https://doi.org/10.1007/s002050000105.[4]
D. P. Bertsekas , Nonlinear Programming , Athena Scientific, 2nd ed., 1999.[5]
F. Bizzarri, A. Brambilla, and G. Storti Gajani , Lyapunov exponents computation for hybridneurons , Journal of Computational Neuroscience, 35 (2013), pp. 201–212, https://doi.org/10.1007/s10827-013-0448-6.
This manuscript is for review purposes only. [6]
S. A. Burden, S. S. Sastry, D. E. Koditschek, and S. Revzen , Event-selected vector field disconti-nuities yield piecewise-differentiable flows , SIAM Journal on Applied Dynamical Systems, 15 (2016),pp. 1227–1267, https://doi.org/10.1137/15M1016588.[7]
B. Chazelle and J. Friedman , Point location among hyperplanes and unidirectional ray-shooting ,Computational Geometry, 4 (1994), pp. 53–62, https://doi.org/10.1016/0925-7721(94)90009-4.[8]
S. H. Collins, M. B. Wiggin, and G. S. Sawicki , Reducing the energy cost of human walking usingan unpowered exoskeleton , Nature, 522 (2015), pp. 212–215, https://doi.org/10.1038/nature14288.[9]
M. de Berg, M. van Kreveld, M. Overmars, and O. Schwarzkopf , Computational Geometry:Algorithms and Applications , Springer, 2000.[10]
M. di Bernardo, C. Budd, A. R. Champneys, and P. Kowalczyk , Piecewise-smooth dynamicalsystems: theory and applications , vol. 163, Springer Science + Business Media, 2008.[11]
L. Dieci and L. Lopez , Fundamental matrix solutions of piecewise smooth differential systems , Math-ematics and Computers in Simulation, 81 (2011), pp. 932–953, https://doi.org/10.1016/j.matcom.2010.10.012.[12]
R. Elandt, E. Drumwright, M. Sherman, and A. Ruina , A pressure field model for fast, ro-bust approximation of net contact force and moment between nominally rigid objects , in IEEE/RSJInternational Conference on Intelligent Robots and Systems (IROS), Nov. 2019, pp. 8238–8245,https://doi.org/10.1109/IROS40897.2019.8968548.[13]
J. Eldering and H. Jacobs , The role of symmetry and dissipation in biolocomotion , SIAM Journal onApplied Dynamical Systems, 15 (2016), pp. 24–59, https://doi.org/10.1137/140970914.[14]
A. F. Filippov , Differential equations with discontinuous righthand sides , Springer, 1988.[15]
R. E. Groff , Piecewise linear homeomorphisms for approximation of invertible maps , PhD thesis, Uni-versity of Michigan, 2003.[16]
A. Hatcher , Algebraic topology , Cambridge University Press, 2002.[17]
J. P. Hespanha , Linear systems theory , Princeton University Press, 2009.[18]
I. A. Hiskens and M. A. Pai , Trajectory sensitivity analysis of hybrid systems , IEEE Transactionson Circuits and Systems I: Fundamental Theory and Applications, 47 (2000), pp. 204–220, https://doi.org/10.1109/81.828574.[19]
A. P. Ivanov , The stability of periodic solutions of discontinuous systems that intersect several surfacesof discontinuity , Journal of Applied Mathematics and Mechanics, 62 (1998), pp. 677–685, https://doi.org/10.1016/S0021-8928(98)00087-2.[20]
M. R. Jeffrey , Dynamics at a switching intersection: hierarchy, isonomy, and multiple sliding ,SIAM Journal on Applied Dynamical Systems, 13 (2014), pp. 1082–1105, https://doi.org/10.1137/13093368X.[21]
A. M. Johnson, S. A. Burden, and D. E. Koditschek , A hybrid systems model for simple manip-ulation and self-manipulation systems , The International Journal of Robotics Research, 35 (2016),pp. 1354–1392, https://doi.org/10.1177/0278364916639380.[22]
R. I. Leine and H. Nijmeijer , Dynamics and bifurcations of non-smooth mechanical systems , vol. 18,Springer Science + Business Media, 2013.[23]
L. Ljung , System identification: theory for the user , Prentice-Hall, 1999.[24]
P. L¨otstedt , Mechanical systems of rigid bodies subject to unilateral constraints , SIAM Journal onApplied Mathematics, 42 (1982), pp. 281–296, https://doi.org/10.1137/0142022.[25]
T. E. Oliphant , A guide to NumPy , vol. 1, Trelgol Publishing USA, 2006.[26]
A. M. Pace and S. A. Burden , Piecewise-differentiable trajectory outcomes in mechanical systems sub-ject to unilateral constraints , in Hybrid Systems: Computation and Control (HSCC), 2017, pp. 243–252, https://doi.org/10.1145/3049797.3049807.[27]
T. S. Parker and L. O. Chua , Practical numerical algorithms for chaotic systems , Springer, 1989.[28]
E. Polak , Optimization: algorithms and consistent approximations , Springer-Verlag, 1997.[29]
L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko , The mathe-matical theory of optimal processes (translated by KN Trirogoff) , John Wiley & Sons, 1962.[30]
Python Software Foundation , Python language reference, version 3.7 , https://docs.python.org/release/3.7.0/.[31]
D. Ralph and S. Scholtes , Sensitivity analysis of composite piecewise smooth equations , MathematicalProgramming, 76 (1997), pp. 593–612, https://doi.org/10.1007/BF02614400.
This manuscript is for review purposes only. -DERIVATIVE OF AN EC r VECTOR FIELD’S
P C r FLOW 23 [32]
S. M. Robinson , Local structure of feasible sets in nonlinear programming, part III: stability and sensitiv-ity , Nonlinear Analysis and Optimization, 30 (1987), pp. 45–66, https://doi.org/10.1007/BFb0121154.[33]
R. T. Rockafellar , A property of piecewise smooth functions , Computational Optimization and Appli-cations, 25 (2003), pp. 247–250, https://doi.org/10.1023/A:1022921624832.[34]
S. Sastry and M. Bodson , Adaptive control: stability, convergence, and robustness , Prentice Hall, 1989.[35]
S. S. Sastry , Nonlinear Systems: Analysis, Stability, and Control , Springer, 1999.[36]
S. Scholtes , Introduction to piecewise differentiable equations , Springer-Verlag, 2012, https://doi.org/10.1007/978-1-4614-4340-7.[37]
S. N. Simic, K. H. Johansson, J. Lygeros, and S. Sastry , Towards a geometric theory of hybridsystems
K. Sreenath, H.-W. Park, I. Poulakakis, and J. W. Grizzle , A Compliant Hybrid Zero DynamicsController for Stable, Efficient and Fast Bipedal Walking on MABEL , The International Journal ofRobotics Research, 30 (2011), pp. 1170–1193, https://doi.org/10.1177/0278364910379882.[39]
A. Tornambe , Modeling and control of impact in mechanical systems: theory and experimental results ,IEEE Transactions on Automatic Control, 44 (1999), pp. 294–309, https://doi.org/10.1109/9.746255.[40]
V. Utkin , Variable structure systems with sliding modes , IEEE Transactions on Automatic Control, 22(1977), pp. 212–222, https://doi.org/10.1109/TAC.1977.1101446., IEEE Transactions on Automatic Control, 22(1977), pp. 212–222, https://doi.org/10.1109/TAC.1977.1101446.