Automated generation of hybrid automata for multi-rigid-body mechanical systems and its application to the falsification of safety properties

ABSTRACT What if we designed a tool to automatically generate a dynamical transition system for the formal specification of mechanical systems subject to multiple impacts, contacts and discontinuous friction? Such a tool would represent an advance in the description and simulation of these complex systems. This is precisely what this paper offers: Dyverse Rigid Body Toolbox (DyverseRBT). This tool requires a sufficiently expressive computational model that can accurately describe the behaviour of the system as it evolves over time. For this purpose, we propose an alternative abstraction of multi-rigid-body (MRB) mechanical systems with multiple contacts as an extended version of the classical hybrid automaton, which we call MRB hybrid automaton. One of the chief characteristics of the MRB hybrid automaton is the inclusion of computation nodes to encode algorithms to calculate the contact forces. The computation nodes consist of a set of non-dynamical discrete locations, discrete transitions and guards between these locations, and resets on transitions. They can account for the energy transfer not explicitly considered within the rigid-body formalism. The proposed modelling framework is well suited for the automated verification of dynamical properties of realistic mechanical systems. We show this by the falsification of safety properties over the transition system generated by DyverseRBT.


Motivation and scope of the problem
This paper proposes a novel computational modelling framework that facilitates the automated generation of transition systems of complex mechanical systems for which a fully and clear description is difficult to obtain. We model mechanical systems with multiple impacts, contacts and discontinuous friction as hybrid dynamical systems, specifically, as hybrid automata [1][2][3][4], which provide a suitable modelling framework for the specification and analysis of complex hybrid systemswhere continuous and discrete, smooth and abrupt parts interact with each other.
The automated generation of the dynamic relationships and discrete transitions of the systems treated here are possible by the proposal of a computational modelcalled the multi-rigid-body (MRB) hybrid automatonthat can accurately describe the behaviour of the system as it evolves over time. A method of converting a certain class of MRB problems into an MRB hybrid automaton is proposed. The derived hybrid automaton expresses the rigid-body system in a standard form with a number of special elements termed computation nodes. The contribution of this paper is to provide an alternative solution to the problem of the formal specification of complex dynamic interactions in MRB systems, with the ultimate goal of the automated generation of a discrete description to be able to simulate the system's dynamic behaviour. By building on previous mechanical models for MRB systems with multiple impacts, contacts and discontinuous friction [5][6][7][8][9][10][11][12][13][14], we offer a higher level of abstraction to specify the multiple transitions (dynamic and discrete) on these systems. These transitions and dynamic relationships grow with an increasing number of rigid bodies and multiple contacts and can get unmanageable and not clearly specified with classical mechanical models. This is overcome with the MRB hybrid automaton and its automated generation.
Although there are restrictions and limitations in the systems treated in this paper (see Section 1.2), they are broad enough to validate this paper as a novel and worthwhile contribution for the hybrid system modelling and engineering communities. For example, the growing interest in formal methods and hybrid systems in the robotics community is manifested in a wide range of recent works [15]. Formal specification is also well suited for autonomous vehicle control [16], robot motion planning and control [17,18] and control of semiautonomous vehicles [19]. Furthermore, the modelling and analysis framework built in this paper is applicable to a broader class of networked complex systems, as it has been recently shown for the case of neuronal networks [20].
Creating computational models of MRB systems is not straightforward. The governing dynamic equations are piecewise smooth and contain discontinuities and jumps in the state space resulting from changing friction and impacts. Moreover, the governing equations themselves may be non-linear and contain transcendental functions. The possible occurrence of multiple contacts further expounds the problem.
Hybrid automaton is a powerful computational-oriented framework for specifying the dynamic behaviour of rigid-body mechanical systems. Thus far, however, work has mainly focused on examples with only single points of contact. Consider for instance, the ubiquitous bouncing-ball example (see Ref. [21] and references therein). Modelling mechanical systems using hybrid automata becomes more difficult when there is the possibility of multiple contacts. Consider two separate objects sharing a set of contact points along their respective surfaces. Now, consider a third object which collides with one of the two original objects. We know that the collision, despite being separate, will affect the behaviour at the site of the other contacts between the first two objects. Physically, we know this to be caused by energy transfer occurring over extremely short timescales. However, when we model mechanical systems, we often simplify the problem by assuming rigid bodies that neglect the microscopic and fast-acting internal physical processes. The transition of energy from one contact site to the next is assumed instantaneous, and the mechanism by which this process occurs is invisible to us. We must instead resort to nonanalytical numerical procedures, such as linear and non-linear complementarity or proximal point methods, to obtain solutions to the contact forces [5][6][7]9,11,13,14,22].
In this paper, we propose to model mechanical systems as hybrid automata by supporting multiple contacts in the system that occurs simultaneously. The hybrid automaton assumes that Newton's impact law (restitution) and Coulomb friction are acting at the points of contact. Multiple contacts are addressed by introducing what we call computation nodes. This is one of the key new elements that we add to the classical hybrid automaton modelling framework. Computation nodes are a collection of non-dynamical discrete states, transitions, guards and reset functions which are used to find valid combinations of contact forces for each configuration of the mechanical system. The combination of all these elements can account for the energy transfer not explicitly considered within the rigid-body formalism. In the computation node, the executions of the hybrid automaton are assumed to operate over zero time as the discrete states are non-dynamical. The physical analogue of this is the fast-acting dynamics acting at a microscopic level which we assume to occur instantaneously in the rigid-body mechanics formalism. The inclusion of the computation nodes gives us the flexibility to encode more complex and even iterative algorithms into the hybrid-automaton framework, with the benefit that we can model more complex mechanical systems with multiple contacts.
The resulting computational-mechanical model is referred to as the MRB hybrid automaton, which consists of two types of discrete locations: dynamical discrete locationsthe classical discrete locations in a hybrid automaton, that is locations with an associated dynamical systemand non-dynamical discrete locationsthe new type of discrete locations introduced in our modelling framework. This automaton provides a complete and precise modelling framework and fully describes the transitions, transition conditions and operating modes present in the class of mechanical systems under study. The preliminary idea of the MRB hybrid automaton was proposed in other work [23]. The MRB hybrid automaton is a modification of the discontinuous dynamical systems (DDSs) hybrid automaton [24] which models DDSs with discontinuous state derivatives (e.g. systems with dry friction)to include jumps in the states (a discontinuity usually associated with systems with impacts). These modifications are inspired by mythical modes in discontinuous systems, originally proposed in Ref. [25] and further expanded to hybrid systems by Refs. [26,27].
The hybrid automaton modelling framework proposed to model mechanical systems with multiple impacts and friction has two main potential applications. First, the generation of event-driven simulations of MRB mechanical systems. Second, the MRB hybrid automaton is well suited for the formal verification of dynamical properties of realistic mechanical systems. Due to the complexity involved, verifying these systems manually is laborious and could be eased by automatic computational techniques.
The automatic generation of MRB hybrid automata is a necessary step in the formal analysis of a complex mechanical system. Moreover, we highlight its useful application in the simulation of MRB systems, in particular with respect to event-driven schemes. For this purpose, we have produced a tool called Dyverse Rigid Body Toolbox (DyverseRBT). The acronym Dyverse corresponds to DYnamically driven VERification of Systems with Energy considerations, and RBT corresponds to Rigid Body Toolbox. DYVERSE is a novel computational-dynamical framework for the modelling, analysis and control of complex systems [28]. DyverseRBT automatically generates the MRB hybrid automaton and its simulation from a mechanical specification given by the user. In other words, DyverseRBT produces a dynamical transition system that defines the computational semantics of the MRB system.
In the system description given, the events or discrete transitions are defined in order to clearly specify all the possible changes in the dynamics. As a consequence, the hybrid automaton proposed may be easily translate to a program or to any other description language, as it is proposed, for instance, in Refs. [29][30][31][32][33]. We note that the hybrid automaton formulation can be interpreted as an event-driven scheme (or event-tracking time-stepping scheme) for the numerical time-integration of the system; however, an improvement of event-driven schemes is out of the scope of this paper. The alternative numerical time-integration techniques to event-driven schemes are the time-stepping schemes (also known as event-capturing time-stepping schemes). These schemes are not based on the accurate detection of events, are quite robust in contact detection and also have convergence proofs. The reader is referred to Refs. [5,14,34,35] for further details and references on time-stepping schemes for mechanical systems.
The approach proposed here follows a similar rationale to the recent work by Mosterman et al. [53], in the sense that we create computational semantics of the dynamical behaviour of a system which is described by differential equations and works in several modes of operation. The ultimate goal is using computer science methods (like model checking) for analysis purposes, in addition to computing simulations of continuous-time-based models. We highlight that Ref. [53] proposes a general design framework for stiff hybrid systems with special emphasis on control synthesis. Our work could be extended to be used for control synthesis, but it is out of the scope of this paper.
The results presented in this paper must be distinguished from the different simulators for hybrid and cyber-physical systems like Stateflow of Simulink® [54], Modelica [55], Ptolemy [47] or CyPhySim [56]. However, since the MRB hybrid automaton fits within the standard automaton framework, it should be possible to implement the derived hybrid automaton on other hybrid system simulators. A by-product of the automated MRB hybrid automaton generation has been the automated creation of Simulink S-functions. These S-functions implement the hybrid automaton and can be dropped into the Simulink models to interact with other components. Given this, a similar approach could be followed for other simulator tools.
To explore the potential application of our MRB modelling framework and the tool DyverseRBT in automated verification, we produce a method to check in an automated way that some properties are not held for a family of MRB mechanical systems which are subject to multiple contacts, impacts and discontinuous friction and interpreted as MRB hybrid automata. Particularly, we follow an engineering-driven approach, with the main goal of finding bugs in the design rather than proving the correctness of the system. Indeed, from a practical viewpoint, engineering and control designers are interested in automated verification techniques that can provide counterexample trajectoriesthat is, trajectories violating a given propertyin order to consider scenarios difficult to extract from simulations. This is typically referred to as falsification (finding errors), which is different to proving that errors do not exist. In this paper, the analysis of properties in MRB mechanical systems will be done through falsification. Dynamical restrictions will be specified as safety properties ('something bad will never happen'), which will be ultimately reduced to check the satisfiability of Boolean formulas.
The idea of falsification is not new in dynamical systems, see for example Refs. [57][58][59][60] where continuous-time non-linear smooth systems are studied. Particularly, in Ref. [59], there is a focus on the validation of controllers. We highlight the works of Refs. [61][62][63][64][65] on falsification of safety properties in hybrid systems, and the Simulink®-based tool S-TaLiRo for hybrid systems with relatively simple dynamics in each mode or subsystem [66].
The implementation of our especially tailored falsification method for MRB systems is done through the pilot software Dyverse Bounded Model Checker (DyverseBMC), especially designed for mechanical systems with multiple impacts and friction. The input of DyverseBMC is a dynamical transition system generated by DyverseRBT. The approach we adopt for falsification in this paper will mainly use bounded model checking (BMC) [67,68]. In BMC, instead of computing the states that satisfy a specification (property), a SAT (Boolean satisfiability) solver is used to find an error, typically known as a counterexample, which is a trace that does not satisfy the specification. Hence, a property can be assessed not to be satisfied if no counterexample is found; specifically, if the SAT procedure cannot generate a propositional formula that is satisfiable for a counterexample. In DyverseBMC, instead of following the traditional approach of using SAT-based BMC [69,70], we use a lazy-SMT (Satisfiability Modulo Theories) solving approach inspired by the recent works of Refs. [71][72][73][74]. SMT-based BMC is becoming increasingly important in practical applications and hybrid systems, for example iSAT [75,76], MathSAT [77], Absolver [78,79], verification of networks of linear hybrid automata [80], bounded reachability analysis for linear hybrid automata [81] and the recent bounded model-checker dReach [82] that uses the SMT solver dReal [83,84]. Whilst some of these tools (iSAT, Absolver and dReach with dReal) can be mainly applied to non-linear dynamical systems and some types of hybrid systems, they cannot be directly applied to our systemsnon-linear hybrid automata, with a high number of different types of discrete locationsin their present form.
Even though the falsification technique proposed and implemented here via DyverseBMC is not optimized and has limitations, still it represents a valuable solution, considering that not many current automated verification and falsification techniques would be able to be applied without important expansions to the type of systems our MRB hybrid automaton is modelling.
In conclusion, this paper proposes a new way of modelling and simulating of dynamical behaviours of mechanical systems with multiple impacts, contacts and discontinuous friction within the framework of hybrid automata, which is applicable to the analysis of dynamical properties in an automated way. We formulate for the first time a class of MRB mechanical systems with friction and impacts as an expanded type of hybrid automata. The consideration of computation nodes within hybrid automata opens the possibility to explicitly include in the model the computation of key parameters and dynamic relationships for complex systems, with a clear application on the automated generation of computational models for these systems. The conceptual idea of this alternative abstraction of hybrid systems has potential to influence not just the modelling of mechanical systems, but also the modelling and analysis of complex systems with complex computations that cannot be clearly encoded in classical hybrid automata. As an evidence of this statement, we have recently applied the idea of the MRB hybrid automaton and the computation nodes for the modelling of complex neuronal networks [20]. The computation nodes are used to calculate and specify the synaptic weights and the synaptic currents for single neurons or populations of neurons. We think that it is important to separate these computations from the membrane dynamics of neurons. The inclusion of the computation nodes gives us the flexibility to encode more complex and even iterative algorithms into the hybrid-automaton framework, which has the benefit that we can model non-trivial neuronal interactions. A summary of the contribution of this paper can be found in Figure 1.

