Vaibhav Donde
email: vdonde@lbl.gov

Ian A Hiskens
email: hiskens@engr.wisc.edu

SHOOTING METHODS FOR LOCATING GRAZING PHENOMENA IN HYBRID SYSTEMS *

Keywords: Grazing, limit cycles, piecewise smooth dynamical systems, shooting methods, continuation methods

Hybrid systems are typified by strong coupling between continuous dynamics and discrete events. For such piecewise smooth systems, event triggering generally has a significant influence over subsequent system behavior. Therefore, it is important to identify situations where a small change in parameter values alters the event triggering pattern. The bounding case, which separates regions of (generally) quite different dynamic behaviors, is referred to as grazing. At a grazing point, the system trajectory makes tangential contact with an event triggering hypersurface. The paper formulates conditions governing grazing points. Both transient and periodic behaviors are considered. The resulting boundary value problems are solved using shooting methods that are applicable for general nonlinear hybrid (piecewise smooth) dynamical systems. The grazing point formulation underlies the development of a continuation process for exploring parametric dependence. It also provides the basis for an optimization technique that finds the smallest parameter change necessary to induce grazing. Examples are drawn from power electronics, power systems and robotics, all of which involve intrinsic interactions between continuous dynamics and discrete events.

Introduction

Many systems exhibit dynamic behavior that is best characterized by strong coupling between continuous dynamics and discrete events. Trajectories evolve smoothly through state-space until satisfying conditions that trigger an event. The event alters the system description and/or induces an impulsive change in the state. Smooth evolution then continues until the next event is triggered. Such systems have come to be known as hybrid systems [START_REF] Van Der Schaft | An Introduction to Hybrid Dynamical Systems[END_REF][START_REF] Liberzon | Switching in Systems and Control[END_REF] or piecewise smooth dynamical systems [START_REF] Di Bernardo | Guest editorial: Special issue on switching and systems[END_REF]. 1 Examples can be drawn from a wide range of application areas, including process control [START_REF] Lennartson | Hybrid systems in process control[END_REF], constrained mechanical systems [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control[END_REF], robotics [START_REF] Spong | Robot Dynamics and Control[END_REF][START_REF] Piiroinen | Recurrent dynamics of nonsmooth systems with application to human gait[END_REF], power systems [START_REF] Hiskens | Power system modeling for inverse problems[END_REF], and power electronics [START_REF] Rajaraman | Nonlinear dynamics and switching time bifurcations of a thyristor controlled reactor circuit[END_REF][START_REF] Yuan | Border-collision bifurcations in the buck converter[END_REF]]. In fact, any physical device that exhibits hysteresis, or control loop with anti-wind-up limits [START_REF] Goodwin | Control System Design[END_REF], is effectively a hybrid system.

Event triggering generally has a significant influence over subsequent system behavior. Therefore, it is important to identify situations where a small change in parameter values alters the event triggering pattern. Figure 1 provides an illustration. For a certain value of parameter θ hit , the system trajectory encounters an event triggering hypersurface at a point x hit . The event occurs, and the trajectory continues accordingly. However for a small change in parameter value to θ miss , the trajectory misses (at least locally) the triggering hypersurface, and subsequently exhibits a completely different form of response. At a parameter value θ g , lying between θ hit and θ miss , the continuous trajectory tangentially encounters (grazes) the triggering hypersurface. Behavior beyond the grazing point x g is generally unpredictable, in the sense that without further knowledge of the system, it is impossible to determine whether or not the event triggers. This bounding case typifies grazing phenomena. The paper proposes a shooting method for locating grazing points x g and the associated critical parameter value θ g .

The significance of grazing can be illustrated by considering power system protection. For a

Parameter values

Triggering hypersurface particular set of parameters, a disturbance may cause the system trajectory to pass close to (but not encounter) a protection triggering characteristic [START_REF] Singh | Direct assessment of protection operation and non-viable transients[END_REF]. The protection does not operate, no equipment is tripped, and the system recovers. However it may require only a small change in parameter values, such as loads, for protection to operate, perhaps leading to cascading outages and system fragmentation [U. S.-Canada Power System Outage Task Force, 2004]. Knowledge of the grazing point would highlight system vulnerability.

In a more general sense, grazing plays a fundamental role in dynamic embedded optimization of hybrid systems. Such optimization problems incorporate the system dynamic description into the constraint set [START_REF] Galán | Dynamic optimization of hybrid systems[END_REF][START_REF] Hiskens | Power system modeling for inverse problems[END_REF]. As [START_REF] Galán | Dynamic optimization of hybrid systems[END_REF] and [START_REF] Piccoli | Hybrid systems and optimal control[END_REF] have shown, even though hybrid system dynamics are nonsmooth, the cost function is smooth away from grazing points. This statement is expressed using different terminology though, with [START_REF] Piccoli | Hybrid systems and optimal control[END_REF] requiring that all flows have the same history, and [START_REF] Galán | Dynamic optimization of hybrid systems[END_REF] requiring the existence of trajectory sensitivities (discussed later). Therefore, away from grazing points, gradient-based optimization techniques are appropriate. However, smoothness and even continuity of the cost function is often lost at grazing points. This introduces the need for branch-and-bound algorithms that search over the grazing-induced partitions of parameter space.

Grazing phenomena have been widely investigated recently, particularly with reference to periodic systems [START_REF] Dankowicz | Low-velocity impacts of quasiperiodic oscillations[END_REF][START_REF] Chin | Universal behavior of impact oscillators near grazing incidence[END_REF][START_REF] Fredriksson | Bifurcations caused by grazing incidence in many degrees of freedom impact oscillators[END_REF]. In that context, grazing is closely related to bordercollision bifurcations [START_REF] Di Bernardo | Grazing and border-collision in piecewise-smooth systems: A unified analytical framework[END_REF][START_REF] Nusse | Bordercollision bifurcations: An explanation for observed bifurcation phenomena[END_REF][START_REF] Yuan | Border-collision bifurcations in the buck converter[END_REF], also known as C-bifurcations [START_REF] Di Bernardo | Local analysis of C-bifurcations in ndimensional piecewise-smooth dynamical systems[END_REF]. Transient grazing is considered in [START_REF] Rajaraman | Nonlinear dynamics and switching time bifurcations of a thyristor controlled reactor circuit[END_REF][START_REF] Jalali | Switching time bifurcations in a thyristor controlled reactor[END_REF], where switching-time bifurcations are analyzed. These investigations have focused largely on classifying the (local) consequences of grazing through normal form analysis, particularly for periodic systems. Computation of actual grazing points has received less attention. With numerical packages such as AUTO [START_REF] Doedel | Auto97: Continuation and bifurcation software for ordinary differential equations[END_REF]] generally unsuited to nonsmooth systems, 2 ad hoc approaches have prevailed. This paper addresses that deficiency by establishing a shooting method that is applicable for general nonlinear hybrid (piecewise smooth) dynamical systems.

The paper is organized as follows. Section 2 considers the modeling of hybrid systems. Section 3 provides background material for trajectory sensitivities and limit cycles. Conditions governing transient grazing are developed in Sec. 4, and adapted to periodic systems in Sec. 5. The grazing point formulation is extended to a continuation process in Sec. 6. Section 7 develops an optimization technique for finding the smallest parameter change necessary to induce grazing. Conclusions are provided in Sec. 8.

Hybrid System Representation

Model

As mentioned previously, hybrid system behavior is characterized by continuous and discrete states, continuous dynamics, discrete events (triggers), and mappings that define the evolution of discrete states at events. Numerous formal models exist for rigorously describing hybrid system dynamics. Examples include petri nets [START_REF] David | Petri Nets and Grafecet[END_REF] and hybrid automata [START_REF] Van Der Schaft | An Introduction to Hybrid Dynamical Systems[END_REF]. However, those representations are not immediately amenable to numerical implementation. A useful, nonrestrictive model formulation should be,

• capable of capturing the full range of continuous/ discrete hybrid system dynamics, • computationally efficient, and • consistent with the development of shooting methods.3

As shown in [START_REF] Hiskens | Trajectory sensitivity analysis of hybrid systems[END_REF][START_REF] Hiskens | Power system modeling for inverse problems[END_REF], these specifications can be met completely by a model that consists of a set of differential-algebraic equations, adapted to incorporate impulsive (state reset) action and switching of the algebraic equations. This DA Impulsive Switched (DAIS) model can be written in the form,

ẋ = f (x, y) + mr j=1 δ(y r[j] )(h j (x, y) -x) (1) 0 = g(x, y) ≡ g (0) (x, y) + ms i=1 g (i) (x, y) (2)
where

g (i) (x, y) = g (i-) (x, y) g (i+) (x, y) y s[i] < 0 y s[i] > 0 i = 1, . . . , m s
(3) and • x ∈ R n are dynamic states, and y ∈ R m are algebraic states; • δ(•) is the Dirac delta. Each impulse term of the summation in (1) can be expressed in the alternative state reset form

x + = h j (x -, y -) when y r[j] = 0 (4)
where the notation x + denotes the value of x just after the reset event, while x -and y -refer to the values of x and y just prior to the event. This form motivates a generalization to an implicit mapping

h j (x + , x -, y -) = 0. • Subscripts r[j], s[i] index elements of y that trig-
ger the jth state reset (impulsive) event and ith algebraic switching event, respectively. (This refers to the ith modeled event, not the event that occurs ith in the time sequence.) • f, h j : R n+m → R n .

• g (0) , g (i±) : R n+m → R m . Some elements of each g (•) will usually be identically zero, but no elements of the composite g should be identically zero. Each g (i±) may itself have a switched form, and is defined similarly to (2)-( 3), leading to a nested structure for g.

