CompileTime Extensions to Hybrid ODEs
Abstract
Reachability analysis for hybrid systems is an active area of development and has resulted in many promising prototype tools. Most of these tools allow users to express hybrid system as automata with a set of ordinary differential equations (ODEs) associated with each state, as well as rules for transitions between states. Significant effort goes into developing and verifying and correctly implementing those tools. As such, it is desirable to expand the scope of applicability tools of such as far as possible. With this goal, we show how compiletime transformations can be used to extend the basic hybrid ODE formalism traditionally supported in hybrid reachability tools such as SpaceEx or Flow*. The extension supports certain types of partial derivatives and equational constraints. These extensions allow users to express, among other things, the EulerLagrangian equation, and to capture practically relevant constraints that arise naturally in mechanical systems. Achieving this level of expressiveness requires using a binding timeanalysis (BTA), program differentiation, symbolic Gaussian elimination, and abstract interpretation using interval analysis. Except for BTA, the other components are either readily available or can be easily added to most reachability tools. The paper therefore focuses on presenting both the declarative and algorithmic specifications for the BTA phase, and establishes the soundness of the algorithmic specifications with respect to the declarative one.
1 Introduction
Reachability analysis for hybrid systems [Alur93] is an active area of development and has resulted in many promising prototype tools. Prominent examples of such tools include CHARON [CHARON], HyTech [HyTech], PHAVer [phaver], dReach [dreach], dReal [dreal], SpaceEx [spaceex], and Flow*[flowstar]. Most of these tools allow users to express hybrid systems as automata with a set of ordinary differential equations (ODEs) associated with each state, as well as rules for transitions between states. In particular, ODEs must be in the explicit form where the left hand side of an equality has to be the derivative of a state variable. Significant effort goes into verifying and correctly implementing those tools. As such, it is desirable to expand the scope of applicability tools of such as far as possible.
1.1 Contributions
With this goal, we present a systematic method to translate an expressive language with partial derivatives and equations to a standard language supporting ODEs, guards, and reset maps. The method can be used to extend reachability analysis tool such as SpaceEx or Flow*. An experimental implementation of the proposed technique is available in the freely available, open source Acumen language implementation [acumenURL]. Examples illustrating the use of these extension can be found in the directory examples/04_Experimental/04_BTA
. Since both partial derivatives and equations are eliminated completely after the compiletime transformation, the user benefits from the added expressivity but the underlying tools do not need to change. The two extensions allow the user to express, among other things, the EulerLagrangian equation, and to capture practically relevant constraints that arise naturally in mechanical systems. Achieving this level of expressivity requires using a binding timeanalysis (BTA) [Jones85, Gomard91, christensenaccurate, Moggi97], program differentiation, symbolic Gaussian elimination, and abstract interpretation using interval analysis. Except for BTA, the other components are either readily available or can be easily added. The technical part of the paper therefore focuses on presenting both the declarative and algorithmic specifications for the BTA phase, and establishes the soundness of the algorithmic with respect to the declarative.
After reviewing related work on compiletime extensions (Section 2), we introduce the syntax and type system for a core differential equation language (Section 3). Then, we present a declarative specification of bindingtime analysis (BTA) and a big step semantics for specialization (Section 4), along with a formal proof of type safety (Theorem 1). We then present an algorithmic specification of the BTA that works by first generating a set of constraints and then attempting to solve them (Section LABEL:section:Implementation), and we show that this algorithmic specification is faithful to the declarative BTA (Lemma LABEL:lemma:faithful) and always produces a unique minimum solution that maps as much of the code as possible to static if an assignment exists (Theorem 2). To illustrate the practical value of the formalism, we present two case studies that have been carried out using the implementation (Appendix LABEL:appendix).
2 Related Work
Bindingtime Analysis (BTA) is a static analysis traditionally supported in the offline partial evaluation of general purposes languages. It works by identifying a twolevel structure in the program being analyzed, where the first level is a computation that can be done at “partial evaluation time” (“compile time” in our case), and the second level must be left as a “residual” that is executed at runtime. BTA has generally been studied for general purpose languages. In our setting, we study it in the context of Domain Specific Languages (DSLs) [Hudak97, Mernik05, Sujeeth14] intended for modeling hybrid systems. It should also be noted that our primary purpose is to use it for extending expressivity. Partial evaluation, in contrast, is only concerned with improving the runtime performance of programs. In what follows, we elaborate on these key points.
A key idea in the work we present in this paper is that there are powerful techniques from the programming languages community that can help make reachability tools more broadly applicable. To do this, this work uses twolevel languages in a novel way. To put the existing related work in context, it is useful to consider several characteristics relating to the language considered and the transformation used, namely, whether the language is domainspecific (or general purpose), whether it supports equations, whether the transformation is done at compiletime (or runtime), whether the tool performs letinsertion (to avoid code duplication), whether the language is statically typed, and whether the tool provide accurate source level error reporting. The systems that we will consider are partial evaluation systems for C, namely, Cmix [CMix] and Tempo [tempo]; template instantiation mechanisms, namely, C++ Templates [metaOCamel] and Template Haskell [templateHaskell]; multistage programming languages, namely, MetaOCaml [metaOCamel] and LMS [LMS]; the Verilog Preprocessor [VPP]; and the hybridization technique [bak].
Static  
checking  Source level  
reporting  Compiletime  Domainspecific  Let  
insertion  Equations  
CMix  Yes    Yes  No  Yes  No 
Tempo  Yes  Yes  Yes  No  Yes  No 
C++ Template  No  No  Yes  No    No 
Template Haskell  Yes  Yes  Yes  No    No 
MetaOcaml  Yes  Yes  No  No  No  No 
LMS  Yes    No  No  Yes  No 
Verilog preprocessor  Yes  Yes  Yes  Yes    No 
Hybridization  Yes    Yes  Yes    No 
This paper  Yes  Yes  Yes  Yes  Yes  Yes 