Assumptions and restrictions
The approach outlined in this paper can, in principle, be applied to the analysis of a broad range of rigid-body systems. However, in the present work, we impose a number of assumptions to simplify the systems under test. We do so as a practical measure, to make the problem solvable and to reduce the corresponding complexity of our analysis tools. We outline these assumptions and restrictions here as follows: Figure 1. Summary of the contribution of this paper and the applications of the proposed modelling framework.
• Mechanical-rigid-body systems with a small number of objects. The method requires the construction of a complete map of all possible contact combinations to completely define the MRB hybrid automaton. This can create very large automata due to combinatorial explosion, which can be impractical to handle. This will be better understood in the examples of Section 3.8.2 and in Table 2, where several multiple contact problems are presented. Thus, rigid-body systems with only a small number of objects (and contact combinations) can be practically studied. • Mechanical-rigid-body systems consisting of spheres. Our definition of the MRB hybrid automaton only accounts for rigid bodies with spherical surfaces. It would be straightforward to extend the definition to other shapes and surfaces; however, we do not do so in the present work. Restricting to spheres keeps the possible number of contact combinations low that is, each pair of objects can share only one possible contactand thus reduces the size of the automata to a more manageable size. • Our hybrid automaton formulation is an event-driven scheme, and as such, is ill-equipped to handle finite accumulations of non-smooth events (like impacts or stick/slip transitions).