A nonautonomous version of the DAIS model follows from explicit inclusion of time t in f , g, and h. Alternatively time can be modeled as a state, with dynamics ẋt = 1 and initial condition x t (0) = t 0 . This monotonic state x t may complicate analysis of periodic systems. If so, periodic x t can be achieved using a reset equation of the form (4),

x + t = x - t -T when x t = t 0 + T,
where T is the period. A compact development of the equations describing grazing phenomena results from incorporating parameters p ∈ R into the dynamic states x. (Numerical implementation is also simplified.) This is achieved by introducing trivial differential equations ṗ = 0

(5) into (1), and results in the natural partitioning

x = x p , f = f 0 , h j = h j p (6)
where x are the true dynamic states, and p are parameters.

Away from events, system dynamics evolve smoothly according to the familiar differentialalgebraic model

ẋ = f (x, y) ( 7) 0 = g(x, y) ( 8)
where g is composed of g (0) together with appropriate choices of g (i-) or g (i+) , depending on the signs of the corresponding elements of y s [i] . At switching events (3), some component equations of g change.

To satisfy the new g = 0 constraints, algebraic variables y may undergo a step change. Impulse events (1) (equivalently reset events (4)) force a discrete change in elements of x. Algebraic variables may again step to ensure g = 0 is always satisfied.

The flows of x and y are defined as

x(t) = φ(x 0 , t) = φ(x 0 , t) p (9) y(t) = ψ(x 0 , t) (10)
where x(t) and y(t) satisfy ( 1)-( 3), along with initial conditions,

φ(x 0 , t 0 ) = x 0 (11) g(x 0 , ψ(x 0 , t 0 )) = 0. ( 12
)
The partitioning of φ in ( 9) is in accordance with (6).

Example 1 (Modeling)

The simple static var compensator (SVC) of Fig. 2 provides a practical illustration of the DAIS model formulation. This example will be revisited later to examine various algorithms for locating grazing phenomena. The parameters of this model are taken from [START_REF] Jalali | Switching time bifurcations in a thyristor controlled reactor[END_REF]. For clarity of presentation we consider only a single thyristor, whereas [START_REF] Jalali | Switching time bifurcations in a thyristor controlled reactor[END_REF] used a back-to-back pair.

Resistor R s and inductor L s characterize the transformer and source impedance, while R r and L r are the resistance and inductance of the inductor coil in series with the thyristor switch. The thyristor conducts current i Lr only in the forward direction. It turns on when triggered by a firing pulse at firing angle α, and continues to conduct while i Lr remains positive. It turns off (commutates) when i Lr becomes zero. The thyristor is assumed ideal, so detailed nonlinearities in the on/off process are neglected. The model consists of differential equations

R r ,L

r R s ,L s C (a) Firing angle, α + - + - L s =0.195 mH R s =0.9 mΩ R r =31.3 mΩ L r =1.66 mH v Lr v Ls i Ls i Lr v c C=1.5 mF V s = sin ωt (b)
d dt           x t i Ls i Lr v c x rr x stat           =                   1 v Ls L s v Lr L r i c C 0 0                   (13)
together with algebraic equations

0 = R s i Ls + v Ls + v c -sin(ωx t ) (14) 0 = i Ls -i c -i Lr ( 15
)
y on = ωx t π -x rr - α 180 (16) v Lr = 0 y off = -1 x stat < 0 (thyristor open circuit) (17) v c -R r i Lr -v Lr = 0 y off = i Lr x stat > 0 (thyristor conducting) (18)
and state reset equations

x + rr = x - rr + 2 x + stat = x - stat + 2 when y on = 0 (19) x + stat = x - stat -2 when y off = 0. ( 20
)
Model behavior can be best explained by referring to Fig. 3. Assume that initial time x t = 0, and the thyristor is in blocking (open circuit) mode, that is x stat = -1. Initially x rr = 0, so according to ( 16), the turn-on indicator y on < 0. Subsequently y on increases linearly with time until a reset event ( 19) is triggered when y on = 0. At that event, the thyristor status x stat toggles from -1 to +1, i.e. from blocking to conducting. Also, a step occurs in x rr , forcing y on < 0 again. While the thyristor is conducting, the turn-off indicator y off monitors the thyristor current i Lr . If i Lr drops to zero, reset event (20) toggles x stat from +1 to -1, and the thyristor self commutates (blocks).

In order to motivate later grazing analysis, Fig. 4 shows thyristor current i Lr for various values of firing angle. It can be seen that as the firing angle increases from α = 120 • to 125 • , the dip during the first conduction period distorts towards zero. Somewhere between 125 • and 130 • , the increasing distortion causes the current to fall to zero, and the thyristor prematurely commutates. The pivotal case, separating the two qualitatively different forms of behavior, corresponds to grazing. As the firing angle increases further, another grazing situation occurs between 130 • and 135 • , with premature commutation occurring during the second conduction period.

Numerical integration

Predicting the dynamic response of a system, for a given set of initial conditions, generally requires numerical integration. The accuracy of the resulting approximate trajectory depends on factors such as integration technique, time step and roundoff error [START_REF] Parker | Practical Numerical Algorithms for Chaotic Systems[END_REF][START_REF] Shampine | Numerical Solution of Ordinary Differential Equations[END_REF]. Gross errors, including numerical instability, can be avoided by using implicit integration with an appropriate time step. However, some inaccuracy will always remain. In cases where event triggering hypersurfaces are encountered transversally, slight inaccuracy in the predicted trajectory will have only a minor effect on the triggering point. 4 However in grazing situations, where the encounter is tangential, behavior will be as sensitive to simulation parameters as it is to true system parameters.

To illustrate, Fig. 5 shows the first cycle of SVC inductor current i Lr , for a firing angle α = 125 • . Trapezoidal integration was used to obtain this trajectory; the responses predicted using two different time steps are shown. The two trajectories are almost identical, except for slight divergence around the peaks and troughs. However referring to Fig. 4, it is the trough that undergoes grazing as the firing angle α is increased. The predicted grazing value of α will therefore depend (slightly) on the integration time step.

There are certainly integration techniques that are more accurate than trapezoidal, so numerical inaccuracies can be further reduced. However the tangential contact implicit in grazing suggests that any inaccuracy will (slightly) alter the grazing point prediction.

Background

Trajectory sensitivities

Shooting method algorithms, which form the basis for locating grazing phenomena, require the 4 The point where the trajectory encounters the triggering hypersurface can be accurately located by adapting the implicit integration solver. Normal progression from time t k to t k+1 = t k + t step , for a given time-step t step , is achieved by solving a (typically nonlinear) implicit algebraic equation of the form,

F int (x k , y k , x k+1 , y k+1 , t step ) = 0 (21)
for x k+1 , y k+1 . When locating a triggering point, t step becomes a variable, and ( 21) is augmented by the triggering hypersurface equation s(x k+1 , y k+1 ) = 0.

sensitivity of a trajectory to perturbations in parameters and/or initial conditions [START_REF] Stoer | Introduction to Numerical Analysis[END_REF]. To obtain the sensitivity of the flows φ and ψ to initial conditions x 0 , the Taylor series expansions of ( 9)-( 10) are formed. Neglecting higher order terms gives

∆x(t) = ∂x(t) ∂x 0 ∆x 0 ≡ Φ(x 0 , t)∆x 0 (22) ∆y(t) = ∂y(t) ∂x 0 ∆x 0 ≡ Ψ(x 0 , t)∆x 0 ( 23
)
where Φ and Ψ are the sensitivity transition matrices, or trajectory sensitivities, associated with the x and y flows, respectively [START_REF] Frank | Introduction to System Sensitivity Theory[END_REF]. Equations ( 22)-( 23) describe the changes ∆x(t) and ∆y(t) in a trajectory, at time t along the trajectory, for a given (small) change in initial conditions ∆x 0 = [∆x T 0 ∆p T ] T . In accordance with the partitioning (6), Φ can be decomposed as

Φ = Φ x Φ p 0 I ( 24
)
where I is the × identity matrix. An overview of the variational equations describing the evolution of these sensitivities is provided in Appendix A. It should be emphasized that Φ and Ψ do not require smoothness of the underlying flows φ and ψ. As shown in Appendix A, they are well defined for nonsmooth and/or discontinuous flows generated by hybrid systems.

Limit cycles

Poincaré maps

Periodic behavior of limit cycles implies that the system state returns to its initial value every cycle. This can be expressed in terms of the flow as

x * = φ(x * , T ) ( 25
)
where T is the limit cycle period.5 For nonautonomous systems, the period T is a known quantity. However it is not known a priori for autonomous systems. The unknown period, or return time, can be found using Poincaré map concepts [START_REF] Parker | Practical Numerical Algorithms for Chaotic Systems[END_REF][START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF]. These concepts are well known; the following summary is provided as it underlies later analysis of grazing bifurcations in periodic systems.

Γ Σ P(x 0 ) x * φ(x 0 ,t) x 0 σ x ~Fig. 6. Poincaré map.
Referring to Fig. 6, let Σ be a hyperplane that is transversal to the flow φ(x 0 , t), and defined by

Σ = {x : σ T (x -x) = 0} (26)
where x is a point anchoring Σ, and σ is a vector normal to Σ. The return time τ r for a trajectory emanating from x 0 ∈ Σ is therefore given by

σ T (φ(x 0 , τ r ) -x) = 0. ( 27
)
Because the flow φ(x 0 , t) transversally encounters Σ at τ r , the tangent to the flow at that intersection point ∂φ ∂t

(τ r ) = f (x(τ r ), y(τ r )) ≡ f | τr (28) satisfies σ T f | τr = 0. ( 29
)
By differentiating ( 27) with respect to x 0 and τ r ,