Table 1 provides an overview of how these different systems related to these key properties. The main observations from the table are as follows. Almost all tools are compiletime (except MetaOCaml), and almost all are statically checked (except C++ Templates). A key feature of static checking is that it facilitates accurate sourcelevel reporting. That is primary reason for choosing an approach based on BTA or some type of static analysis. Compiletime program specializers, such as CMix and Tempo focus on automatically specializing a program through a well understood set of transformations to produce a program that is faster than the original one. There are no fundamental reasons why specialization (and twolevel languages) need to be limited to general purpose languages. In fact, as this paper shows, they can be quite useful as they can be used to increase expressivity. Letinsertion was invented in the partial evaluation community, and is adopted by automated tools by LMS (but not be explicit tools like MetaOCaml). It is quite critical when there are significant compiletime computations, as is the case when we are trying to eliminate nontrivial constructs like partial derivatives or performing substantial manipulations to turn equations into formulae. However, none of these works address the question of supporting equations, that is, allowing the user to write constraints in equational form, and then translating them directly into “formula” form for directed evaluation.
Syntax  
Constant  
Variable  
Type  
Type Environment  
Function  
Expression  
Equation  
Our work is comparable to that of HyST [hyst], which is a tool that aims to facilitate interchange of models between different tools. This way, HyST facilitates sharing of models and comparing solvers. In contrast, our work explores another dimension for reuse, namely, how these tools can be extended to support a more flexible and expressive modeling formalism.
3 A Differential Equation Language (Acumen17)
Fig. 1 introduces the syntax and type system for a core differential equation language called Acumen17. We use the following notational conventions:

Writing denotes a vector . We will occasionally omit the superscript and write when the range of is clear from context.

Writing denotes a set , and we write for when we require that .
The set Names is a finite countable set of names, and we use to denote elements of this set. We use to range over natural numbers, to range over rationals, and to range over booleans.
Similarly, we introduce the natural number drawn from the set of natural numbers , rational from rational number set and lastly boolean values from , denoted by . Variables are either a name or a name followed by a number of primes (). Type terms represent naturals, reals, Booleans, and products, respectively. A type environment is a partial function from variables to type terms. We treat environments as graphs of functions or as functions.
Function names are drawn from a fixed set containing basic operators. Expressions include constants, variables, vector expressions, vector indexing, function application, time derivatives, and partial derivatives. Derivatives can be applied to both expressions and variables. The time derivative on a variable, for example , has special status, in that it can both be used in expressions to mean the value of the derivative at a given time and can also be equated to a value. When there is a constraint that equates a time derivative of a variable to a value, the effect is that integration is used to compute the value of the variable itself. In principle, in an equational language, if a symbolic expression for the variable is known, the derivative variable can be determined from that expression. In practice, it is generally rare that a closed form expression for the result of a simulation is known. Instead, it is more common to have the value of the derivative known, and then numerical integration is used to compute the value of the variable itself. The partial derivative () is an operator that takes two expressions and returns the result of the first expression differentiated with respect to the second expression. The ASCIIbased syntax is . For two scalars, the result is simply the first expression partially differentiated with the second one. If one expression is a scalar, and the other a vector, the operator is applied componentwise. Allowing arbitrary expressions instead of just variables in partial derivatives, allows us to express things like the EulerLagrange equation directly.
The first type of equations is a continuous equation. In processing such equations, we distinguish between two cases, one where the left hand side is a variable, and the other when the left hand side is not a variable. This will be used to illustrate that the formalism is able to accommodate languages where equations need not be directed. The second type of equation is the discrete assignment. “ is reset to e”. Discrete assignments are essential for modeling hybrid systems, where instantaneous changes of a value (resets) can occur in juxtaposition to continuous dynamics. The third type of equations is a conditional equation, which allows us to express the choice between which of two equations holds depending on the boolean condition given as an expression. The fourth kind of equations is a universal quantification, and it provides a concise way of describing the dynamics of a system that has a family of state variables. The variable introduced by this construct may only be an unprimed name. The last construct is a set of equations .
3.1 Type System
A Acumen17expression has type under environment when the judgement is derivable according to the rules presented in Fig. 1. The rules for natural, real, and boolean constants are straightforward. The rule for unprimed names is simple environment lookup. The rule for primed variables, however, requires that both primed and unprimed variables have type real. The rule for vector construction is also straightforward. Vector indexing is a bit more interesting, as it treats the case when the index is a literal as a special case, allowing elements to have different types. This makes it possible to use vectors both for tuples and for (homogeneous) vectors. Function applications assume that we have a function that determines the arity of the function , and a function that determines the type of the th argument to the function. Partial derivatives have straightforward rules. The rules for equations are straightforward. Finally, environment extension of environment with the binding , written is an environment defined as follows:
3.2 Example: A Lagrangian Model