Motivating example: automatic generation of the MRB hybrid automaton by DyverseRBT
An example is introduced to facilitate the description of our MRB hybrid automaton construction. We will return to this example at each stage, building up the complexity as we progress. The example here is the Interception Game which is shown in Figure 2. The game involves two balls. The first ball (ball 1, the upper ball in the figure) is unconstrained and is subject only to gravity. The second ball (ball 2, the lower ball) is driven by a proportional feedback controller which constrains the ball to the horizontal plane (y ¼ 0) and tracks the position of ball 1 along the x and z axis. The objective of the game is to use ball 2 to intercept ball 1 as it falls towards y ¼ 0 and bounce it back into free space. If ball 2 does not intercept, then ball 1 falls through the horizontal plane. If ball 1 falls through the plane within some bounded time, then the game is lost.
The system starts in a non-contact state and is governed by the following dynamic equations: (1) Figure 2. Interception game.
where the double dot denotes second derivative with respect to time t, g is the force generated by gravity and K is the proportional feedback gain. The accelerations for the remaining states is set to zero. We also define a set of initial conditions: where the dot denotes derivative with respect to time.
The MRB hybrid automaton is automatically generated by DyverseRBT. The package DyverseRBT uses a Matlab-based program to generate an MRB hybrid automaton of an MRB mechanical system from a minimal set of information on a rigid-body system: a list of entities and their interacting forces. The package can then use the automaton to generate an S-function of the mechanical system for Simulink simulations. This is a potentially powerful domain-specific application, and to our knowledge, nothing similar exists to date.
A schematic of the package is shown in Figure 3. The user inputs a text file consisting of information on a set of rigid bodies, the forces interacting between entities, and external inputs as shown in Figure 4 for the example of the interception game. This file has the extension .drbt. The user also creates a Simulink model containing an empty S-function block. The external inputs in the text file become the inputs to the S-function block in the Simulink model. A function called 'DyverseCreate' reads the input file, extracts the relevant information and generates a hybrid automaton in the form of a Matlab structure, following the formalism that will be defined in the next sections. The hybrid automaton is then converted into C code for an S-function which is compiled using the on-board MEX compiler. Thus, the package generates an S-function for the Simulink model (sim.dll), and a text file to display the automata (display.txt). The S-function can be dropped into a Simulink model and used as an ordinary block. More details are given when we come back to this example in Section 3.8.
The programs that form DyverseRBT and the input files for the examples used in this paper are available at http://staff.cs.manchester.ac.uk/navarroe/research/dyverse/DyverseRBT_BMC.zip.