σ T Φ(x 0 , τ r )dx 0 + σ T f | τr dτ r = 0, (30)
and taking account of the transversality condition (29), it can be seen that the implicit function theorem [START_REF] Fleming | Functions of Several Variables[END_REF] guarantees the existence of a function τ r (x 0 ) which (locally) satisfies ( 27),

σ T (φ(x 0 , τ r (x 0 )) -x) = 0. ( 31
)
The flow φ and hyperplane Σ together describe a Poincaré map P : Σ → Σ, defined by

P (x 0 ) = φ(x 0 , τ r (x 0 )) ( 32
)
where τ r (•) is given (implicitly) by ( 27). Therefore from (25), a limit cycle of an autonomous system must satisfy

x * = P (x * ) = φ(x * , τ r (x * )). ( 33
)
The corresponding limit cycle is labeled Γ in Fig. 6. For a nonautonomous system, the Poincaré map P : R n → R n becomes simply P (x 0 ) = φ(x 0 , T ), (34) with limit cycles given by x * = P (x * ).

Locating limit cycles

Limit cycles can be located by solving (25) for nonautonomous systems or (33) for autonomous systems. In the case of nonautonomous systems, rewriting (25) gives

F l (x * ) = φ(x * , T ) -x * = 0. ( 35
)
The solution x * can be obtained using a shooting method,6 which solves the iterative scheme

x k+1 = x k -(DF l (x k )) -1 F l (x k ) ( 36
)
where

DF l (x k ) = Φ(x k , T ) -I ( 37
)
and I is the n-dimensional identity matrix.

For autonomous systems, rewriting (33) gives

F l (x * ) = φ(x * , τ r (x * )) -x * = 0. ( 38
)
Again (36) can be used to find the solution x * . In this case

DF l (x k ) = Φ(x k , τ r (x k )) + f | τr(x k ) dτ r dx (x k ) -I. ( 39
)
The term dτ r /dx in (39) can be obtained from (30), giving

DF l (x k ) = Φ(x k , τ r (x k )) -f | τr(x k ) σ T Φ(x k , τ r (x k )) σ T f | τr(x k ) -I (40) = I - f | τr(x k ) σ T σ T f | τr(x k ) Φ(x k , τ r (x k )) -I. ( 41
)
Recall from Appendix A that the sensitivity transition matrix Φ in (37) and ( 41) is well defined for hybrid (piecewise smooth) systems. Therefore, the proposed shooting method is suitable for piecewise smooth limit cycles. This will be illustrated in later examples.

Parameters within the state

It can be seen from ( 9) that the last states within x, those corresponding to parameters p, always remain constant. Therefore the last equations in (35) or ( 38) are always satisfied by default, and hence are redundant. Discarding those equations results in the under-constrained problem,

F l (x) = 0 (42)
where F l : R n → R n-. Furthermore, (24) indicates that the last rows of Φ are identical to the corresponding rows of the identity matrix I. Therefore from (37),

DF l = DF l 0 ×n ( 43
)
where 0 ×n refers to the × n matrix of zeros. Likewise, it is straightforward to show that DF l in (41) also has the structure (43).

It follows that solutions of ( 35) or ( 38) are well defined only when = 0, i.e. no parameters are incorporated into x. In that case F l ≡ F l . It will be shown in Sec. 5.1 though that periodic grazing requires = 1, and in Sec. 6 that continuation requires = 2. Closest grazing point concepts, developed in Sec. 7, require ≥ 2.

Transient Grazing

Mathematical description

Grazing is characterized by a trajectory (flow) of the system touching a triggering hypersurface tangentially. Let the target hypersurface be des

[u T v T ] T that satisfy [b x b y ] u v = 0. ( 45
)
The vector [ ẋT ẏT ] T is directed tangentially along the flow, so it must satisfy (45) at a grazing point. Furthermore, differentiating (8) and substituting (7) gives,

0 = ∂g ∂x ẋ + ∂g ∂y ẏ (46) ⇒ 0 = g x f (x, y) + g y v (47)
where for notational convenience v replaces ẏ.

A single degree of freedom is available for varying parameters (and initial conditions) to find a grazing point. This single degree of freedom can be achieved by parametrization x 0 (θ), where θ is a scalar. Grazing points are then described by combining together the flow definition (9) (appropriately parameterized by θ), algebraic equations ( 8), target hypersurface (44), and tangency conditions (45), ( 47), to give

F g1 (x g , θ, t g ) := φ(x 0 (θ), t g ) -x g = 0 (48) F g2 (x g , y g ) := g(x g , y g ) = 0 (49) F g3 (x g , y g ) := b(x g , y g ) = 0 (50) F g4 (x g , y g , v) := b x b y g x g y (xg,yg) f (x g , y g ) v = 0. ( 51
)
The grazing encounter occurs at time t g along the trajectory, and its state-space location is given by x g , y g . This set of equations may be written compactly as

F g (x g , y g , θ, t g , v) = F g (z) = 0 (52)
where

F g : R n+2m+2 → R n+2m+2 and z = [x T g y T g θ t g v T ]
T . Solution of ( 52) can be achieved using a shooting method, as discussed in the following section.

Shooting method

Algorithm

Numerical solution of (52) using Newton's method amounts to iterating on the standard update formula

z k+1 = z k -(DF g (z k )) -1 F g (z k ) (53)
where DF g is the Jacobian matrix

DF g =           -I 0 Φ dx 0 dθ f 0 g x g y 0 0 0 b x b y 0 0 0 f t b xx + b x f x + v t b yx f t b xy + b x f y + v t b yy 0 0 b y f t g xx + g x f x + vt g yx f t g xy + g x f y + vt g yy 0 0 g y           (54) with f =         f f f . . . f         ∈ R mn×m , v =        v v v . . . v        ∈ R m 2 ×m ( 55
)
g xx =             ∂ 2 g 1 ∂x 2 ∂ 2 g 2 ∂x 2 . . . ∂ 2 g m ∂x 2             ∈ R mn×n , g yx =             ∂ 2 g 1 ∂y∂x ∂ 2 g 2 ∂y∂x . . . ∂ 2 g m ∂y∂x             ∈ R m 2 ×n (56) g yy =              ∂ 2 g 1 ∂y 2 ∂ 2 g 2 ∂y 2 . . . ∂ 2 g m ∂y 2              ∈ R m 2 ×m , g xy =              ∂ 2 g 1 ∂x∂y ∂ 2 g 2 ∂x∂y . . . ∂ 2 g m ∂x∂y              ∈ R mn×m . ( 57
)
The matrices g xx , g yx , g xy and g yy are usually extremely sparse. It has been found that often the error introduced into DF g by ignoring them has negligible effect on convergence. However, situations can arise where these terms are vital for reliable convergence. This is the case, for example, when the trajectory has multiple turning points (peaks and troughs) in the vicinity of the target hypersurface. The example of Sec. 4.4 provides an illustration. Two approaches have been used to obtain the second derivative terms:

1. Numerical differencing. Many simulators provide direct computation7 of g x and g y . Numerical differencing of g x and g y is straightforward, but not particularly efficient for high dimensional systems. 2. Direct computation. By utilizing an object oriented modeling structure [START_REF] Hiskens | Systematic modeling and symbolically assisted simulation of power systems[END_REF], second derivative terms occur only within components. There are no terms introduced by inter-component dependencies. Explicit formulae for second derivative terms can be established for each component model. The sparse matrices can then be efficiently constructed.

Care must be taken in evaluating the terms of ( 52) and ( 54) that relate to trajectory solution. The flow term φ(x 0 (θ k ), t k g ) in ( 48) evaluates, via numerical integration, to the value of x at time t k g along the trajectory that has initial value x 0 (θ k ). Likewise, the terms Φ and f in the first row of DF g should also be evaluated at time t k g along that trajectory. All other terms in DF g should be evaluated at x k g , y k g .

Initialization of variables

As with all iterative procedures, solution of (53) requires a good initial guess z 0 . In terms of the original system variables, initial (approximate) grazing point values of x 0 g , y 0 g , θ 0 , t 0 g , and v 0 are required. These can be obtained from regular simulation.

Referring to Fig. 1, parameter values that are near the critical value result in trajectories that either, (i) encounter the target hypersurface, or (ii) just miss the hypersurface. Therefore, the trajectory induced by parameter θ 0 should be monitored for, 1. a point where b(x, y) = 0, i.e. an intersection with the target hypersurface, or 2. an appropriate local minimum of b(x, y), i.e.

a point where the trajectory passes close by the target hypersurface. This point is given by db/dt = 0, which implies b x ẋ + b y ẏ = 0. Substituting for ẏ from (46), and using (7), gives

db dt = (b x -b y g -1 y g x )f (x, y) = 0.
In both cases, the identified point directly provides initial values for x 0 g , y 0 g , and t 0 g . The corresponding value of v 0 can be obtained from (47) as v 0 = -g -1 y g x f (x 0 g , y 0 g ) with all partial derivatives evaluated at x 0 g , y 0 g .

Example 2 (Switching-time bifurcation)

Referring to Fig. 4, variation of the thyristor firing angle α suggests that a switching-time bifurcation [START_REF] Rajaraman | Nonlinear dynamics and switching time bifurcations of a thyristor controlled reactor circuit[END_REF][START_REF] Jalali | Switching time bifurcations in a thyristor controlled reactor[END_REF] occurs between 125 • and 130 • . The triggering hypersurface in this case is given by i Lr = 0. The shooting method (53) was used to determine the conditions that gave rise to grazing. For illustration purposes, an initial parameter value of α = 120 • was chosen, even though it was known that the actual grazing value α * lay in the range 125 solution process converged in three iterations, with convergence progress given in Table 1.