For a variety of technical reasons, researchers working on novel robotic systems tend to make extensive use of the Lagrangian method. It is especially useful in the case when the system being described has more than one state variable. Then modeling using Lagrange employs families of equations, which are written as one equation but really represent a collection of different equations derived by instantiating certain indices. Figure 2(a) provides the Acumen17 model of a second order nonlinear system shown in Figure 3. It consists of a pendulum hanging from a mass, which is attached through a spring to a wall.
As the system has two degrees of freedom, and , the example introduces a vector of state variables . The EulerLagrange equation that appears at the end of the example is expressed by the family of equations. In Figure 2(a), the quantifier is used to introduce the index variable for a family of equations. In the ASCIIbased syntax, the keyword
foreach
represent this quantifier. The intent is to express as concisely and as close to what would typically appear in a mechanics textbook:
This notation generally has a syntactic interpretation, that is, the contained in the th element of the vector is looked up. In other words, this family of equations literally represents the following two equations:
The offline partial evaluation strategy enables us to support family of equations and partial differentiation by utilizing the two most important components, namely the bindingtime analysis (BTA) and specialization. A successful BTA annotates the model with instructions for performing certain part of the computation early and other part for later processing. The annotated model for the pendulum/spring example is illustrated in Fig. 2(b). In this illustration, computations that remain for further processing are shaded grey, whereas computations can be performed immediately in the next specialization phase appear in a white background. The value of is marked known or static as it is a value, and the BTA also annotates variable known for that both and are known variables. A more interesting case is the indexing operator . Although the state variable vector being marked unknown, in fact we need to solve for the values of state variables and in the simulation, we can still perform this operation statically for the reason that the size of and the index variable are known.
The step which performs the work that a BTA schedules is called specialization. The result of specializing our running example is presented in Fig. 2(c). Computing the value of is simple rational arithmetic. The instantiation of a family of equations is essentially a type of iteration, which also replaces by and by vector lookup. In the same time, symbolic time and partial differentiation are performed using the chain rule. Solving multiple implicit ODEs to explicit form equations is achieved using an analog of symbolic Gaussian elimination. For our running example, the result of this step is presented in Fig. 2(d). Abstract interpretation with interval analysis is used to ensure the pivot expression is non zero. To control the system, for example, stabilizing the position and the angular displacement, one can add two PD controllers. The modification to the original model in Fig. 2(a) are as follows:
The Acumen implementation supports an enclosure simulation semantics that produces rigorous overapproximations (guaranteed upper and lower bounds) for all simulations [adam]. Previously, this implementation only supported a formalism that worked with hybrid ODEs. With the work we presented here, this implementation can now handle models such as the pendulum spring mass model presented above. The plot of controlled system are as follows:
3.3 A Cam and Follower Example
We further demonstrate the expressiveness of the proposed language using the following two case studies. Transforming rotational motion into any other motions is often conveniently accomplished by means of a cam mechanism. A cam is defined as a machine element having a curved outline, which by its rotation motion, gives a predetermined motion to another element, which is often called follower. Fig .4(a) shows such a cam mechanism, the curved outline of cam is a function of rotational angle , defined as below:
In the study of various aspects of the follower motion, the velocity and acceleration of the follower are needed. To get the correct form, the modeler usually has to manually derive the partial derivatives. Fig. LABEL:fig:armcode in the Appendix shows the mathematical model and the corresponding Acumen17program. Clearly, supporting partial derivatives in the language greatly simplifies the modeling task, and can save the modeler much tedious and errorprone work.
3.4 A Compass Gait Biped Example
The Compass gait biped model [cornellBiped, Aaron2D] is a two dimensional unactuated rigid body system placed on a downward surface inclined at a fixed angle from the horizontal plane. A diagram of the model is shown in Fig .4(c) with its physical parameters. The configuration of this twolink mechanism can be described by the generalized coordinates , where is the angle from the vertical line to the stance leg and is the angle between two legs. It is a hybrid model featuring two phases. At the start of each step, the system is governed by its continuous dynamics until the swing leg hits the ground. The discrete event can be modeled as an inelastic collision conserving angular momentum. The stance and swing legs switch instantaneously during the collision and go into the next step after.
Continuous Dynamics and Discrete Event
The continuous phase of this system can be modeled using the same Lagrange method shown earlier. Let point denote the position of centralized masses shown in Fig .4(b), form which its easy to define the kinetic and potential energy of the system. Applying the same Lagrange equation shown in Equation 1 with , we have the dynamic equations of the system during the swing phase. The perpendicular distance from the walking surface to the tip of the swing leg is given by
Where is the slope of the ground. Impact occurs when the tip of the swing leg hits the walking surface in a downward direction, which can be describe as follows: . Using conservation of angular momentum [goswamicompass], the explicit solution of post impact velocities can be determined. Fig. LABEL:fig:armcode in the Appendix shows the mathematical model and the full Acumen17model. This example shows the proposed formalism can support a direct mapping from mathematical model to simulation code for a hybrid system model with complex dynamics.
Binding Time  