Mechanics
The primary element in our modelling framework based on the MRB hybrid automaton is the rigid-body object or 'entity', which will include the position, velocity and acceleration in the contact reference frame, in addition to geometrical and physical properties of each body. The notation used is primarily drawn from [13] and will, as far as possible, be kept the same. In the rest of the paper, the subscripts N and T will denote normal and tangential. T as a superscript on a vector or matrix will denote transpose. Definition 2.1: (Entity) An entity E i is a representation of a rigid body consisting of the collection, is a vector containing the position " q i 2 R 3 of the centre of rotation of the entity in a global Cartesian coordinate frame, and the rotation of the entity about the xaxis (α i ), about the y-axis (β i ) and about the z-axis (γ i ) of the global frame. • m i is the total mass of the entity. • J i : R 6 ! R is a strictly convex differentiable function, called the surface function, which defines the surface of the entity as J i q i ð Þ ¼ 0.
The condition of strict convexity of J i means that any pair of entities can only share one contact point. This simplifies the synthesis of the MRB hybrid automaton in the next section. Note that this restriction does not preclude the multiple contact case. An entity can have more than one contact point provided that each of the contact points is with a different entitythat is, each pair of entities shares only one contact point. The condition of differentiability is necessary to ensure that we can calculate a unique solution for the surface normal at any point along the surface. This is critical for our treatment of contact problems.
The conditions on the function J i are heavily restrictive on the types of rigid-body problems we can consider. For instance, the conditions imply that we are only able to address rigid bodies or entities with shapes such as spheres or ellipsoids, but not those with corners or concave shapes. However, it is possible, in principle, to extend what follows to relax these conditions, for example by allowing an entity to consist of multiple surface functions, having rounded corners etc. We do not do so in the present work in order to maintain simplicity.
From now and so on, Á j j denotes the absolute value for any real number, or the cardinality of a set, and jÁj is the 2-norm on the Euclidean space R n . Furthermore, to ease the notation, in the ordinary differential equations (ODEs) presented, we will not write the dependency on time of the variables.
Consider the case of two entities E 1 ; E 2 about to come into contact as shown in Figure 5. Denote this point of contact the ith point of contact: that is, we assume that there are already i À 1 points of contact present in our system of rigid bodies. Denote p 1 as a point in the global reference frame that lies on the surface of E 1 , and p 2 similarly as a point on the surface of E 2 . The points are chosen so that when collision occurs, p 1 ¼ p 2 . We can define a gap function g N;1 which returns the distance between the two points in a direction normal to the surface of E 1 , where n 1 is a vector normal to the surface of E 1 at the point p 1 . We can put the relative velocity and acceleration between the surface points in terms of entity coordinates q 1 ; q 2 as, and w a;N; Â andr i , for all i, is a skew-symmetric matrix. We shall use N in the subscript from here on to indicate that the term refers to the normal component. Note that This can be generalized for an MRB system consisting of M entities E 1 ; E 2 ; . . . ; E M f gwith I ¼ 1; 2; . . . ; i; . . . ; P f gpoints of contact between them occurring at any one time.  Figure 5. Two entities about to contact.
The following terms are introduced: w a;N q; _ q ð Þ ¼ w a;N;1 ; w a;N;2 ; . . . ; w a;N;P À Á T : Then, for the whole system consisting of M entities and P contacts: The resulting Equations (2), (7) and (8) are useful for expressing the constraint that no two rigidbodies can occupy the same space, i.e. penetration cannot occur. We now consider the constraints acting tangentially to the contact surface. Returning to the single contact in Equation (2), define two vectors t x;1 ; t y;1 orthogonal to one another and to the normal vector n 1 . The combination of the three vectors forms a coordinate reference frame at the site of contact. We refer to this frame as the contact frame. The position of point p 2 relative to p 1 projected onto the tangential plane is Following the same procedure as for the normal component, we can generalize for the ith contact and obtain expressions for the relative tangential velocity and acceleration, For the whole system, we obtain where w a;T ¼ w a;T;1 ; w a;T;2 ; . . . ; w a;T;P À Á T : The subscript T from here-on refers to the tangential component. There are two forces acting at the point of contact. The first is the reaction force, which prevents any two rigid bodies from occupying the same space, that is, penetration between entities cannot occur. The second force is friction, which resists motion along the surface of the point of contact.
The reaction force λ N;i acts along the surface normal to the point of contact and is always sufficiently large to enforce the unilateral constraints g N;i ¼ 0, _ g N;i ¼ 0 and € g N;i ¼ 0 at the point of contact. We refer to these as the non-penetration constraints or Signorini's contact law. The reaction force can be expressed by the complementarity condition on the acceleration level as, & which yields to the so-called contact linear complementarity problem (LCP) [8]. We also refer the reader to Chapter 4 of Ref. [85] for an explanation of complementarity systems under the framework of hybrid dynamical systems.
Friction acts on the plane tangential to the point of contact the direction of g x;i and λ y;i is in the direction of g y;i . Here, we use the Coulomb friction given by the friction law in Ref. [86]. The friction force is bounded such that λ T;i < μ i λ N;i where μ i is the coefficient of friction and can either be a constant or a function of the relative velocity. The contact can be in one of three modes: stick, trans or slip mode. In stick mode, the friction force λ T;i is any value within the bound such that the relative acceleration along the surface of the contact point is constrained to zero, that is, the bilateral constraint € g T;i ¼ 0 is enforced. In slip mode, the friction force is on the bound, acting in a direction opposing the velocity. The transitional-slip mode is the same as the slip mode but where the friction force acts in a direction opposing acceleration rather than velocity. Hence, the friction force is determined by, Impact energy transitions can be modelled using Newton's restitution law [11]: i are the post-impact velocities and _ g À N;i ; _ g À T;i are the pre-impact velocities. The terms e N;i ; e T;i are the normal and tangential restitution coefficients, respectively.
We can create complementarity conditions similar to the contact force laws. Define Λ N;i as an impulse along the surface normal to the ith contact. Then, based on the non-penetration constraint, we have & and in the tangential plane, ( This impact law has, in general, poor prediction qualities when comparing with experimental observation. However, the law is simple and is very numerically and computationally tractable [12]. In this paper, we will use the lower case λ interchangeably to refer to both contact forces and impulses (Λ). It will be obvious to the reader which we are referring to, force or impulse, from the context of its use and the supporting discussion. By choosing λ to represent both force and impulse, we remove redundant notation and avoid the unnecessary duplication of some of the equations we are about to introduce.

The MRB hybrid automaton
In this section, we will describe the components of our modelling framework, mainly the elements of the MRB hybrid automaton, the contact list and the contact combination graph. These latter two components are essential for the definition of the possible configurations and contacts that may occur in the system. The elements of our hybrid-automaton-based framework and their interrelationships that will be described in the following subsections can be interpreted as an event-driven scheme for the numerical time-integration of a type of mechanical systems with multiple impacts, contacts and discontinuous friction. Here, we invite the reader to check well-known event-driven schemes that have been specifically designed for the numerical time-integration of non-smooth mechanical systems, especially Chapter 8 of Ref. [5].
The MRB hybrid automaton is based on the hybrid model given in Ref. [3,4] and the DDSs hybrid automaton proposed in Ref. [24]. For the sake of simplifying the notation, we will not consider inputs and outputs for the basic MRB hybrid automaton.

The general MRB hybrid automaton
The MRB hybrid automaton has the following general definition. The details of its elements are explained in the subsequent sections.
Consider a system consisting of M entities with P possible contacts. The MRB hybrid automaton is a collection, with: • S a finite set of discrete locations. We will have two types of discrete locations: dynamical discrete locations (the classical discrete locations in a hybrid automaton, that is, locations with an associated dynamical system) and non-dynamical discrete locationsthe new type of discrete locations introduced in our modelling frameworkwhich are grouped into computation nodes. The reader is referred to Sections 3.2, 3.3 and 3.4 for more details. The set of dynamical discrete locations is denoted by S d , and the set of non-dynamical discrete locations is denoted by S non , and S ¼ S d S S non . • X R 12Mþ3P is the continuous state space. The continuous state vector is a generalized coordinate vector q 2 R 6M and a generalized velocity vector _ q 2 R 6M for all the M entities in the system, plus λ N ¼ λ N;1 ; λ N;2 ; . . . ; λ N;P È É 2 R P and λ T ¼ λ T T;1 ; λ T T;2 ; . . . ; λ T T;P n o 2 R 2P as the continuous vectors of contact forces, acting at each of the P possible contacts. That is, for x 2 X; x ¼ q; _ q; λ N ; λ T ð Þ T . • E S Â S represents the set of discrete transitions or edges in H MRB , which is finite. The discrete transitions are deterministic and are of different types, being between dynamical discrete locations, non-dynamical discrete locations and combinations of both; they are specified in Section 3.5. • Dom : S d ! 2 X are the dynamical discrete location domains. Dom assigns a set of continuous states to each dynamical discrete location of S d , thus, for s 2 S d , DomðsÞ X. These domains are formally specified in Section 3.6. • F represents the continuous dynamics: F ¼ f s ðxÞ : s 2 S d f gis a collection of vector fields such that f s : X ! X describes the dynamics in each dynamical discrete location. Each f s ðxÞ with x 2 X is Lipschitz continuous on X in order to ensure that in each dynamical discrete location, s the solution exists and is unique. • Init S d Â X is a set of initial states. • G represents the guard maps for each discrete transition. G : E ! 2 X . G assigns to each edge e 2 E a set of continuous states (GðeÞ & X) which enables transitions along that edge. The guards are of different types depending on the type of discrete transition they are associated with, being between dynamical discrete locations, non-dynamical discrete locations and combinations of both. More details are given in Section 3.5. • R denotes the reset maps, which define the re-initialization of the continuous state every time a discrete transition takes place. All the resets associated with discrete transitions between dynamical discrete locations do not change the continuous state. The resets between non-dynamical discrete locations within computation nodes depend on the numerical procedure which uses the dynamical discrete location models. More details are given in Section 3.7.
The formal definitions of a hybrid time trajectory and an execution of the MRB hybrid automaton on such a trajectory are similar to the ones given for classical hybrid automata and can be found in Ref. [3]. In general terms, the consideration of computation nodes does not change the definition of the hybrid time trajectory and the execution of the hybrid automaton, since time does not evolve in the computation nodes. However, these nodes are considered as nodes that perform some necessary computations, and roughly speaking, 'they guide the executions' of the hybrid automaton to the relevant dynamical discrete locations based on the information of the contact list and the contact combination graph (described in the next section), and consequently, they might be considered as part of the dynamical discrete locations of the MRB hybrid automaton. The influence of the computation nodes on the evolution of the dynamics is formally specified in the sections below.

Contact list and contact combination graph
In this section, we will define the contact list and the contact combination graph, which are necessary concepts to later define the dynamical and non-dynamical discrete locations in Sections 3.3 and 3.4.
The set of all gap functions g N;i for every possible contact that could occur in the system is . . .  We assign to each gap function in the set Γ N an index or label and define " I as the collection of all these labels. Thus, i 2 " I , g N;i 2 Γ N . Note that the set " I naturally labels the elements of Γ N , _ Γ N and € Γ N . That is, for some i 2 " I, we have g N;i 2 Γ N , _ g N;i 2 _ Γ N and € g N;i 2 € Γ N . As an example, take the three-ball problem in Figure 6. There are three possible contacts and thus three gap functions: g N;ab between balls a and b, g N;bc between balls b and c and g N;ac between balls a and c. Thus, " I ¼ ab; bc; ac f gand Γ N ¼ g N;ab ; g N;bc ; g N;ac È É . The contact combination graph is a structure that contains information on the possible configurations the system of entities could take, and whether it is possible for the system to evolve from one configuration to another. Each vertex of the contact combination graph represents some unique combination of closed and open contacts which could possibly occur at some time in our system. A closed contact is when two surfaces are touchingthey are actually in contactand an open contact is when they are not. Each edge represents a possible transition from one set of closed and open contacts to another. Then, the contact combination graph is denoted by Ω ¼ V; E Ω ð Þ, where V is a set of vertices and E Ω a set of edges. For each vertex V k 2 V, we define I k 2 " I as the set of indices of all closed contacts when the contact combination graph is in the kth vertex, with q as defined in Section 2. Each set of indices will form a configuration of the system. If it is possible for the system to directly evolve from configuration I k , to a configuration I l , then there is an edge between V k and V l , that is k; l ð Þ 2 E Ω . From here on, we will refer to any specific contact combination by its index. That is, when the contact combination graph is in the kth vertex, we say that the system is in the kth contact combination. Further, as a convention, we will reserve k ¼ 0 for the contact combination where no contacts occur I k ¼ ; f g ð Þ . Finally, we denote N k ¼ I k j j as the number of contacts occurring in the kth contact combination.
Returning to the three-ball example, let us say that V 1 is associated with the configuration I 1 ¼ ab; ac f g, shown on the left of Figure 7, V 2 is associated with I 2 ¼ ab; ac; bc f g , shown on the upper right, and V 3 with I 3 ¼ ab; bc f g, shown on the lower right. There is an edge between V 1 and V 2 because it is possible for the three balls to change from one configuration to the other. For the same reason, there is an edge between V 2 and V 3 . However, the three balls cannot change directly from V 1 to V 2 without being in some intermediate configuration; therefore, there is no edge between these two.

Dynamical discrete locations and associated dynamics
For each vertex V k 2 V, there is a set of 3 N k dynamical discrete locations: S k ¼ s k;1 ; s k;2 ; . . . ; s k;3 N k È É . Each dynamical discrete location represents a different arrangement of stick, slip or transition mode at each of the contact points. For example, the vertex V 1 in our three-ball example has two contacts. These contacts could both be in stick mode, or both in slip mode, or one in trans mode, the other in stick mode, and so on. Each combination (stick-stick, slip-stick etc.) is represented by a different dynamical discrete location and the ensuing dynamics are governed by a different continuous dynamical equation.
To each location s k;i , we associate the following index sets: • I st;k;i the indices of all contacts where friction is in stick mode: S k . We also note that I st;k;i [ I tr;k;i [ I sl;k;i ¼ I k and I st;k;i \ I tr;k;i \ I sl;k;i ¼ 1 f g. For each dynamical discrete location s k;i , we have a continuous dynamical system of the form: where MðqÞ is the mass matrix, hðq; _ qÞ is a force vector containing all conservative forces, input forces, and Coriolis terms, and The contact forces λ N ; λ T are computed using the following equation, The þ indicates the pseudo-inverse. If there are no redundant constraints, then the matrix on the right-hand side becomes non-singular and the pseudo-inverse can be replaced by a conventional matrix inversion.

Non-dynamical discrete locations and computation nodes
The MRB hybrid automaton includes a new type of discrete locations, referred to as non-dynamical discrete locations. The collection of these discrete locations, and the edges, guards and resets between them, is termed a computation node, whose main elements are given in Table 1. The discrete locations in the computation nodes do not have any continuous-time dynamics.
The computation nodes are used to model the numerical methods that affect the dynamics of the system, but which cannot themselves be represented as dynamical elements. To clarify, in our MRB hybrid automaton, we use the computation nodes to model the numerical procedures used in calculating contact forces at impact, or the state of the contact forces (stick, slip etc.) when in sustained contact. In the framework of rigid-body dynamics, we often assume that changes in contact forces are instantaneous. Otherwise, we must model the fast-acting microscopic dynamics, which is computationally wasteful and largely unnecessary in a macroscopic analysis. Thus instead, we use iterative or LCP-based numerical procedures to calculate the contact forces. LCP stands for linear complementarity problem.
The computation nodes are inspired by the well-known mythical modes, originally proposed in Ref. [25] and later expanded to hybrid systems by Ref. [26,27]. They have to be understood in this sense, just a set of transition intermediate states associated to a discontinuous change. They do not have any continuous dynamics. Further, one might consider the physical analogue of the computation nodes as being the fast-acting dynamics operating at the microscopic scale which we conventionally assume to act instantaneously in the framework of rigid-body mechanics.
For each vertex V k 2 V k 2 V : N k > 0 f g , there is a set of non-dynamical locations: I k ¼ I k;en ; I k;1 ; . . . ; I k;ex È É ; C k ¼ C k;en ; C k;1 ; . . . ; C k;ex È É ; Table 1. Main elements of the computation nodes.
Purpose Elements Types To compute contact forces Set of non-dynamical discrete locations Entrance discrete location (all edges going into the computation node end here) Exit discrete location (all edges to locations outside the node leave from here) Edges and guards between non-dynamical discrete locations Reset functions to find contact forces on transitions Impact computation nodes Contact computation nodes which are called the impact computation node and the contact computation node, respectively. Both computation nodes have similar form. The sets I k and C k consist of an entrance discrete location with subscript en, an exit discrete location with subscript ex and a set of discrete locations in between which are used for the computation of the contact forces. The number of these discrete locations varies in number depending on the algorithm used to calculate the contact forces. The edges going into the computation node are received by the entrance discrete location, and all edges leaving the computation node leave from the exit discrete location. The edges connecting the intermediate locations are explained in Section 4. We highlight that every different computation node will have a different number of non-dynamical discrete locations.
The algorithm we use to compute the contact forces is the successive over-relaxation proximal point method (SORPROX) as described in Ref. [14]. The SORPROX method operates by computing some new estimate of the contact force at some point of contact based on some previous initial guesses. If it is a feasible contact force, that is, it satisfies the complementarity conditions in Section 2, then the new estimate is accepted and becomes the next iteration of the contact force. Otherwise, the next iteration is chosen to be the feasible point nearest the new estimate. A new estimate is then produced for the next contact force, and so on, repeatedly cycling through and iterating the contact forces until they converge to a solution.
An outline of the overall computation node is shown in Figure 8. For simplicity, let us say we are in the ith impact node with the contact index set I i ¼ 1; 2; 3; . . . ; n f g . Recall from our earlier definition that the contact index is an index of all closed contacts in a given configuration. In this case, there are n contact points, indexed 1; 2; 3; . . . ; n.
When an impact occurs, the entrance location I i;en becomes active. The execution of the computation node traverses each of the intermediate locations, computing new contact forces, until the nth node I i;1 ! I i;2 ! I i;3 ! Á Á Á ! I i;n . The location I i;n has an edge going to the exit node and an edge going to the entrance node. If the all the constraints are satisfied, then a transition I i;n ! I k;ex occurs, and the state vector q is reset. If the constraints are not satisfied, the transition I i;n ! I i;en takes place and the execution of the computation node is repeated to determine a new set of contact forces. In this paper, we use computation nodes based on the SORPROX method [14] for calculating contact forces. We have previously reported details of the computation node for SORPROX in Ref. [23]. For the sake of brevity, we will omit repeating these details in this paper and refer the interested reader to the referenced work.

Edges and guards
The set of edges of H MRB consists of the union of the following: • For all k such that V k 2 V : N k > 0 f g : • For all k and j such that ðk; jÞ 2 E Ω and N k > N j^Nj > 0 À Á : • For all k and i such that V k 2 V : N k > 0^I tr;k;i ¼ 0 È É : • For all k, j and i such that V k 2 V : N k > 0^I tr;k;i > 0^I tr;k;j ¼ 0^I tr;k;i [ I sl;k;i ¼ I sl;k;j À Á È É : s k;i ; s k;j À Á ; . . . È É : • For all k, j and i such that ðk; jÞ 2 E Ω , N k < N j À Á and s k;i 2 S k : s k;i ; I j;en À Á ; . . . È É : • For all k, j and i such that ðk; jÞ 2 E Ω , N k > N j^Nj > 0 À Á and s k;i 2 S k : s k;i ; C j;en À Á ; . . . È É : • If a discrete location S 0 such that N 0 ¼ 0 exists, then for all ðk; 0Þ 2 E Ω f g : and for all i such that s k;i 2 S k : The guards for each of the edges in H MRB are as follows: G I k;ex ; C j;en À Á ¼ x 2 X : "n 2 I j ; _ g N;n ¼ 0 "n 2 I k nI j ; _ g N;n > 0 n o ; G I k;ex ; C k;en À Á ¼ x 2 X : "n 2 I k ; _ g N;n ¼ 0 G s k;i ; s k;j À Á ¼ x 2 X : "n 2 I tr;k;i ; _ g T;n > 0 n "n 2 Å k nI k ; g N;n > 0 _ _ g N;n > 0 "n 2 I k ; λ N;n ! 0 À Á É ; G s k;i ; I j;en À Á ¼ x 2 X : "n 2 I j ; g N;n 0^_ g N;n 0 n "n 2 Å k nI j ; g N;n > 0 _ _ g N;n > 0 o ; where λ N;j ; λ T;j are jth rows of the vectors λ N ; λ T given in equation (15). Finally, the following guards are specific only to S 0 , These formula for the edges and guards can be better understood by applying them to specific examples. We refer the reader to Section 3.8 where some examples are presented, and to check the programs that form DyverseRBT and the input files for the examples used in this paper, which are available at http://staff.cs.manchester.ac.uk/navarroe/research/dyverse/DyverseRBT_BMC.zip.

Location domains of dynamical discrete locations
The location domains Dom are only considered for dynamical discrete locations, s k;i 2 S d : "j 2 I sl;k;i _ g T;j > 0 o :

Resets
We highlight that all the resets associated with transitions between dynamical discrete locations do not change the continuous state. The resets between non-dynamical locations within computation nodes depend on the numerical procedure which uses the dynamical discrete location models, and we refer to this type of resets as non-general resets, because they cannot be described with a general expression. For brevity, we do not include a description of the computation node we use in the present work. However, it is straightforward to derive itincluding their resetsfrom the SORPROX method. We refer the reader to Ref. [23] for more details.
For the sake of clarity and to ease the specification of the model, we are taking some liberties with notation and do not employ the typical nomenclature of resets in hybrid automata.
Although we suggest that the resets in the computation nodes are non-general, we make two exceptions within the MRB hybrid automaton for which we give an expression for the resets. For an impact computation node consisting of the collection of the non-dynamical locations I k;en ; I k;1 ; . . . ; I k;N ; I k;ex È É , there is a reset on the edge I k;N ! I k;ex , Similarly, for a contact computation node consisting of locations C k;en ; C k;1 ; . . . ; C k;N ; C k;ex È É , there is a reset on the edge C k;N ! C k;ex , 3.8. The interception game example and multiple contact problems 3.8.1. MRB hybrid automaton for the interception game Now, we return to our motivating example introduced in Section 1.3. The MRB hybrid automaton for the interception game is shown in Figure 9. The guards, resets, and vector fields are omitted for clarity. The automaton features four dynamical locations: S 0 for the free-state (i.e. the balls are not in contact with one another), s 1;1 when the balls are in contact and in slip mode (i.e. they are Figure 9. Hybrid automaton of the interception game. slipping across each other), s 1;2 for contact and transition mode (or trans mode, see Section 3.3), and s 1;3 for contact and stick mode. The impact and contact computation nodes are shown with their non-dynamical locations. The computation nodes are shown in the form for resolving contact forces using the SORPROX method. The first transition I 1;en ! I 1;1 (or equivalently C 1;en ! C 1;1 ) sets an initial guess for the two contact forcesthen normal contact force, and the tangential friction force. The remaining transitions make new approximations of the contact forces and then test to see whether they have converged. If they have not converged, the cycle is repeated (I 1;5 ! I 1;1 or C 1;5 ! C 1;1 ); otherwise, the appropriate resets are applied on the state vector and the computation node is exited.
An MRB hybrid automaton is constructed for the interception game example using the DyverseRBT tool. The input file given by the user contains two entities (ball 1 and ball 2), each with a state vector, a chosen mass and a surface function for a sphere. This file was displayed in Figure 4. The input forces u1, u2 and u3 are external forces which can be input manually.
The implementation generates the contact list, the contact combination graph, locations, and computation nodes of an MRB hybrid automaton for the interception game model. It also generates the vector fields and domains of each location, and the edges with their respective guards and resets.
The hybrid automaton is then converted into C code for an S-function which is compiled using the on-board MEX compiler. The S-function can be dropped into a Simulink model and used as an ordinary block. The external inputs (u1, u2 and u3) given in the text file of Figure 4 become the inputs to the S-function block in the Simulink model. All this process is automated and the result is shown in Figure 10.

Multiple contact problems
We also consider the automatic construction of the MRB hybrid automaton from rudimentary rigid-body models consisting of 1, 2, 3 and 4 balls using the DyverseRBT tool. The case of 3 balls was considered in Section 3.2. The tool first constructs a list of possible contact situations that occur in the model. Then from this, it derives the number of the different types of locations (contact, impact and dynamical) and assigns the possible edges to each of the locations. Finally, the tool derives the dynamic equations governing the motion for each location from the contact situation assigned to that locationfor example, whether the contact is open or closed, stick, trans or slip etc. Finally, the guards and resets are determined. The end result is a complete and precise description of the rigid-body system, as a Matlab structure, in the formulation of an MRB hybrid automaton. Table 2 shows the number of locations, edges and possible contacts determined by DyverseRBT for each model. Clearly, the number of locations and number of edges between locations rapidly increase with the number of possible contacts. This combinatorial explosion limits the practical size of any model which can be completely described in this way, even with relatively simple rigidbody entities.

Application of the MRB hybrid automaton to falsification of safety properties
In this section, we briefly describe DyverseBMC: a procedure for falsifying a safety property over finite time intervals for a rigid-body system against our hybrid automaton abstraction. DyverseBMC has to be understood as a valuable application of the wider modelling and simulation framework that we propose in this paper, which includes the specification of the computational semantics for an MRB system. This specification is given by the MRB hybrid automaton formalism and was explained in the previous section. The automatic generation of the MRB hybrid automaton and its simulation are implemented through the tool DyverseRBT. The main elements of our framework were given in Figure 1.

Outline of the falsification procedure
The problem to solve is formally specified as • For a multi-rigid body mechanical system H MRB with continuous states xðtÞ 2 X and physical parameters p 2 R n p , the property ϕðxðtÞÞ must hold for all xðtÞ such that xð0Þ 2 I R n , p 2 P & R n p and t 2 T & R, with n; n p 2 N some natural numbers, with ϕðxðtÞÞ a safety property, I a set of initial conditions, P a set of possible values for the physical parameters (mass, damping etc.), and T some time interval. It is well worth reminding that the safety problem will be reduced to the satisfiability of Boolean formulas, and the approach to follow will be that of falsification.
Verification of hybrid systems is generally a very difficult problem. There are a number of practical restrictions to our method in its current form and in its current implementation: • Rigid-body dynamical equations should have low stiffness. The method of falsification contains a numerical integration component to evaluate ODEs. The consequence is a buildup of error in the continuous states. In order to keep this error low and to ensure that numerical integration method is stable and efficient, it is necessary that the rigid-body system has relatively low-stiffness dynamics. • Solution uncertainty and conservative safety properties. The numerical integration error is also significant when finding state-space trajectories which violate the safety property. It is uncertain whether a trajectory which runs very close to the boundary of the region in the state space defined by the safety property is safe or unsafe. The integration error may cause, or prevent, the trajectory from crossing the boundary resulting in an incorrect result. To avoid this problem, it is necessary to choose conservative safety properties such that the numerical error will never be large enough over the bounded time interval to cause the bounded model checker to imply safety when the system is unsafe. This still leaves the opposite possibility that the system model checker declares unsafe when the system is safe. However, this can be checked by analysing the counter example and by successive refinement of the safety property. However, the process of successive refinement may never terminate, something well proven even for linear continuous-time systems and linear hybrid systems [57,58,62]. This problem might be solved with the use of alternative numerical time-integration techniques: event-capturing time-stepping schemes (for instance, Ref. [34]), which are not based on the accurate detection of events and are quite robust in contact detection, unlike our event-driven scheme. Note that the work in Ref. [34] proposes a scheme which stabilizes the constraints, which implies that the constraint drift phenomenona well-known phenomenon in the systems treatedis avoided. • Short-time intervals and small number of location transitions. The error in the continuous states caused by the numerical integration increases with time. Although we have not characterize this error, we can conclude that the bounded model checker should only be used to check short intervals to keep this error small. We also restrict to small number of transitions from location to location in the MRB hybrid automaton. For example, transitions caused by impacts or changes in friction (stick to slip, for example). The complexity of the problem rises significantly with each transition. For practical solving times, we therefore only consider one or two transitions at the most.
Our falsification approach is based on a methodology akin to BMC described in Section 1. A logical formula is constructed consisting of Boolean variables and constraints on real variables. This formula is satisfied only by a certain set of trajectories across the continuous state space over finite time intervals. A lazy SMT solver [74,79] is used to find a solution to this formula which falsifies the safety property. If no falsifying solution exists, then we can conclude that there are no trajectories amongst that set which are unsafe over the time interval. We then move to a different set of trajectories and repeat the search for unsafe trajectories.
The overall procedure is shown in Algorithm 1. The algorithm uses a depth-first-search method to build possible discrete evolutions, which are then checked for safety. More specifically, the discrete evolution is checked for the existence of a continuous trajectory that starts in the initial conditions, visits each discrete location in order and, at some point, violates the safety property.
For instance, the discrete evolution s 1 ! s 2 ! s 3 is unsafe, if there exists a continuous trajectory which starts from a set of initial conditions in the domain of the discrete location s 1 , then enters the domain of s 2 , followed by the domain of s 3 , and at some stage enters a region of the state space specified as unsafe.
If the discrete evolution is determined to be safe, then a new location is added to the evolution following the depth-first-search method. This new and longer discrete evolution is then checked for safety, and so on, until an unsafe solution is found, or until some chosen maximum depth is reached.
The safety verification task is conducted using the function ExploreDynamicalLocation. The functions CreateRoot and ExploreCompNode are used to aid in constructing formulas for trajectories that traverse between discrete locations and domains. The purpose and details of these functions are given in the 'Supplemental online material' file. In this file, we also explain the main limitations of our method. One of the most significant limitations is the need to use numerical integration in the Lazy-SMT solver to compute ODE functions.
As we pointed out in Section 1, even though DyverseBMC is not optimized and has limitations, still it represents a valuable application of the proposed alternative modelling framework for the type of systems the MRB hybrid automaton is modelling.

Convex and concave constraints
It is important to give a note on the convex and concave constraints of our falsification procedure. The constraints which appear in the clauses of the formula in our method (see the 'Supplemental online material' file) are not restricted to be convex. In the case of our case study in the next section, both convex and concave constraint surfaces appear regularly in the formulas the Lazy SMT solver is required to solve. This represents a problem. If the Lazy SMT solver returns that the formula is infeasible, how are we to know that this truly is the case? It is possible that the minimization may have become trapped in a local but infeasible minimum, where a feasible minimum also exists. In which case, our determination of infeasibility is dependent on the initial conditions and the result is unsound. Global optimization methods are often statistically based and thus can only assert probable infeasibility.
To solve this problem, we are fortunate in the definition of our MRB hybrid automaton. In Section 2, we imposed the restriction that the contact surfaces are convex. Consequently, as the constraint surfaces in our formula are derived from these contact surfaces, we limit the types of constraint surfaces to either a convex surface or a concave surfacenot both in different regions.
For convex surfaces, infeasibility can be known, and the result is sound. For concave surfaces, it follows that the minimal points are on the edge of the allowable range of the variables in our minimization problem. Thus, we in effect know how many local minima exist, and approximately where they are, and can cycle through each one by choosing appropriate initial conditions until one becomes feasible. If none are feasible, then we can return the result infeasiblewith the minimum number of infeasible constraintswhich we now know to be sound.

Implementation
DyverseBMC is implemented in Matlab (Mathsworks Inc., USA). Boolean satisfiability problems are solved using Matlab's on-board SAT solver. Numerical integration of ODE problems is performed by Matlab's ode45 function. Optimization problems are solved using the NAG toolbox nag opt nlp2 solve function for minimizing arbitrary smooth objective functions subject to smooth linear or non-linear constraints [87].
The programs that form DyverseBMC and the input files for the examples used in this paper are available at http://staff.cs.manchester.ac.uk/navarroe/research/dyverse/DyverseRBT_BMC.zip.

The interception game revisited
Recall our motivating example introduced in Section 1.3. The aim is to prove that, for all K over an interval K 2 1; 3 ½ , where K is the feedback controller gain of ball 2, the game will not be lost within the first two seconds of play. Thus far, we have: • An MRB hybrid automaton of the systemsee Section 3.8.1. • A set of initial conditions, I, as defined in Section 1.3. • A property which must hold, y 1 ðtÞ > 0; where x 1 ; x 2 are the state vectors for ball 1 and ball 2, respectively, as defined in Section 1.3, t is time in seconds and the dot denotes derivative with respect to time. By checking this property, we find whether the game can be sustained for the first 2 s and, by counter example, the solution of K if the game is lost.
The MRB hybrid automaton is constructed for the interception game using the DyverseRBT tool and checked to falsify the safety property using the DyverseBMC tool. The external forces input in the Simulink model (see Figure 10) are added to the automaton manually in the form of symbolic expressions.
We look at how the size of the problem changes as the model checker checks for safety further into the future of the model, and how this change in problem size can affect the performance. Note that the current implementation has not been optimized for performance, and we therefore offer no benchmarks with other SMT-type packages. Furthermore, other SMT-type packages do not support numerical integration which is necessary to generate trajectories. We focus instead on how the method of posing the safety property and physical model in the form described in this paper will affect the SMT solver.
The simulation time is the length of the possible trajectories in time, or in other words, the amount of time that has been checked in the model of the interception game. Figure 11 compares the simulation time with the number of calls to the Lazy SMT solver. We point out that the Lazy SMT solver is invoked by the function Satisfy, in Algorithms ExploreDynamicalLocation and ExploreCompNode given in the 'Supplemental online material' file. The CreateRoots function adds to the formula a possible range for the simulation time. This divergence of simulation time is represented by a dark-grey region. Each edge transition invokes a call to CreateRoots, which diverges this range further. The plateau in the light-grey region is caused by calls to satisfy while exploring the computation node, i.e. Satisfy is called from Algorithm ExploreCompNode. Computation nodes are considered to act over instantaneous moments in time and, therefore, do not advance the simulation time.
Each iteration of the outer-most while-loop in Algorithm ExploreDynamicalLocation increments the simulation time by a time step. During an iteration, the Lazy SMT solver is asked to solve at least N þ 2 formulas, where N is the number of edges of the dynamical location. Specifically, the Lazy SMT solver is called N times to check if a trajectory exists which enters any of the N guard sets (it is called once for each of the N guards), then once to check if a trajectory is in the domain, and if it is, again to check for unsafe trajectories.
When exploring a computation node, the Lazy SMT solver is called at each iteration of the contact force to test if the new iteration has any non-convergent solutions. After all possible contact-forces have converged, the solver is called a further N times to check for solutions that enter a guard, where as before N is the number of edges leaving the computation node.
The procedures of DyverseBMC are premised on adding new clauses to a formula as new time steps are added and contact situations change. Figure 12 shows the size of the formulas given to the Lazy SMT solver with respect to the simulation time. The size is judged in terms of the number of clauses with one and two literals in the formulas. The number of clauses is broadly stable while exploring dynamical locations. The sudden increase occurs when the computation node is being explored. Thereafter, when exploring the dynamical location, the size of the formulas remains stable. It is primarily the occurrence of contact in the model which increases the difficulty of the problem. The difficulty of finding a solution to a given formula grows as the number of clauses increases. The number of single-literal clauses is less significant than clauses with multiple literals. Single-literal clauses only add new constraints for the constraint solver. However, clauses with multiple literals (e.g. a _ b) increase the number of possible Boolean abstractions which might satisfy the formula. Figure 13 shows the number of iterations in Algorithm Satisfy (see the 'Supplemental online material' file) required to solve (or prove infeasible) a formula with respect to the number of one-literal and two-literal clauses in that formula. Evidently, the number of iterations increases more dramatically with the number of two-literal clauses compared to single-literal clauses.
In the 'Supplemental online material' file, we also include a multiple contact problem to explore the limitations of DyverseBMC.

Conclusions
In this paper, we have presented a novel computational model referred to as the MRB hybrid automaton to fully describe the transitions and operation modes present in a general class of mechanical systems with impacts and friction. Significantly, our computational model can accommodate multiple contacts. This extends the already existing results related to the modelling of systems with impacts by using a class of hybrid automata referred to as MRB hybrid automata. One of the chief characteristics of our model is the inclusion of computation nodes to calculate the contact forces. Each computation node consists of a set of non-dynamical discrete locations, discrete transitions and guards between these locations, and resets on transitions, which can account for the energy transfer not explicitly considered within the rigid-body formalism. Based on the MRB hybrid automaton formalism, we propose a specification of the computational semantics for a class of MRB systems. The automatic generation of the MRB hybrid automaton and its simulation has been implemented through the tool DyverseRBT. The proposed MRB hybrid-automaton-based modelling framework is well suited for the automatic generation of simulation models and the formal verification of dynamical properties of realistic mechanical systems. Due to the complexity involved, building these models or verifying these systems manually is laborious and could be eased by automatic computational techniques as the ones presented in our work. Note that we assume here that j is incremented from 1 to P and that the normal forces are computed before the tangential ones. The contact forces and impulses must be part of the feasible set implied by the complementarity conditions defined in Section 2. A proximal point projection function is used to enforce this. The function returns the closest point in the feasible set and assigns it to our current iteration of the contact forces. Using the specific contact laws in this paper, we have for all P contacts.