Figure 7 shows the behavior of i Lr for the initial guess α = 120 • and the grazing value α * = 129.43 • . It can be seen that grazing occurred at time t * g = 0.007946 sec, when the dip during the first conduction period dropped to i Lr = 0. Behavior is ill-defined at that point; i Lr could continue as shown in Fig. 7, or the thyristor could commutate and i Lr remain at zero. Exact behavior would be governed by unmodeled effects.

Figure 4 also indicates that a switching-time bifurcation can occur during the second conduction period, for a value of firing angle in the range 130 • < α * < 135 • . The shooting method (53) was again used to find this grazing value. An initial guess of α = 120 • was used, with solution progress given in Table 2. Even though the initial guess for α was (deliberately) quite poor, convergence was obtained in three iterations. Grazing occurs at t * g = 0.02499 sec, with a firing angle α * = 132.66 Performance specifications often place bounds on transient excursions of system quantities. Determining parameter values that ensure the specifications are satisfied can be formulated as a grazing problem (52). The single machine infinite bus power system of Fig. 8 will be used to illustrate this application.

For this example, the generator was accurately represented by a sixth-order machine model, 8 and the generator excitation system was modeled according to Fig. 9. Note that the output limits on the field voltage E fd are anti-wind-up limits, while the limits on the stabilizer output V P SS are clipping limits [START_REF] Goodwin | Control System Design[END_REF]. Therefore even though this example utilizes a simple network structure, it exhibits nonlinear, nonsmooth, hybrid system behavior. Larger systems are no more challenging. A single phase fault was applied at the generator terminal bus at 0.05 sec. The fault was cleared, without line tripping, at 0.28 sec.

Generators are susceptible to over-voltage protection operation if their terminal voltage rises too high. This may occur during transients following a large disturbance. The field voltage maximum limit E fdmax has a large influence on transient over-voltages. Therefore, this example considers the maximum value of E fdmax that ensures the initial terminal voltage overshoot does not rise above a AVR Exciter PSS

Ṽt

V ∞ E fd Fig. 8. Single machine infinite bus system.

8

The sixth-order machine model consists of two axes, with two windings in each axis [START_REF] Sauer | Power System Dynamics and Stability[END_REF]].

V max V min 1+sT 1 1+sT 2 sT w 1+sT w K PSS ∆ω V PSS 1 1+sT R 1+sT C 1+sT B K A 1+sT A Σ E fd V ref V t - + + E fdmax E fdmin 1+sT 3 1+sT 4
Fig. 9. Excitation system (AVR/PSS) representation.

specified value of 1.2 pu. The target hypersurface in this case is therefore V t -1.2 = 0. Results of the iterative process are given in Table 3, and presented graphically in Figs. 10 and 11. Convergence of the shooting method was achieved in four iterations. This is an encouraging result, as an onerous test condition was chosen. Referring to Fig. 10, it can be seen that the original voltage trajectory is quite flat over the first extended peak. The grazing formulation ( 52) not only describes peaks, but also troughs and points of inflection. It turns out that for this example, (52) actually has three solutions in close proximity. Accordingly, the Jacobian DF g is quite ill-conditioned. It was found that if the second derivative terms in (52) were ignored, the shooting method converged, but to the wrong solution. This occurred for initial conditions over most of the extended peak. Clearly the directional information provided by the second derivative terms is important in cases such as this, where the encounter between the trajectory and the border is not unimodal.

It is evident from Figs. 10 and 11 that this system exhibits quite nonsmooth behavior. In fact fifteen events occur over the initial 2 sec transient, primarily V P SS banging on maximum and minimum limits. Discrete events clearly exert a strong influence on system dynamics. However because trajectory sensitivities, and hence the Jacobian DF g , take those events into account, shooting method convergence is unaffected.

Notice in Fig. 10 that enforcing the performance specification at the first peak has a detrimental effect on the second peak. This reflects the fact that the grazing formulation ( 52) solves for local encounters between the trajectory and the tangent hypersurface. However the specification could be enforced at both peaks by freeing a second parameter, and coupling together two sets of equations of the form ( 48)-( 51).

Periodic Grazing

Algorithm

For a grazing bifurcation to occur in a periodic system, the conditions governing limit cycles, ( 35) or ( 38), and grazing (52) must be jointly satisfied. This requires more degrees of freedom than was the case for transient grazing. In order to locate a limit cycle, the initial point x 0 must be free to vary. Also, grazing point solution requires one free parameter. Computational and notational advantages follow from a consistent treatment of all variables. Therefore, throughout subsequent developments, parameters p ∈ R will be incorporated into x 0 in accordance with (6).

To allow for the extra degrees of freedom, (48) should be restated, F g1 (x g , x 0 , t g ) := φ(x 0 , t g )x g = 0.

(58) Also, recall from Sec. 3.2.3 that with parameters p incorporated into x, the equations describing limit cycles take the form (42). Assembling the full set of equations governing grazing and limit cycles gives

F gl (z) ≡ F gl (x g , y g , x 0 , t g , v) :=         F g1 (x g , x 0 , t g ) F g2 (x g , y g ) F g3 (x g , y g ) F g4 (x g , y g , v) F l (x 0 )         = 0 (59) where z = [x T g y T g x T 0 t g v T ]
T and F gl : R 2n+2m+1 → R 2n+2m+2-. Point solutions of ( 59) require a match between the numbers of equations and variables. This is achieved when = 1, i.e. a single free parameter. When > 1, (59) describes a continuum of solutions. That situation arises in Secs. 6 and 7.

As in the case of transient grazing, numerical solution of (59) using Newton's method amounts to iterating on the standard update formula

z k+1 = z k -(DF gl (z k )) -1 F gl (z k ) ( 60
)
where now DF gl is the Jacobian matrix

DF gl =           -I 0 Φ f 0 g x g y 0 0 0 b x b y 0 0 0 f t b xx + b x f x + v t b yx f t b xy + b x f y + v t b yy 0 0 b y f t g xx + g x f x + vt g yx f t g xy + g x f y + vt g yy 0 0 g y 0 0 DF l 0 0           (61)
with f , v, and the second-order partial derivatives defined by ( 55)-( 57). The term DF l is defined by ( 43).

As noted in Sec. 4.2.1, care must be given to evaluating terms of ( 59) and ( 61) that relate to trajectory solution. Referring to (59), the flow term φ(x k 0 , t k g ) in F g1 evaluates, via numerical integration, to the value of x at time t k g along the trajectory initialized at x k 0 . The associated Jacobian entries Φ and f in the first row of DF gl should also be evaluated at t k g along that trajectory. Similarly, F l in (59) involves the flow term φ(x k 0 , T ) or φ(x k 0 , τ r (x k 0 )) for nonautonomous or autonomous systems respectively. The former evaluates to x at time T along the trajectory initiated from x k 0 . The latter gives x at the point where the trajectory, starting from x k 0 ∈ Σ, re-encounters Σ.9 The associated Jacobian entry DF l , which involves Φ, should be evaluated at time T , or the return point, accordingly.

Example 4 (Nonautonomous limit cycles)

The SVC of Sec. 2.2 was previously used (in Sec. 4.3) to illustrate transient grazing. In steady state, the SVC exhibits (nonautonomous) limit cycle behavior. 10 Figure 12 shows steady state operation for a range of firing angles α. These limit cycles were located using the shooting method (36). It is apparent from Fig. 12 that the limit cycles for α = 100 • and α = 105 • have structurally different forms. A grazing bifurcation occurs somewhere between these two cases. The shooting method of Sec. 5.1 was used to determine the value of α * corresponding to this grazing bifurcation. Initial conditions were given by the limit cycle with α = 100 • . Convergence progress is reported in Table 4. The bifurcation occurs at α * = 102.16 • , and is shown in Fig. 13, along with the initial condition limit cycle.

Figure 13 also shows the post-bifurcation limit cycle that develops when α is incrementally perturbed above the bifurcation value. Interestingly, this single-trough limit cycle coexists with the double-trough version over the range 97.73 • < α < 102.16 • . Both coexisting limit cycles are stable, 10 The SVC model is nonautonomous, due to the time-dependent source.

though the region of attraction for the single-trough one is quite small. At α = 97.73 • , the single-trough limit cycle undergoes its own grazing bifurcation and vanishes.

Example 5 (Autonomous limit cycles)

Robot model

A model of the compass gait biped robot is discussed in detail in [START_REF] Goswami | A study of the passive gait of a compass-like biped robot: Symmetry and chaos[END_REF]. A summary is included here for completeness. The biped robot can be treated as a double pendulum, with point masses m H and m concentrated at the hips and legs, respectively. Figure 14 provides a schematic representation and identifies other important parameters: lengths a and b, and incline angle γ. The robot configuration is described by the support angle θ s and the nonsupport angle θ ns .

The dynamic equations describing the robot can be written

M (θ) θ + N (θ, θ) θ + 1 a g(θ) = 0 (62)
where θ = [θ ns θ s ] T . The state vector is therefore x = [θ ns θ s θns θs ] T ∈ R 4 . The matrix coefficients of (62) are given by, An event occurs when the nonsupport (swinging) leg collides with the ground. This establishes the triggering condition

M (θ) = β 2 -(1 + β)β cos θ d -(1 + β)β cos θ d (1 + β) 2 (µ + 1) + 1 N (θ, θ) = 0 ( 1+β)β θs sin θ d -(1 + β)β θns sin θ d 0
θ ns + θ s + 2γ = 0. ( 63
)
At the event, the nonsupport leg becomes the support leg, and vice-versa. Velocities θns and θs undergo step changes to ensure conservation of momentum through the collision. The resulting reset equations can be written

θ θ + =   0 1 1 0 0 0 Q + (θ - d ) -1 Q -(θ - d )   θ θ - ( 64
)
where