Expression  
Equation  
Binding Time type  
Binding Time Environment 
4 BTA and Specialization for Acumen17
This section presents a declarative specification of the bindingtime analysis (BTA) and specialization process for Acumen17, and proves the correctness (typesafety) of the BTA with respect to the specialization process.
4.1 Bta
BTA is the analysis performed in an offline partial evaluation system to determine, given some early or “static” inputs to a program, which of the program’s computation can be done at an early stage [Jones85]. Fig. 5 gives a declarative specification of the BTA. There are two binding times S and D, representing “static” and “dynamic” computations, respectively. Static is for compiletime computations that are done before the simulation starts, and dynamic is for computations done during the simulation proper. Expressions, equations, types, and type environments are all annotated with binding times.
The changes to the derivation rules are largely straightforward. Essentially, binding times are propagated with types. In addition, when multiple subexpressions occur, their binding times are combined using the least upper bound operator which returns static only if all arguments are static, otherwise returns dynamic. However, for vector indexing, when the index expression is static and the subexpressions dynamic, we will still perform the look up operation. Finally, the rule for primed variable requires the unprimed variable with the same name to be dynamic. And in the rule for vector indexing with a literal, where the literal is annotated as static, the binding time of the expression is the same as the corresponding entry.
4.2 Specialization
Fig. LABEL:fig:semantics in Appendix presents the bigstep semantics for the specialization process. Values include constant with static annotation, dynamic expression and vector of values. Normal form equations are straightforward, with the absence of universal quantification equation. The first auxiliary function is used to compute static function application. The last two are for eliminating total and partial derivatives using the chain rule. Function returns free variables in an expression and function extracts the left hand side variable of a directed equation.
Relation essentially specializes all subexpressions to values then combines them according to their binding times. Function application with static binding time returns the evaluation result of the corresponding function. Vector indexing with static index performs look up operation statically, even if the vector itself has dynamic binding time. Total and partial derivatives get eliminated statically using different rules depending on the what their subexpressions specialized to. The rules for relation are similar. However, they all require that the equation to be specialized does not contain free variables that are defined in the equations following it. For directed equations, in addition to specializing the right hand side expression, the rule also substitutes the result into the rest of the equations. Both relations can also generate err terms, which will be propagated to top level for error reporting. For example, the static index may be specialized to a natural number that is bigger than the size of the vector. However, one type of error we do not catch is the case of partial derivative , when can not be specialized to a variable . It is analogous to the traditional division by zero error.
4.3 Type Safety
Definition 1.
The erasure relation for , and is defined as follows:
Lemma 1 (Erasure preserves typablity).
Proof.
By induction on the derivation of and , respectively. ∎
Lemma 2.
Substitution type preservation
.
Proof.
By induction on the derivation of and respectively. ∎
Lemma 3 (Type Preservation).