Q + (θ d ) = β(β -(1 + β) cos θ d ) (1 + β)((1 + β) -β cos θ d ) + 1 + µ(1 + β) 2 β 2 -β(1 + β) cos θ d Q -(θ d ) = -β -β + (µ(1 + β) 2 + 2(1 + β)) cos θ d 0 -β .
Equation (64) matches the form (4). The event triggering state y r follows from (63) as y r = θ ns + θ s + 2γ.

Results

To illustrate the shooting algorithm, a triggering hypersurface (border) was established by b(x, y) = x 3 -2.5 = 0 where x 3 ≡ θns . This could be interpreted as a maximum allowable nonsupport leg velocity. (Perhaps the robot falls apart at higher velocities!) Walking motion that just satisfied this constraint was achieved by varying the incline angle γ. Solution progress is given in Table 5, and illustrated in Fig. 15. The figure shows the initial trajectory, and final grazing limit cycle. The Poincaré hyperplane Σ is defined by normal vector σ = [0 0 1 0] T and anchor point x equal to the initial guess for x 0 . This gave Σ = {x : x 3 = 0.1}. A projection of this hyperplane is shown in Fig. 15. The initial trajectory started from, and returned to, this hyperplane. It was not a limit cycle. The final solution started from a point on the hyperplane and returned to that same point. Along the way it grazed the surface b(x, y) = 0. The characteristic multipliers for the grazing limit cycle are: -1.549, -0.1873 and 0.0970. Because one of the eigenvalues lies outside the unit circle, this limit cycle is unstable. Even so, shooting method convergence was fast and reliable. On the other hand, locating this limit cycle by repeated simulations would be extremely difficult.

-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 -2.5 -2 -1.5 -1 -0.5 0

Continuation Method for Grazing

Limit Cycles

Algorithm

When a single parameter is free to vary, (59) describes isolated grazing limit cycles. If two parameters are free, then F gl : R 2n+2m+1 → R 2n+2m , so F gl (z) = 0 describes a 1-manifold, or curve, of solutions. Such curves are useful for exploring parameter dependencies, but cannot be computed directly. Rather, a continuation process is required to generate successive points along the curve. This can be achieved robustly using an Euler homotopy algorithm [START_REF] Garcia | Pathways to Solutions, Fixed Points and Equilibria[END_REF]. A brief description of this predictor-corrector process follows.

Assume the algorithm starts from a point z 1 on the curve.11 The first step of the algorithm, as illustrated in Fig. 16, is the prediction of the next point on the curve. This is achieved by determining the unit vector that is tangent to the curve at z 1 , and moving along that vector at predefined distance κ. This κ is a (scalar) control parameter that effectively determines the distance between successive points along the curve. In regions of high curvature, κ may need to be small. When the curve is almost linear, a large value of κ would suffice. The unit vector η ∈ R 2n+2m+1 that is tangent to the curve at z 1 is given by,

DF gl (z 1 )η = 0 (65) η 2 = 1 (66)
where DF gl is the Jacobian (61) which, because of the two free parameters, now has dimension (2n + 2m) × (2n + 2m + 1).12 The prediction of the next point on the curve is given by,

z p = z 1 + κη. ( 67
)
Having found the prediction point, the next step is to correct to a point z 2 on the curve. The Euler homotopy method does this by solving for the point of intersection of the curve and a hyperplane which passes through z p and that is orthogonal to η. Points z on this hyperplane are given by, (zz p ) T η = 0 (68)

or alternatively

(z -z 1 ) T η = κ. ( 69
)
Either ( 68) or ( 69) can be used. The point of intersection of the curve and the hyperplane is then given by,

F cont gl (z) = F gl (z) (z -z 1 ) T η -κ = 0. ( 70
)
Note that in (70), z 1 , η and κ are fixed, with z being the only unknown. The first 2n + 2m equations ensure the point is on the desired curve. The last scalar equation ensures the point is on the hyperplane. Altogether (70) forms a set of 2n + 2m + 1 equations in the same number of unknowns. They can be solved for z 2 using a standard Newton-based shooting method, of the form (60), with

DF cont gl = DF gl η T . ( 71
)
The complete predictor-corrector process is illustrated in Fig. 16.

After the second point z 2 on the curve has been determined, an approximate tangent vector can generally be used for obtaining successive points. The approximate tangent vector at the ith point, which is used to calculate the (i + 1)th point, is given by,

η i = z i -z i-1 z i -z i-1 2 . ( 72
)
Obtaining the approximate tangent vector involves much less computation than finding the exact tangent vector using ( 65), (66). However the approximation may not be adequate in regions of high curvature.

Example 6 (Continuation)

Illustration of the continuation process again utilizes the SVC model in Sec. 2.2. The aim is to explore the relationship between values of L r and C that give rise to grazing. In Sec. 5.2, an isolated grazing limit cycle was found by allowing the firing angle α to vary, while holding all other parameters constant. That case, corresponding to α = 102.16 • , L r = 1.66 mH and C = 1.5 mF, provides an initial point for the continuation process. Now though, the firing angle is fixed at α = 102.16 • , and L r and C are allowed to vary. The Euler homotopy algorithm provides a sequence of points along the resulting 1-manifold. Figure 17 shows the outcome of this continuation process. Notice that the initial point is located quite close to a sharp turn in the curve. Even so, the continuation process proceeded reliably.

Grazing trajectories typically separate different forms of dynamic behavior. The curve of grazing points in Fig. 17 therefore separates regions where behaviors are structurally different. To illustrate, Fig. 18 shows limit cycles corresponding to parameter values C = 1.51 mF, and L r = 1.3 mH, 1.4 mH and 1.55 mH. The parameter space locations of these three cases are identified in Fig. 17 that the parameter set (L r = 1.4 mH) from the region inside the curve gives rise to a limit cycle that differs greatly from the two cases that lie outside the curve.

Closest Grazing

It has been established previously that when ≥ 2, i.e. two or more parameters are free to vary, F gl : R 2n+2m+1 → R 2n+2m+2-describes an ( -1)manifold, or hypersurface. Let P denote the projection of that manifold onto parameter space. As indicated in Fig. 18, limit cycles corresponding to parameters p 0 / ∈ P do not undergo grazing. However the parameter space "distance" between p 0 and P provides an indication of the closeness, or robustness, of that limit cycle to grazing phenomena. A small distance would suggest that the system was vulnerable to grazing. Therefore identifying points on P that are (locally) closest to a nominal parameter set p 0 provides important design information.

Mathematical formulation

In order to establish the closest grazing minimization problem, it is helpful to distinguish the parameters in (59). It follows from the partitioning (6), which specifies x T = [x T p T ], that (59) can be restated A denotes the A-norm (pp 0 ) T A(pp 0 ). For this to be a valid norm, A must be (symmetric) positive definite. It is common for A to be diagonal, with entries taking account of scaling differences between parameters.

F gl (z, p) = 0 (73) where z = [x T g y T g x T 0 t g v T ] T ∈ R k and k = 2n + 2m + 1 -.
The Lagrangian corresponding to (74)-( 75) is

L(z, p, λ) = 1 2 p -p 0 2 A + λ T F gl (z, p), ( 76
)
and therefore optimal solutions are given by

∇L(z, p, λ) =         ∂F T gl ∂ z λ A(p -p 0 ) + ∂F T gl ∂p λ F gl (z, p)         = 0 (77)
where ∂F gl /∂ z is the (k + 1) × k matrix constructed from DF gl in (61) by removing the columns corresponding to parameters p. Those columns form the (k + 1) × matrix ∂F gl /∂p. Well defined solutions of (77) require rank{∂F gl /∂ z} = k.

To understand the significance of (77), consider the hyperplane that is tangent to the manifold described by (73). Vectors (dz, dp) on that tangent hyperplane satisfy

∂F gl ∂ z dz + ∂F gl ∂p dp = 0. ( 78
)
Premultiplying by the Lagrangian multiplier λ T results in

λ T ∂F gl ∂ z dz + λ T ∂F gl ∂p dp = 0. ( 79
)
But it follows from (77) that the coefficient of dz in (79) equals zero. Therefore λ T (∂F gl /∂p)dp = 0, implying that the parameter space vector (∂F gl / ∂p) T λ ≡ N (p) is normal to P at (z, p).

It may be concluded that (77) describes points on the manifold (73) where the normal vector N (p) aligns with A(pp 0 ).

Algorithm

An iterative procedure for finding the closest (in parameter space) grazing point can be established by utilizing the fact that vectors A(p-p 0 ) and N (p) align at that closest point. The following algorithm is adapted from [START_REF] Dobson | Distance to bifurcation in multidimensional parameter space: Margin sensitivity and closest bifurcations[END_REF]. It converges locally if P is convex (at least locally) or if it is only slightly concave.

Parameter values along the ray emanating from p 0 and in the direction of unit vector ρ ∈ R are given by p = p 0 + θρ, where θ is a scalar. As parameters vary along that ray, the grazing hypersurface P will be encountered at a point described by F gl (z, p 0 + θρ) = 0. This is a slightly modified version of (59), and can be solved using a shooting method similar to (60). The following iterative procedure seeks to minimize ||pp 0 || A .

1. Choose an initial direction ρ 1 . For example, this initialization could be achieved by freeing a single parameter and holding other parameters fixed. 2. Given ρ i , compute the first grazing point in the direction ρ i ; that is, compute θ i such that p i = p 0 + θ i ρ i ∈ P. 3. Compute the vector N (p i ) = (∂F gl /∂p) T λ normal to P at p i . The Lagrangian multiplier vector λ is given by the first equation in (77), and can be obtained directly from the factorization ∂F gl /∂ z = QR. The last column of Q is the desired λ. Computation of N (p i ) is then straightforward.

4. Set ρ i+1 = A -1 N (p i ) and normalize. Iterate until convergence of ρ i to ρ * . The (locally) closest grazing point occurs at p * = p 0 + θ * ρ * ∈ P.

The parameter space direction ρ * of a closest grazing point is aligned with the normal vector N (p * ), and is the fixed point of the iterations. Note that in the special case of P being a hyperplane, iterations converge in a single step.

Example 7 (Closest grazing point)

The SVC example again provides a useful illustration. The algorithm in Sec. The initial direction was chosen as ρ 1 = [1 0] T . Successive points generated by the algorithm are plotted in Fig. 19. (The figure is scaled to match the scaling introduced by A.) It can be seen that convergence was fast, even though a poor initial direction was specified. The vector p *p 0 is shown as a dashed line in the figure. In the scaled coordinate system, this vector is normal to the continuation curve P at the point p * .

Conclusions

Hybrid (piecewise smooth) dynamical systems are characterized by intrinsic interactions between continuous dynamics and discrete events. Such systems are susceptible to grazing phenomena. Grazing refers to situations where the system trajectory tangentially (rather than transversally) encounters an event triggering hypersurface. Behavior beyond the grazing point is often not well defined.

Grazing points can be formulated as a boundary value problem, consisting of a set of nonlinear, algebraic equations that incorporates the system flow. Iterative solution via Newton's method requires numerical integration of the system trajectory, and therefore has the form of a shooting method. The Jacobian required by Newton's method incorporates trajectory sensitivities, which can be efficiently computed along with the trajectory. The shooting method is therefore practical for arbitrarily large hybrid systems.

As the number of parameters increases, the grazing point formulation gives rise to a continuum of solutions. A predictor-corrector continuation process has been developed to follow 1-manifolds, or curves, of grazing solutions. Furthermore, the grazing point formulation underlies the development of an algorithm for finding the smallest (in an A-norm sense) parameter change necessary to induce grazing. Such information is valuable in determining the robustness of a design (set of parameters) to grazing phenomena.

The paper has developed and illustrated shooting, continuation and optimization methods for both transient and periodic grazing. Examples are drawn from power electronics, power systems, and robotics.

Let (x(τ ), y(τ )) be the point where the trajectory encounters the hypersurface s(x, y) = 0, i.e. the point where an event is triggered. This point is called the junction point and τ is the junction time. Assume that the trajectory encounters this triggering hypersurface transversally.

Just prior to event triggering, at time τ -, x and y are given by x -≡ x(τ -) = φ(x 0 , τ -) y -≡ y(τ -) = ψ(x 0 , τ -) where g -(x -, y -) = 0.

Similarly, x + , y + are defined for time τ + , just after the event has occurred. [START_REF] Hiskens | Trajectory sensitivity analysis of hybrid systems[END_REF] show that the jump conditions for the sensitivities Φ are given by Φ(x 0 , τ + ) = h *

x Φ(x 0 , τ -)f +h * x f -τ x 0 (A.6) where h * x = h xh y (g - y ) -1 g - x τ - τ x 0 = -s xs y (g - y ) -1 g - x τ -Φ(x 0 , τ -) s xs y (g - y ) -1 g - x τ -f - f -≡ f (x(τ -), y -(τ -)) f + ≡ f (x(τ + ), y + (τ + )).

The sensitivities Ψ immediately after the event are given by Ψ(x , τ + ) = -(g + y (τ + )) -1 g + x (τ + )Φ(x 0 , τ + ).

Following the event i.e. for t > τ + , calculation of the sensitivities proceeds according to (A.1)-(A.2), until the next event is encountered. The jump conditions provide the initial conditions for postevent calculations.

Fig. 1 .

 1 Fig. 1. Trajectory grazing triggering hypersurface.

Fig. 2 .

 2 Fig. 2. Static var compensator (SVC). (a) One-line diagram. (b) Circuit representation.

  Fig. 3. SVC model behavior.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Thyristor current waveforms for various firing angles.

  b : R n+m → R. Vectors that are normal to b are therefore given by ∇b = [∂b/∂x ∂b/∂y] T ≡ [b x b y ] T , and the tangent hyperplane is spanned by vectors

Fig. 10 .

 10 Fig.10. Terminal voltage, V t .

Fig. 11 .

 11 Fig. 11. Generator field voltage, E fd .

Fig. 12 .

 12 Fig. 12. SVC steady state behavior over a range of firing angles.

Fig. 13 .

 13 Fig. 13. SVC grazing bifurcation.

Fig. 14 .

 14 Fig. 14. Compass gait biped robot.

Fig. 15 .

 15 Fig. 15. Nonsupport leg angle θns versus velocity θns.

  Fig. 16. Predictor-corrector process.

Fig. 18 .

 18 Fig. 18. Limit cycles for C = 1.51 mF and three values of Lr.

  7.2 was used to find values for parameters p = [L r C] T that give rise to grazing, and that are closest, in the A-norm sense, to the nominal values p 0 = [1.2 1.5] T . Notice from Fig.17that variations in L r are about ten times greater than for C. Therefore, for the sake of the illustration, this scaling difference was compensated via the A-norm, by choosing A =

Fig. 19 .

 19 Fig. 19. Parameter space view of closest grazing iterations.

Table 1 .

 1 Shooting method convergence, first conduction period.

		Grazing Values
	Iteration	Firing Angle, α	Time, tg
	0	120.00	0.007900
	1	128.28	0.007954
	2	129.39	0.007947
	3	129.43	0.007946

• < α * < 130 • . The

Table 3 .

 3 Shooting method convergence.

		Grazing Values
	Iteration	Param, E fdmax	Time, tg
	0	5.80	1.14
	1	3.16	1.28
	2	4.37	1.19
	3	4.72	1.20
	4	4.78	1.21

Table 4 .

 4 Shooting method convergence to grazing bifurcation.

		Limit Cycle
		Initial Point, x 0	Grazing Values
	Iteration	i Ls	vc	Firing Angle, α Time, tg
	0	3.8462 -0.5853	100.00	0.01050
	1	4.6489 -0.9602	102.13	0.01043
	2	4.6628 -0.9640	102.16	0.01043
	3	4.6627 -0.9640	102.16	0.01043

Table 5 .

 5 Shooting method convergence to grazing bifurcation.

		Limit Cycle Initial Point, x 0	Grazing Values
	Iteration	θns	θs	θns	θs	Slope (deg), γ Time, tg
	0	-0.3500 0.2100 0.1 -1.0000	3.00	0.2800
	1	-0.4076 0.2253 0.1 -1.1408	4.66	0.2742
	2	-0.4107 0.2223 0.1 -1.1278	5.07	0.2762
	3	-0.4095 0.2220 0.1 -1.1215	4.99	0.2764

  . Notice

		2
	1.5
		1		L r =1.3mH
	0.5
		0	0	0.002	0.004	0.006	0.008	0.01	0.012	0.014	0.016
	(A)	4
	Thyristor current, i Lr	0 1 2 3	0	0.002 =1.4mH L r	0.004	0.006	0.008	0.01	0.012	0.014	0.016
		2
	1.5
		1		L r =1.55mH
	0.5
		0	0	0.002	0.004	0.006	0.008	0.01	0.012	0.014	0.016
							Time (sec)
								1.54
								1.535
								1.53
							Capacitance, C (mF)	1.51 1.515 1.52 1.525
												Original parameter values
								1.505
								1.5
									Cases considered in Figure 18
								1.495
								1.1	1.2	1.3	1.4	1.5	1.6	1.7	1.8
											Inductance, Lr (mH)

Fig. 17. Curve relating Lr and C, given by the continuation process.

Shooting methods couple Newton's method with numerical integration, see[START_REF] Stoer | Introduction to Numerical Analysis[END_REF]. Such methods are discussed in detail later in the paper.

It is assumed throughout that the algebraic equations g(x, y) = 0 are always satisfied.

Reformulation as a multiple shooting method[START_REF] Stoer | Introduction to Numerical Analysis[END_REF] is straightforward. Likewise for the grazing shooting methods discussed later.

These quantities are required for implicit numerical integration.

This occurs at return time τr(x k 0 ), though that time is not explicitly determined.

This point would be obtained by solving (59), as discussed in Sec. 5.1.

A numerically robust method of obtaining η follows from the factorization DF T gl = QR, where Q is orthogonal and R upper triangular. The last column of Q is exactly the desired η.

* Research supported by the National Science Foundation through grant ECS-0332777. 1 community is more familiar with piecewise smooth systems.

2

Appendix A Variational Equations

Away from events, where system dynamics evolve smoothly, the sensitivities Φ and Ψ are obtained by differentiating ( 7)-( 8) with respect to x 0 . 13 This gives

where f x ≡ ∂f /∂x, and likewise for the other Jacobian matrices. Note that f x , f y , g x , g y are evaluated along the trajectory, and hence are timevarying matrices. It is shown in [START_REF] Feehery | Efficient sensitivity analysis of large-scale differentialalgebraic systems[END_REF][START_REF] Hiskens | Trajectory sensitivity analysis of hybrid systems[END_REF][START_REF] Li | Sensitivity analysis of differential-algebraic equations: A comparison of methods on a special problem[END_REF] that the solution of this (potentially high order) linear, time-varying DAE system can be obtained as a byproduct of solving the original DAE system ( 7)-( 8). Initial conditions for Φ are obtained from (11) as Φ(x 0 , t 0 ) = I where I is the identity matrix. Initial conditions for Ψ follow directly from (A.2), 0 = g x (t 0 ) + g y (t 0 )Ψ(x 0 , t 0 ). Equations (A.1)-(A.2) describe the evolution of the sensitivities Φ and Ψ between events. However at an event, the sensitivities are often discontinuous. It is necessary to calculate jump conditions describing the step change in Φ and Ψ. For clarity, consider a single switching/reset event, so the model ( 1)-( 4) effectively reduces to the form ẋ = f (x, y) ( A . 3 ) 0 = g -(x, y), s(x, y) < 0 g + (x, y), s(x, y) > 0 (A.4)

x + = h(x -, y -), s(x, y) = 0. (A.5) (The switching hypersurface s has been given an explicit form so as to elucidate its role in the jump conditions.)

13 Keep in mind that parameters p are incorporated into x 0 .


Daniel Shevitz

Brad Paden

Lyapunov Stability Theory of Nonsmooth Systems

This paper develops nonsmooth Lyapunov stability theory and LaSalle's invariance principle for a class of nonsmooth Lipschitz continuous Lyapunov functions and absolutely continuous state trajec• tories. Computable tests based on Filipov's differential inclusion and Clarke's generalized gradient are derived. The primary use of these results is in analyzing the stability of equilibria of differential equations with discontinuous right-hand side such as in nonsmooth dynamic systems or variable structure control.

I. INTRODUCTION

There are many systems which have nonsmooth dynamics. Examples include systems with Coulomb friction, contact interactions, and variable structure systems where control inputs are allowed to be discontinuous. It is essential to rigorously analyze these systems and address such issues as the existence of equilibria, their stability, and quaiitative dynamics. As important as nonsmooth systems are in practice, techniques are still lacking for their analysis. All classical existence theorems for ordinary differential equations require vector fields which are at least Lipschitz continuous. The aforementioned examples, and many others, fail this requirement. With respect to these classical techniques, one cannot even define a solution, much less discuss existence of equilibria and stability.

What is needed is a set of tools which allow the analysis of differential equations with discontinuous right-hand sides. The seminal contribution in this area was made by Filipov [START_REF] Filipov | Differential equations with discontinuous right-hand side[END_REF] who developed a solution concept for differential equations whose right-hand sides were only required to be Lebesgue measurable in the state and time variables. Using this framework, theorems were proved for existence, uniqueness, and continuous dependence on initial conditions. One area missing from this program is the stability analysis of equilibria using nonsmooth Lyapunov functions.

Yoshizawa [START_REF] Yoshizawa | Stability Theory by Lyapunov's Second Method[END_REF] developed the Lyapunov theory for Lipschitz potential functions, but this work assumed a continuous vector field and smooth trajectories. In his book [START_REF]Differential Equations with Discontinuous Right Hand Sides[END_REF], Filipov studies the equilibria of differential equations with discontinuous right-hand sides, but deals with smooth Lyapunov functions.

There has been work on Lyapunov stability theory in variable structure systems; see [START_REF] Utldn | Sliding Modes and Their Applications in Variable Structure Systems[END_REF], [START_REF] Decar1o | Variable structure control of nonlinear multivariable systems: A tutorial[END_REF] for a detailed list of references. In variable structure systems one tries to pick a sliding surface and show that nearby trajectories converge first to this surface, then once on the surface, to an equilibrium point. The control is chosen in such a way that the dynamics of the closed-loop system never pass through the surface and all analysis can be done from a single side of the switching surface, and then once the trajectory is on the switching surface and smooth dynamics are restored, ordinary Lyapunov theory is adequate. As a result of this piecewise approach, the nonsmoothness of the trajectory is unimportant and differentiable Lyapunov functions suffice in many cases. For our purposes, the "kinks" are an essential part of the dynamics, and this one-sided analysis is inadequate.

The nonsmooth Lyapunov analysis of equilibria is present in the differential inclusions literature [START_REF] Aubin | Differential Inclusions[END_REF], [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF], [START_REF] Frankowska | Lower semicontinuous solutions of Hamilton-Jacobi-Bellman equations[END_REF], of which Filipov's differential inclusion is a special example. Theorems are developed which rely on the notions of set valued maps and derivatives. An alternative formulation using Dini derivatives was developed by Roxin [START_REF] Roxin | Stability in general control systems[END_REF]. Although very general, both these techniques lack a calculus for computing derivatives of nonsmooth Lyapunov functions. We resolve this problem by exploiting the additional structure ofFilipov's solution to develop more useful tools in this context. Additionally, the motivation for studying general differential inclusions is different in that the fqcus is on integral funnels, in which one searches for monotone trajectories which tend to the equilibrium.

The purpose of this work is to provide some tools of generalized Lyapunov analysis with which the stability properties of nonsmooth dynamic systems can be determined. These tools will turn out to be computable and easily applicable making nonsmooth Lyapunov analysis no more difficult than its smooth counterpart. Specifically, our contribution is weakening the restriction of differentiability to a broad class of nonsmooth Lipschitz continuous Lyapunov function!f; the trajectories are only required to be absolutely continuous: This is important because we will show there are nonsmooth dynamic systems whose equilibria cannot be proved stable using continuously differentiable Lyapunov theory. In addition, nonsmooth Lyapunov functions are natural for nonsmooth dynamic systems.

MATHEMATICAL FRAMEWORK

In this section we review the Filipov solution concept for differential equations with discontinuous right-hand sides, the nonsmooth analysis of Clarke's generalized gradient, and develop a connection via a new chain rule for differentiating regular functions along Filipov solution trajectories.

Filipov Solutions

We consider the vector differential equation

:i: = J(x, t) (l)
where f: Rn x R ____. Rn is measurable and essentially locally bounded. We must first define what it means to be a solution of this equation.

Definition 2.1 ( Filipov ):

A vector function x ( •) is called a solution of (1) on [to, t 1 ] if x(•)
is absolutely continuous on [to, tl] and for almost all t E (to , h]

x E K[f](x, t) (2)
where

K[f](x, t) := n n co f(B(x, 6)-N, t ) (3) 8 > 011-N =O
and n i-'N=O denotes the intersection over all sets N of Lebesgue measure zero. An equivalent definition [START_REF]Differential equations with second members discontinuous on intersecting surfaces[END_REF], [START_REF] Paden | A calculus for computing Filipov' s differential inclusion with application to the variable structure control of robot manipulators[END_REF] is: there exists

Nf C Rm, J1NJ = 0 such that for all N C R m, 11N = 0 K[f](x) = co{!imf(x;) I x; ____. x, x ; rf. Nf UN}. ( 4
)
The content of Filipov' s solution is that the tangent vector to a solution, where it exists, must lie in the convex closure of the limiting values of the vector field in progressively smaller neighborhoods around the solution point. It is important in the above definition that we discard sets of measure zero. This technical detail allows solutions to be defined at points even where the vector field itself is not defined, such as at the interface of two regions in a piecewise

f[B(x,o)-N]

Fig. I. The limiting procedure used to calculate K[f)(x, t) .
defined vector field. Fig. 1 helps illustrate the limiting procedure used to define I<[!] ( x, t). This figure shows the vector images of a small neighborhood around the base point x. The interface is a neglected set of measure zero where the vector field is not defined. The set

I<[f](x, t
) reduces to the convex hull of two limit vectors as the neighborhood becomes vanishingly small. Example: A nontrivial example of a nonsmooth dynamic system is the following defined for almost all x E R 2 ( 0 -1)

x = 1 0 Y'llxlh [START_REF]Differential equations with second members discontinuous on intersecting surfaces[END_REF] where llxlh = lx1l + lx2l• The gradient in [START_REF]Differential equations with second members discontinuous on intersecting surfaces[END_REF] is not defined on the x1 and x2 axes, sets of measure zero. Off the axes, the Filipov set

I<[/] ( x) is a singleton { {(-1, +1f} for x E quadrant 1 K[/](x t) = {( -1, -1f} for x E quadrant 2 ' {(+1, -1)T} for x E quadrant 3 {(+1, +1f} for x E quadrant 4. (6)
On the x1 and x2 axes, K[f](x, t) is the convex hull of each of the vectors in (6) corresponding to the quadrants which the point x borders. For example, on the positive x1 axis K[fJ(x) = eo{ ( -1, 1 f, (1, 1)r }. Any Filipov solution to the differential equation (5) traces a square. observe that the trajectories in this example move along level sets of llxll1, a natural (nonsmooth) Lyapunov function for this system.

In most applications the calculation of K[f](x, t) involves functions f(x, t) expressed as sums, products, and compositions of other functions. Hence a calculus is needed for computing the Filipov set. This calculus was derived in our earlier work [10].

We do not require Lipschitz continuity of f(x, t) in x or t and therefore cannot in general expect uniqueness and continuous dependence on initial conditions, although we will assume both in what follows. Readers interested in more detailed discussions are referred to the original literature [4], [5] .

Generalized Gradients

As shown in the example above, nonsmooth Lyapunov functions arise naturally in the stability theory of differential equations with discontinuous right-hand sides. In the application of the machinery of nonsmooth analysis, Clarke's generalized gradient [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF] is particularly useful in simplifying proofs.

Definition 2.2 (Clarke's Generalized Gradient): For a locally Lipschitz function V: Rn x R-+ R define the generalized gradient of V at (x, t) by 8V(x, t) = co{lim V'V(x, t) I (x;, t;)-+ (x, t), (x;, t ;) fi Dv} [START_REF] Frankowska | Lower semicontinuous solutions of Hamilton-Jacobi-Bellman equations[END_REF] where Dv is the set of measure zero where the gradient of V is not defined.

The gradient V' includes the derivative with respect to time ( 8 j 8t). In this definition, Lipschitz means Lipschitz in (x, t) (discontinuities in tare not allowed). For notational brevity, if a function V(x, t ) has no explicit t dependence, we shall adopt the convention of dropping the last component of av which is identically zero.

One way to view the generalized gradient at a point x is a set valued map equal to the convex closure of the limiting gradients near x. For example, the function

V(x) = lx l, x ER
has a generalized gradient which equals

av c x) = { -1} x E R- = {+1} X ER+ = [-1, 1] X= 0. (8) (9)
The next lemma states that the generalized directional derivative, which is defined presently, is the support function in the sense of convex analysis for the generalized gradient. Both the definition and result are due to Clarke [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF].

Definition 2.3: The generalized directional derivative is defined

f o( ) 1 .
f(y + tv)f(y )

x; V = 1msup . y~x , ttO t ( 10
)
Lemma 2.1: Let f be Lipschitz near x, then

rcx; v) = max{(~, v) I~ E 8f(x)}. ( 11
)
Our chain rule will be for a useful class of functions, called regular functions [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF].

Definition 2.4: f(x, t ): Rm X R-+ R is called regular if 1)
for all v, the usual one-sided directional derivative !' (x; v) exists,

2) for all V, j'(x; v) = r(x; v).

Examples of regular functions include smooth functions and functions which can be written as the pointwise maximum of a set of smooth where

V(x, t): = n ~T (K[f](~(t), t)).
eEBV(x(t) , t ) [START_REF] Vidyasagar | Nonlinear Systems Analysis[END_REF] Proof That V(x(t), t) is absolutely continuous and (dfdt)V(x(t), t) exists almost everywhere is a consequence of the fact that the composition of a Lipschitz function, V ( •), with an absolutely continuous function, ( x ( t), t), is absolutely continuous (see, for example, [9]). At a point where x(t) and V(x(t) , t) are both differentiable (this is true almost everywhere)

d dt V(x(t), t) (14)
= lim V(x(t +h), t +h)-V(x(t), t) hto

h (15) = lim V(x(t) + xh, t +h)+ o(h)-V(x(t), t) htO h (16) (17)
(18)

= V'((x, tl; (x, 1f) = V 0 ((x, tl; (x, 1f) by regularity = max{(~i(x, 1)T)i~ E aV(x, t)}.

a similar argument shows lim _V_,_( x__,_(_t +_h:..:.. ) ,_t _+::-h_,_) _-_V -'('-x~( t:..:..) ,-'-t)

h!O h = min { (~i(x, 1)T )i~ E aV(x, t). Therefore :t V(x(t), t) =(~I (x, 1)T) V~ E aV(x, t). Now, x is a Filipov solution so that x(t) E K[f](x(t)), a.e. ( 20
) (21) (22) (23)
Thus almost everywhere, V = e (i) for all~ E aV(x, t) and some Tf E K[f](x, t). Equivalently

V e•e.Y = n ~T (K[f] (~(t), t)).

{ E8V(x(t), t) D (24)

Example: Let V(x) = llxll1 and f(x) as in the example, (5), above. Let x(t) be the solution passing through (1, O)T at time t = 0. We have

aV(x(O), 0) = C-1, 1 +1]) (25)
(where we have employed our convention of dropping the last term in aV(x, t) when V is independent of time) and

K[f](x(O) , 0) = e-1 • 1 +l]). (26)
Let (1, 6), 6 E [-1, 1] be an arbitrary element of aV(x(O), 0)

then [; 2 ] K[j](x(O), 0) = 6 + [-1, +1] = [6 -1, 6 + 1] (27) implies V(x(O) , 0) = n [6-1, 6 + 1] = 0. (28) { 2E[ -1, 1]
'fhis does not guarantee that V = 0 (or even exists) at this point, but V ( x) ~ 0 by which we mean v < 0 for all v E V, does guarantee stability. The theorems of the next section formalize this.

Ill. STABILITY THEOREMS

In this section we state two existing Lyapunov theorems (uniform stability ~d uniform asymptotic stability) in terms of the set valued map V. The proofs are omitted because they are identical to their smooth counterparts except for some relations holding "almost everywhere" instead of everywhere. See Khalil and Vidyasagar [8], [START_REF] Vidyasagar | Nonlinear Systems Analysis[END_REF] for the smooth versions. Proofs exist for the other versions of stability, such as nonuniform, global, etc., but for the sake of brevity we have omitted them. An application of this theorem is made, and the stability of the spring-mass-coulomb friction system is proved. Theorem 3.1: Let x = f(x , t) be essentially locally bounded and 0 E K[f](O, t) in a region Q :J {x E Rn l ll x ll < r} X {tito ~ t < oo} . Also, let V: Rn x R --+ R be a regular function satisfying [8] for a definition of class K functions.) Then 1) V(x, t) ~ 0 in Q implies x(t) = 0 is a uniformly stable solution.

V(O,t)=O (29) and (30) in Q for some V 1 , V 2 E class K. (See
2) If in addition, there exists a class lC functions w ( •) in Q with the property

V(x, t) ~ -w(x) < 0 (31)
then the solution x( t) = 0 is uniformly asymptotically stable.

Example ({6], [10]): Let R(x, t) be a matrix which is continuous and uniformly positive definite when symmetrized and

x = -R(x(t), t) \7llxlll• (32)
Then x(t) = 0 is asymptotically stable. Proof: Choose V(x) = llxlh-Then V(x , t) = n e K[-R(x , t) \7llxlhJ(x, t ) (33) {E&IIxlll (since V is time independent) and by the calculus for J( [10], we have

V(x, t) = n -~T R (x, t )aii xlll • (34) { Eoll xll1
Since R is uniformly positive definite when symmetrized, there exists

p > 0 such that aJJ xlll is convex so satisfies ~o(x, t) = argmin e (R( x , t) + RT (x, t) )~ eE&IIxlh 2
e ( R(x , t) ~ RT (x, t) )~a (x, t )

;:::: ~o(x, t)T ( R(x, t) ~ R T(x, t) ) eo(x , t)

;::::

P li~o (x , t)ll 2 (35) (36)
(37) for all ~ E a 11 X Il l. Equation ( 37) is a consequence of the fact that for a convex domain D and a smooth function f at the point xo = arg min f ( x ), the point so V ~ -p almost everywhere. This implies x(t) ---+ 0 asymptotically (converges in finite time in fact by the finite time stability theorem in [10]). The proof just demonstrated shows, with minor modifications, the example in the last section where

R(x, t) = (~ ~1 ) (41)
has zero as a stable equilibrium point. We now argue that this equilibrium cannot be shown stable by any smooth time independent Lyapunov function.

The trajectories of this dynamic system are squares. Any nonincreasing function along a closed curve must be a constant, therefore the level sets of any proposed Lyapunov function must be squares. The level sets of any smooth function with nonzero gradients are smooth precluding the possibility of any smooth Lyapunov function for this example. A similar argument also shows the impossibility of finding a continuously differential time dependent Lyapunov function.

We now prove a nonsmooth version of LaSalle's theorem. Theorem 3.2 (LaSalle): Let n be a compact set such that every Filipov solution to the autonomous system x = f(x), x(O) = x(to) starting in 0 is unique and remains in 0 for all t 2:: to. Let V: 0 ---+ R be a .time jndependent regular function such that v ~ 0 for all

v E V (if V
x = f(x), xo
x(to). Since s(t, xo, t9) lies in L+, (d/dt)V(s(t, xo, to)) = 0. Since V(x(t)) E V(x) a.e., in t, we have s(t, x 0 , t 0 ) E S a.e. Thus s(t, x 0 , t 0 ) E S for some t E [to, t 0 +b) and all b > 0. Since s(t, x 0 , t 0 ) is continuous (absolutely) we have that there are points in S arbitrarily close to x 0 implying x 0 E S. Hence L + C S. See Fig. 2 for a phase portrait.

IV. CONCLUSION

In this paper we have extended basic Lyapunov stability theorems to the nonsmooth case using the Filipov solution concept and Clarke's generalized gradient. The result is a theory applicable to systems with switches for which natural Lyapunov functions are often only piecewise smooth. This machinery should find application in variable structure control theory, the analysis and control of mechanical systems, and the analysis of pulse width modulated control systems.

  x 0 satisfies \7f(xo)•(y-xo);::::O forallyED .At .r = 0, 8llxll1 = [-1, +1t. Since llxll1 is convex, 8llxll1l:r:;to n (-1, +1)" = 0 (this follows from [2, Proposition 2.2.9]

0

  The following theorem is a special case of a viability result in the study of differential inclusions (see[1, p. 180] for a general statement and complete proof). Theorem 3.4: If M is an invariant set in a smooth k-dimensional manifold S, then TmS n K[/](m) I 0 (43) for all m E M.

2 3 XFig. 2 .

 232 Fig. 2. The phase portraitofmx+bsgn (:i:)+kx = 0 with m= b = k = 1.

  It follows by the lemma below that L + C S. Since L + is contained in the largest invariant set in S the theorem is proved. 0 Let xo E £+ and s(t, xo, to) be the solution to

	Lemma 3.3: Under the conditions of the theorem above
	(42)
	Proof•

is the .empty set then this is trivially satisfied). Define s = {x E n I 0 E V}. Then every trajectory inn converges to the largest invariant set, M, in the closure of S.

Proof• Let x(t) be a Filipov. solution starting in 0. Since V(x(t)) is absolutely continuous, V is bounded below zero, and V is bounded above zero, V(x) tends to a constant, a, as t ---+ oo.

Uniqueness of trajectories implies continuous dependence on initial conditions so the positive limit set, L +, of x ( t) is an invariant set. Moreover, L+ C 0 because 0 is closed. By continuity of V, V (p) = a for all p E L +. Since L + _is invariant we have V = 0 on L + , hence by Theorem 2.2, 0 E V ( x), a. e. in t.

The authors are with the Mechanical Engineering Department, University of California at Santa Barbara, Santa Barbara, CA 93111 USA.