Summer 2002 Report about Hybrid Systems Modelling

Simon Lacoste-Julien
School of Computer Science, McGill University
Montréal, Québec, Canada

August 30, 2002

Abstract

The three main parts of my researches this summer were to

• get familiar with Hybrid Systems and implement some simple physical models with the 522 Simulator¹;

• extend the 522 Simulator with a DAE solver written in Fortran (DASSL and DASRT);

• add a graphical interface for the modelling environment using AToM³

Contents

1  Introduction and Motivation

Why studying the modelling and simulation of Hybrid Systems? Well, one reason was that this would involve mathematics, physics and computer science, the three fields I'm majoring in; so this would be interesting for me. But as important was that Hybrid Systems are useful to model a lot of different physical systems: in Chemistry for waste management systems; in Biology for feedback control systems; in Mechanical Engineering for hydromechanical systems controlling the direction of an aircraft (say); etc. The ability to model those complex physical systems and to test them (using simulation) can be greatly enhanced with a proper modelling and simulation environment. So this explains why the development of modelling tools like the 522 Simulator, as well as even more general meta-modelling tools like AToM³ [which will be explained in section 4], is useful. For an introduction on Hybrid Systems, I refer the reader to a talk I've given in the MSDL this summer: MSDL Talk: Intro to Hybrid Systems Modelling [2]. It will be assumed in this paper that the reader is familiar with Hybrid Systems.

The goal for me this summer was to get a feel for multidisciplinary research and to get in touch with some of the many issues about Hybrids Systems. The approach was originally on the theoretical side, studying which could be the best formalism to model Hybrid Systems and which could be the best formalism to simulate them efficiently. But gradually, we went back to a more pragmatic view and the researches were oriented in improving the 522 Simulator. This gave impressive results, and at the same time gave us more insights about the theoretical sides of Hybrid Systems [as will be explained in section 4.5].

I can separate my summer work mainly in three parts. The first part was to read the literature about Hybrid Systems, to study the 522 Simulator and also implement some physical models with it in order to find its strengths and weaknesses. This part is described in section 2. One of the conclusions of this part was that the modelling of complex physical systems would be a lot simplified by the support of Differential Algebraic Equations (DAEs) for the specifications of the system. So this motivated the second part of my summer work, which was the extension of the 522 Simulator with a standard DAE solver, DASSL² written by Linda R. Petzold and available on netlib. I didn't find any Python version of DASSL, so I had to do the extension manually. The fact that DASSL was written in Fortran made the task even more challenging (and interesting). During this part, I learnt a great deal of low-level stuff, but also useful stuff for a physicist who wants to use older libraries written in Fortran; or for a computer scientist who wants to glue applications together... Finally, the third step was to build a graphical user interface for the modelling environment of the 522 Simulator, using AToM³. With this tool, it is now easier to ponder about ways to improve the modelling of Hybrid Systems, as will be described in section 4.

2  The 522 Simulator

2.1  One way to model Hybrid Systems: FSA + ODE

The goal of this control system is to keep the speed of the car between v_max and v_min. You can see the FSA structure as well as the ODEs present in each state. The guards (like 'v > v_max') can be implemented easily with a monitoring function like 'v_max - v' which triggers a transition when it passes from positive to negative values.

2.2  Overview of the 522 Simulator

The way the interfacing is done between a user model and the 522 Simulator is by creating classes with specific content, which should inherit from standard classes provided with the simulator. The required content of those classes is described briefly in section 4.3.1. For screenshots and examples of the 522 Simulator, see my talk.

2.3  Improving the 522 Simulator

On the other hand, a feature that seemed to lack in the 522 Simulator was the possibility to model the continuous behaviour of the system using Differential Algebraic Equations (DAEs) instead of only ODEs. Indeed, most laws in Physics are already expressed in a DAE format (like conservation of energy for two point masses, for example: m₁ (d[r\vec]₁/dt)²/2 + m₂ (d[r\vec]₂/dt)²/2 = E = const.). It is often painful to transform those equations to explicit ODEs (and sometimes impossible), so it would be nice if the simulator could do the job at the place of the modeller. And this would permit the human modeller to focus more on physical concepts rather than on the algebraic solving of the equations. For those reasons, we have decided to include a DAE solver in the 522 Simulator, step which will be described in the section.

3  Interfacing DASSL with Python

The solver we wanted to include in the 522 Simulator was DASRT, a modified version of DASSL, which also supported the monitoring function formalism. In other words, we could provide DASRT with a DAE as well as a vector of zero-crossing functions and it would detect if any zero-crossing functions had a zero during a given interval; and find the zero if there was one. Since it took a while before I succeeded to find the source file ddasrt.f for the solver on the internet, all the first tests for this project were done with DASSL [but they have almost the same interface, so this was no big deal].

I now give an overview of what was needed in order to extend the 522 Simulator with DASSL/DASRT. First of all, we wanted all the user modelling to be done in Python, so this meant that the DAE was a python function. This python function could be passed to DASSL in order to be solved. DASSL was written in Fortran, so this meant that I needed to find a way to call a Fortran routine (DASSL) from Python and also that this Fortran routine had the possibility to call-back a Python function (the DAE). Since Python provides an interface with C only and not Fortran, a C wrapper needed to be build between Fortran and C. This whole project of extension with call-back mechanism has been separated in 5 steps, which are all documented in the next subsections:

1. Learning the basics of extension in Python with a call-back mechanism and test this with a simple Euler solver.

2. Learning how to mix C and Fortran in a program.

3. Link step 1 and 2 by implementing a Euler solver in Fortran and using it from Python.

4. Now that I have the basic principles to be able to extend Python with DASSL manually, try to do it automatically with f2py [with the Euler solver first].

5. Finally, use all this knowledge to bring DASSL and DASRT in Python.

3.1  Extending Python with C with call-back mechanism

3.1.1  Fundamentals of Python extensions

First of all, let's clear out the vocabulary: extending Python with a C module means to call and use a C program from Python, which is the main program. On the other hand, embedding Python in a C module means to use Python from a C program, which is now the main program. The first use (extension) is easier and is now easily automated using SWIG (Simplified Wrapper Interface Generator). The other use (embedding) is harder and there is yet no standard tool to automate the process. This is due mainly to the fact that there lacks a conventional way of doing embedding, because C is a lot less ``standardized" language than Python. Finally, extending Python with C with a call-back mechanism means that Python will use a C program which can ``call-back" a Python function. This is the method we needed in our project. It is thus some kind of mixing of extension of Python (Python calls C) with embedding (C calls Python). But one important difference with embedding is that since Python is still the main program in this case, there is no need to start a new Python interpreter as is needed in the pure embedding case.

3.1.2  The C wrapper interface

Wrapper functions

You need one wrapper function per function that can be called directly from Python. This wrapper function takes care of translating Python arguments into C pointers; then calling the appropriate C functions which execute what the function is supposed to do; then finally another translation of the C returned value to a Python return value.

Registration table

This is a C struct that will give the association between the names of the methods called from Python and their appropriate C wrapper function that handles the method definition.

Initialization function

This is the function which is called by Python the first time the C module is imported in Python. It takes care of all the proper initialization. This is the only function which is directly called from Python so it is the only one which should need external linkage. The other wrapper functions should thus be declared static for better encapsulation.

3.1.3  Call-back

Call-back pointer

This is a global variable of type Py_Object which will hold the address of the Python function passed by Python as an argument when calling some C function which needs a call-back. The C wrapper of this function will initialize the value of the call-back pointer when parsing the arguments given by Python.

Call-back wrapper function

This will be the function called by the C code when wanting to call the Python function. It will first translate the C arguments into Python arguments. Then it will call the Python function (using the call-back pointer) using the PyEval_CallObject C/API function which will return the result of the Python object as a Py_Object pointer. Finally, the wrapper function will translate back this result into C types and return this to the caller. So in brief, it does exactly the opposite of what the normal wrapper functions do...

3.1.4  Euler example

3.1.5  Compiling the extension module as a shared library

In this project, we use dynamic linking since it is a lot more flexible (I can compile different versions of the same extension module and test them one after the other without having to recompile the whole Python interpreter) and thus faster to develop. The basic principle of a shared library is that it postpones the symbol linkage step to when the library will be imported into a running program (the Python interpreter here), by storing the required information in the shared object file.

The building of a shared library is quite system dependant. The best way to do it portably for Python extensions is to use the Python tool Distutils and to make a setup.py script. But before explaining this tool, I will give an example of how to build Euler.c on Unix and Windows, just as a reference.

Unix   For a Python extension on Unix, you need to build a shared object file (.so). To do this, you would first do:
gcc -c -fPIC -I/usr/include/python2.2 Euler.c -o Euler.o
where the -fPIC option is to produce Position Independent Code, and the path appended to the include path option -I is the path to the Python include directory which contains the Python header files (Python.h, and others). You would then do the final linking (here, there is none since I use a single file) with:
gcc -shared Euler.o -o Euler.so
where the important option here is -shared to indicate that we produce shared object code on Unix. Of course, we could have combined those two steps here since we were only using one file. More information about shared libraries on Linux can be found here.

Windows   On Windows, we build a special kind of a Dynamic Link Library (DLL) which has extension .pyd for Python. I give here what has worked on my Windows XP machine with the MinGW GNU compiler, but the safest way would be to use Distutils... I used the MinGW compiler distribution since it also had the g77 Fortran compiler and was one of the few Windows compilers supported by Distutils. More information about DLLs can be found here.

The building steps are the following. First, we build the object file (almost the same as Unix):
gcc -c -IC:\Python22\include Euler.c -o Euler.o
Then, we build the export list: a list of functions which will be accessible from the importing program. In our case, we only need to export initEuler. We put the list in a .def file, with the following simple format:

LIBRARY Euler.pyd
EXPORTS
initEuler

Distutils   And now to do all this in a portable way, we write a setup.py script which uses the module disutils. This is the Python equivalent of a Unix Makefile. The minimal and simple setup.py file I have used in this test is given in appendix B.2 or can be found here . More information about how to build a setup.py file for an extension module can be found in section 3.3 of the online Python documentation Distributing Python Modules. To use the setup.py script, you simply type:
python setup.py build
On Windows, you will need to add -compiler=mingw32 at the end of the command in order to specify which compiler you want distutils to use. All the options and library paths will be automatically given by distutils. That's why this tool is portable. The shared library will be automatically produced in a folder build/lib<platorm-name>/.

3.1.6  Using the extension module in Python

3.2  Mixing C and Fortran

As an example, I had created a program in Fortran which called C; and another one in C which called Fortran. For both cases, I was using the MinGW compiler on Windows XP (but the same commands and files also worked on Linux). Let's consider the first case.

3.2.1  Fortran calling C

         PROGRAM fCallsC
            EXTERNAL HELLO
            CALL HELLO()
         END PROGRAM fCallsC

void hello_() {
    printf("Hello world! I'm in C!!");
}

3.2.2  C calling Fortran

extern void hello_();
int main() {
    hello_();
    return 1;
}

         SUBROUTINE HELLO()
            PRINT *, 'Hello world!  I''m in Fortran!'
         END SUBROUTINE HELLO

3.2.3  Differences between C and Fortran

Arguments passing

Arguments in Fortran are all passed by reference; arguments in C are all passed by value. This means that the C code should be modified accordingly to use the right signature to call a Fortran routine (or changes accordingly its signature in order to be called from Fortran). All non-pointers types should send their address (using the ampersand &) instead of their value. This will be shown clearly in an example I've built in section 3.3.

Multidimensional arrays

Arrays indexing start at 0 in C, but at 1 (or bigger if the user has required it) in Fortran. Also, multidimensional arrays in C are stored in row-major order (row by row); whereas they are stored in column-major order in Fortran. This means that the value A[i][j] in C is the same as A(j,i) in Fortran for an array A in memory. So basically, to solve this issue, you simply inverse the indices in Fortran or in C.

3.3  Using Euler function in Fortran from Python

3.3.1  Building the Fortran extension module

In order to bypass this limitation, we compile the required Fortran file as a library that will be included by distutils. The code for this setup_fEuler.py file is given in appendix C.3 or is available here . The main difference with the setup file for the Euler solver in C is that this one uses library inclusion. So to build the Fortran extension in a portable way here, we do the following:

1. We compile the Fortran solver:
   g77 -c Euler.f

2. We create a library for this object file:
   ar -cur libEuler.a Euler.o

3. We finally build the extension with distutils:
   python setup_fEuler.py build [--compiler=mingw32]

3.4  Automatic building of C/Fortran extension for Python

Fortunately, there exists automatic C wrapper interface generation tools. To extend Python with a C module automatically, you can use SWIG, as mentioned in section 3.1.1. This tool also supports the handling of call-back functions as is explained in this section of the online documentation. But the interface to handle call-back mechanism is quite low-level (and we even need to handle some memory manually). On the other hand, the interface used by f2py, a Fortran to Python interface generator written by Pearu Peterson, is a lot cleaner. And this tool can even be used (with some tricks) to extend Python with C modules [see this email in the f2py discussion list for more info].

3.4.1  Introduction to f2py

F2py works by using a signature file [.pyf] written in a format similar to Fortran 90 and which describes the signatures of the functions which should be available from Python. F2py can generate automatically this signature file by reading the Fortran source code, and the user can modify it for the fine-tuning of the specifications. Unfortunately, the user guide currently [August 30, 2002] available on the f2py web site is quite outdated⁴. But Pearu Peterson is currently working on a new version of the documentation, which is available here. We'll describe in the next section how to use f2py to extend Python with our Euler.f example. But before, here are several notes to take in consideration:

• I was using version 2.13.175-1250 of f2py.

• On Windows, I was using the mingw gcc compilers, which works nicely with f2py. Depending of which version of f2py you use, it could happen that you need a patch to make it works with mingw. See the end of this email for the patch.

3.4.2  Using f2py to extend Python with Euler in Fortran

F2py is quite intelligent to guess the signature of the functions, and even of the call-back functions if those were called explicitly in the source file. But in the case of arguments which will be changed in-place (like in most of old Fortran 77 code), the user has to manually indicate the intent(inout) of those variables in the signature file. The main attributes used by f2py in the signature file are well-explained in the old user guide. The intent(<intentspec>) attribute specifies how a variable is going to be used in a routine. The main ones needed for basic uses are: intent(in) which means that the variable is an input variable and is assumed not to be changed by the function; intent(out) which means that the variable is only an output variable (f2py will then generate a different Python signature than the Fortran one for this function, see the old user guide for more details); and intent(inout) which means that the variable is used as input AND output, and can be changed in place by the Fortran routine (or the Python one when in the specification of a call-back function). When there are only a few modifications needed for the signature file, those can be directly included in the Fortran source file using the cf2py tag at the beginning of each line. In the case of our Fortran Euler solver (which was already described in section 3.3), we only needed to add the specification that two variables were used as arguments which could be changed in-place; all the rest of the signature was correctly detected by f2py. So the following two lines were simply added to the file Euler.f (to yield the f2pyEuler.f file which is available here ):
cf2py intent(inout) res
cf2py intent(inout) newx
With this modification, we could build directly the extension module from the source file, without using the intermediate .pyf file:
f2py -c f2pyEuler.f -m f2pyEuler
We still needed to specify the module name for the extension module (using the -m option) since this wasn't specified in the source file. This command will produce a f2pyEuler.so shared library which can then be used as before. A test script for this extension is available here . Note that the numpy package was needed since f2py uses only Numeric arrays.

3.5  The extension of Python with DASSL/DASRT

3.5.1  Building DASSL with f2py

3.5.2  How DASSL works in Fortran

RES

This is a subroutine provided by the user which represents the DAE system to be solved. More info about its format will be given in section 3.6, when I will explain how to use DASRT from Python.

neq

An integer representing the number of equations of the system.

t

The current value of the independent variable (time). This variable is modified in-place by the solver as it solves the equation (and time advances). This means that it should be declared as intent(inout) in f2py.

Y, YDOT

These arrays contain the solution components at t and its derivatives (respectively). They are also modified by the solver as time advances (this is how the user can find back the solution at t).

tout

The time at which the solution is desired.

INFO

This is an integer array of size 15 which gives information on how the solving task must be carried out. For a default use, all components should be equal to 0.

rtol, atol

Those quantities represent the relative and absolute error tolerance which are provided by the user. They can be modified by the solver.

idid

This scalar communicates to the user what the solver did.

RWORK, IWORK

These are real and integer work arrays for the solver, respectively of length lrw and liw.

RPAR, IPAR

These are real and integer parameter arrays respectively which can be used by the user to communicate between the calling program and the RES and JAC subroutines.

JAC

This is a subroutine which can be provided by the user to compute the Jacobian matrix (of partial derivatives) of the DAE. If INFO(5) is set to 0, then the Jacobian is computed automatically by DASSL; if it is set to 1, then DASSL uses the JAC subroutine provided instead.

3.5.3  Notes about f2py

Secondly, because f2py is still in development process, it could happen that there are bugs with your generated extension module. To try to debug your extension module, you can use the -debug-capi function while building it with f2py: this will make the C wrapper interface to be very verbose about what it is doing. You can then post the bug to Pearu Peterson on the f2py discussion list; he usually answers email very rapidly. In my case, once I had figured how to do the extension properly, everything has been working fine. But before I got there, I passed through a lot of debugging and segmentation faults (the omnipresent bug when working with extension modules)... I have explained how to debug shared libraries using gdb in appendix A, because I think it is very useful to understand the working of extension module. This was the only way in which I have succeeded to find the source of my segmentation faults. These came from a subtle point, so it's worth mentioning here. Fortran supports happily unbounded arrays (i.e. array without any bound declared), and those were used heavily in ddassl.f. On the other hand, to translate Fortran arrays into Python arrays, you need to specify the bounds. So that's why that the bounds of any array used in a call-back function need to be declared. This wasn't checked by the version of f2py I was using so it produced happily an extension module which could crash. On the other hand, you can use unbounded arrays in the signature of a routine called by Python since Python keeps the size information with its arrays. So that's why I could use unbounded arrays (using dimension(*)) in the signature of the ddassl routine given in ddassl.pyf (see appendix D), but all the arrays used in the call-back routines signature (RES and JAC) had their bounds declared. In my test Euler solver in Fortran (see section 3.3), I had passed the size of the state vector (which corresponded to the number of equations in the ODE) as an argument of the call-back routine. There was no such argument in the signature of the call-back functions used in ddassl.f, so a trick was needed in order to know what was the size of the arrays. The IPAR argument in the call-back functions was an integer parameters array which could be used by the user to communicate between the caller of DASSL and the call-back routines. I used as a convention that the three first elements of IPAR represented the dimensions of the needed arrays (see my modified ddassl.pyf file).

In the case that you would want to study the C wrapper file generated by the f2py, you can produce it by executing f2py without the -c option on the .pyf file:
f2py ddassl.pyf
This will produces the C wrapper file dasslmodule.c (since the module name given in in ddassl.pyf was dassl).

If you want to reinstall a new version of f2py (this could happen if you're unlucky enough to download a buggy version of f2py), you can simply delete the content of <path to your python>/lib/site-packages/f2py2e, and then reinstall the new version normally.

Finally, note that the interface I have built with DASSL is very primitive. F2py provides very powerful features to modify the functions signature in a clean way (for example, it could automatically provide the work arrays; it could copy the input arrays that need to be modified in-place and return the copy [so without disturbing the original arrays, as shall be expected in a Python way of passing arguments]; etc.). With these changes, the way to call DASSL from Python could be simplified a lot. But for two reasons, I have decided to use the simplest way of interfacing DASSL with Python (simply recopying the Fortran signature in Python). First, most of the f2py features which could accomplish that were no more documented, so it was tricky to use them. Secondly, I was not sure how I was going to use DASSL from Python yet, and it would be easy to build a wrapper in Python around the imported ddassl extension module in order to give a cleaner interface. Since I was building only a prototype, the primitive interfacing was enough. An example of question that could arise when trying to develop an efficient implementation of the solver in Python could be: is copying the input array automatically by f2py (as mentioned above) really a good thing? For example, if we run DASSL in intermediate mode, the control could pass very often between Python and ddassl.f, and each time f2py would do a copying of the input array. On the other hand, doing the copying only once in Python when starting the solver (and not each time the solver is called) would avoid all this overhead. This kind of question explains why I have postponed the clean interfacing of the solver after the development of the prototype.

I have also to mention that using intent(inout) arguments in my interface is against the philosophy of Pearu Peterson, who suggested to use always intent(in, out) instead. This issue is discussed in this email. From what I have understood, using intent(in,out) has two main effects: it will append the variable as a return value in a tuple for the function; it will make the changes to the variable more robust than if intent(inout) had been used (for example, it seems that when intent(inout) is used, it could happen that changes to the variable made by the Fortran code could not be transmitted to Python under certain conditions). But because I wanted to mimic the Fortran interface in Python for simplicity, I decided to still use intent(inout) for the arguments of my interface function in DASSL.

3.5.4  Building DASRT

3.6  Inclusion of those solvers in the 522 Simulator

The code for the 522 Simulator with the new solvers can be found here . This comes with all the solvers already built on Windows and Linux (.pyd and .so). In case they don't work on your machine, you can rebuild all of them by following the instructions given in the readme file in each Solvers_build/ subfolder, and then move the rebuilt solvers in the root directory (Simulator522/). Each subfolder also contains a little Python test script to test the solver (the same that were given in this document). In order to run the 522 Simulator with the new solvers, you will need numpy. And you will of course need f2py if you want to rebuild some of the solvers. To start the simulator, you simply execute the module GUI_new.py.

3.7  Conclusion

4  Building a GUI modelling environment for the 522 Simulator

• A GUI is usually more user friendly than raw code and it permits the user to focus on the model rather than how to implement it for this particular Simulator.

• To build a good GUI, you need precise specifications; so this pushed us to really think about what we wanted to include in our Hybrid Systems formalism. Those specifications (in the meta-model which will be explained later) also permit model checking at the syntax level (using constraints in AToM³).

• By adding another layer of interface between the simulation environment and the modelling level, we can obtain a yet cleaner interface. For example, there was some repetition required in the Python code for the models in the 522 Simulator (as for the variable initialization for example), and you often needed to add self. at many places. In the GUI I have made, you don't need all this redundancy (which was error prone and time consuming): you only need the (almost) ``minimal" information required to specify the model.

The building of a GUI for the 522 Simulator was done in two steps: first of all, the specifications for the modelling environment for Hybrid Systems were given to AToM³ which could then generate automatically a modelling environment tool (see section 4.1). Secondly, I wrote a code generation tool (using graph grammars in AToM³) which could transform a Hybrid Systems model instance built in my GUI to a model instance understandable by the 522 Simulator (i.e. in Python code with the original interface of the 522 Simulator). Before explaining those two steps, I'll give a brief introduction to Meta-Modelling with AToM³.

4.1  Introduction to Meta-Modelling and AToM³

This being said, we can now present AToM³. AToM³ stands for ``A Tool for Multi-formalism and Meta-Modelling''. It was developed originally by Juan de Lara Jaramillo and Hans Vangheluwe at the MSDL of McGill University. AToM³ is a visual tool to write meta-models and generate automatically the (code for the) modelling environment corresponding to the meta-model; and also to write and execute formalism transformations (using graph grammars, since (almost) everything in AToM³ is represented internally by graphs). A good tutorial about meta-modelling with AToM³ was written by my colleague Marc Provost. The following paper [4] also gives a good introduction on how to create modelling environment with AToM³. In this case, the example presented is the building of a modelling environment for simple Causal Block Diagrams (CBD) and the generation of code in the Object-Oriented Continuous Simulation Language (OOCSMP)⁶. My approach in this project is thus very similar: I write a meta-model for Hybrid Systems and a code generator for the 522 Simulator language.

4.2  Meta-model of Hybrid Systems in AToM³

4.2.1  General description of the meta-model

I will first give a high-level description of the meta-model. The meta-model was directly inspired from the informal description of Hybrid Systems that I had given in section 2.1. It consists of a mix of FSA, ODE and DEVS. The two entities in the meta-model are used to model the hybrid behaviour of a HS: the Mode entity represents a possible continuous mode for the system (e.g. free_fall, accelerating, braking, etc.) (and so contains the ODE); the Event_HS⁷ entity represents a possible discrete event which can happen (like start_braking, start_accelerating, collision, etc.) (and so represent the events in the DEVS). A Mode node will contain a ODE describing the state of the system in this mode. The monitoring functions used in this mode are represented by State_Guard relationships connecting Mode to Event_HS, the latter identifying which event is triggered when a zero has been detected in the given monitoring function. The Event_HS node will contain action code which will represent what is happening when the event is triggered. The event can also schedule new events using a Schedule relationship; and can cause a mode transition using the Mode_Transition relationship (this represents the transition for the FSA). The variables and parameters used in the model are entered as global model attributes. The file used to save the meta-model was HybridSystems_ER_mdl.py⁸. It can be found in the HybridSystems folder which contains all the files needed for the new 522 Simulator (a link to it will be given in section 4.4). Note that for the remaining of this part, all the paths will refer to subdirectories of the HybridSystems folder.

Note that with the modelling environment specified by this meta-model, you can express a lot of different Hybrid Systems models. You can use the event scheduling component to implement if-constructs and while-loops, so this means that you can express almost any computable function. I think that anything that could be added to the modelling environment could be considered simply as `syntactic sugar' to facilitate the modelling (like implementing a real while-loop construct).

4.2.2  Generated modelling environment

4.2.3  Bouncing ball example

The Check_Collision State_Guard is connecting the Free_Ball Mode to the Collision Event_HS, with the meaning that when a zero is detected in the monitoring function given in the given State_Guard attributes (here, y as is shown graphically) in the correct direction (here, from + to -), then the given Event_HS is triggered, i.e. its action code, mode transition and schedule events are executed. In this case, Collision contains only the line v = - k * v as action code, which reverses the speed of the ball with a coefficient of restitution k (will be smaller than 1 for inelastic collision; equal to 1 for elastic collision; and greater to 1 for superelastic collision). The Mode Transition to the Free_Ball mode indicates that after this (instantaneous) event, the system is put back to the Free_Ball continuous mode. But there is also a Schedule relationship which connects the Collision event to the getStuck event. The condition for the scheduling to occur is shown graphically with IF t > t_stuck, i.e. the event getStuck will be scheduled in 0 second (because of the AFTER 0) if t (time) is greater than the parameter t_stuck. The event getStuck simply sets the speed to 0 and then puts the system in the Stuck mode.

So an informal denotational semantic for this model in the (informal) ODE + control flow formalism could be the following:

initial conditions:
y = y_0
v = v_0
initial ODE:
dy/dt = v
dv/dt = -g
when y = 0:
    if t > t_stuck then:
        v = 0
        use empty ODE for the remaining of the run (dy/dt = 0; dv/dt = 0)
    else:
        v = -k*v

4.2.4  Detailed description of the meta-model

Global attributes:  

Model_Name

A string for the name of the model.

Author

The author(!).

Description

A short description for the model. It will be included in the generated code.

Documentation

More documentation for the model (as what each variable represents; where special care should be taken with parameters, etc.). It will also be included in the generated code.

Variables_List

A list of strings giving the names of the variables used in the model.

Parameters_List

This is a list of tuples (see figure 6¹⁰ ) of strings, to define the global parameters used in the model. The first element (string) of the tuple represents the parameter name; the second string gives what is the initial value for the parameter.

Import_Modules

This is a list of strings representing Python modules that the user want to use in commands inside the model. For example, the random module can be used to produce random numbers (see the Train example, section 4.2.5).

Mode attributes:  

ODEs

A list of tuples of strings representing the ODE for the continuous mode: (<variable on LHS of ODE>, <allowed python expression for derivative>). The first element of the tuple is the variable name that is differentiated with respect to time (LHS of the equation); the second element is a (Python) expression giving the RHS of the equation. For example, dv/dt = -g is encoded with a tuple ("v","- g"). Notice the space between the '-' and the 'g'; this is because ANY model variable name or parameter name used in the expression needs to be separated by a space for the correct parsing of the expression when generating code.

Only variable names already defined in the global attribute Variables_List are allowed as the first element of the tuple. What is allowed in the expression string in the second element of the tuple is any Python expression, with the restriction that variable or parameter names be separated by spaces. You can use external module functions <module_name>.<method> provided that <module_name> was defined in the Import_Modules global attribute. Also, to define non-autonomous ODE, you can simply use the special variable 't' which represents time.

Local_Parameters

A list of tuples of strings representing declaration of local parameters to the mode: (<parameter_name>, <init_expression>). If the <parameter_name> is the same as one already declared at the global level (in Parameters_List), then it will override the latter. The <init_expression> is the initialization expression of the local parameter, and have the same syntax and restriction as the derivative expression in the attribute ODEs. But the meaning of a state variable in this expression is the value of the variable when entering the mode. For example, having (k,v) with v a state variable would mean that the local parameter k is set to the value that v has when the system is entering the mode. Similarly, 't' will mean the time when the system is entering the mode. Finally, you may (for convenience) use other local parameters in the expression of the RHS with the condition that you are careful with the order in the list: you can only reference local parameters which were defined before in the list. This can be useful if you want to build complex expression with lots of constants in steps.
NOTE: there is no way the value of the local parameters can be changed from outside the mode (for example, by an event). If you want to change parameters with events, you have to define them as global parameters by putting them in the Parameters_List global attribute.

Name

The name of the mode (must be unique).

Documentation

Documentation for the mode (will be added at the corresponding place in the generated code).

Event_HS attributes:  

Name

The name of event.

Documentation

Its documentation (will be put as a comment in the generated code in the definition of the respective event handler).

is_Start_Event

This is a boolean to indicate if the given event is the Start Event, i.e. the first event that will be executed when simulating the model (you should have exactly one Start Event in the model; this is usually where you also do all the state variable initialization). The value of this boolean influences the appearance of the event in the model.

Action_Code

This is a list of strings which represents a list of actions to execute when the event is triggered. Each string represent a line of a Python statement. You can use again model variable names and global parameters if they are separated by spaces. You can use 't' for time. You could even do loop and stuff by formatting it properly (e.g.: build the list ("for x in [1,2,3]:", "    print x") ).
TRICK: In order to make reference to a local parameter in the current mode (i.e. the mode of the system WHEN the event is triggered), you can use the special construct self.diffeqs.<param_name>.

State_Guard attributes:  

Name

The name of monitoring function.

Documentation

Its documentation (will be put as a comment in the generated code in the definition of the monitoring function).

Monitoring_fct

This is a string which represent a Python expression which gives the return value of the monitoring function. The syntax and restriction for the expression is the same as the one for the derivative expression in the ODEs attribute in the Mode entity.

Direction

This is an enum in which you choose which direction of crossing that will be monitored by this function ('+-' is from positive to negative; '-+' is for the opposite; 'any' is for any direction).

Schedule attributes  

Condition

This is a string representing a boolean expression in Python (0 is false, non-zero is true). It can use model variables, parameters, time 't', external module functions, etc. as in the ODEs attribute of the Mode entity (except it can't make direct reference to local parameters).

Schedule_in

This is a string representing an arithmetic Python expression which will tell in how long the connected event will be scheduled in. For the syntax, it has the same possibilities as for Condition.

Mode_Transition attributes   It has none.

4.2.5  Other examples

Train Example   This was the main example that Hans had used to motivate the requirements for the 522 Simulator, since it contained everything that was wanted for the simulator: discrete behaviour, continuous behaviour and event scheduling. The model built with the HS formalism is shown in figure 7. The file containing the model is models_HS/TrainSchedule_HS_mdl.pyc.

This is a model of a train which travels in one dimension, and which takes passengers at a station and then leaves them at another station. The starting situation is the train at rest at the first station, waiting for 10 passengers to arrive. The passengers arrive at a random time interval between 0 and 10 seconds (using the function random.uniform). When the 10 passengers have arrived, the train departure is scheduled in 5 seconds. The train then starts accelerating. Its speed is controlled by bang-bang control between v_max and v_min (when it is not accelerating, air resistance and friction makes it decelerate (this is encoded in the Mode FrictionODE)). When the train reaches the position stopping_x, it starts braking in order to be at rest when reaching the second station. When it becomes at rest, the passenger starts to leave at constant interval of 5 seconds. All of this can be easily understood by following carefully the model in figure 7.

Heating Pot   This is an example that Hans had used extensively in talks he gave about the use of multiple formalisms. The detailed description of the heating pot example can be found in p.14 of his notes about the Modelling and Simulation Process for the course 308-522A. Briefly, it represents a pot of water which can be warmed at the rate of W Joules per second, and filled with water at temperature T_in at the rate of phi cm/s (represents the level variation). The temperature of the pot needs to be kept between two thresholds, as well as its level. This is why you have four monitoring functions. This model is a good example of how you can use local variables, and how you can use parameters to change the state of a system (since at first sight, we could wonder where is the discrete change in this model because there is only one Mode). The main differential equation is

dT dt =  1.0 l (heat - phi ·filling ·(T - Tin))

where

heat =  W c ·rho ·A

where c is the specific heat of the liquid, rho is its density, A is the cross-section surface of the pot, l is the level of the liquid, filling is a parameter set to 1 when filling, 0 when emptying the pot. Also, you can notice the line 'self.W_init = W' in the action code of the Initialize_Model event. This defines a variable W_init which is used to remember what was the initial heating rate. Then, this rate can be set to 0 when the heating is stopped (in event stopWarming) or reset to its initial value when the heating is restarted (in event startWarming). This example also demonstrates the poor display feature that AToM³ has for now: the equations are not displayed completely. Another way to build this model will be explained in section 4.5.

4.3  Code generation using graph grammars

There are two main tasks that I needed to accomplish to build a code generation tool for HS in AToM³. The first one was to understand well the structure of the model in the HS formalism, compare it with the one needed for the 522 Simulator and then create a corresponding code generation algorithm. This will be explained in section 4.3.1. The second task was to implement this algorithm, and to do this, I needed to find how to navigate through the AToM³ data structures and how to get information about the model. This will be explained briefly in section 4.3.2.

4.3.1  Code generation algorithm

• One class per continous mode (they are subclasses of Model.FirstOrderDEs), which gives the specifications of the ODE in each mode.

• One model class (subclass of Model.Model) which does the following:
  • defines the list of state variables and parameters using the attributes modelvarnames and params respectively;

  • defines the start function which is used to initialize the model: it initializes the state variables using certain parameters as initial conditions; it initializes the diffeqs variable to the starting ODE class (possibly an empty ODE); it initializes the zerocrossingfuncs variable to the current list of zero-crossing functions in use;

  • defines all the zero-crossing functions needed in the model;

  • defines all the event handlers functions needed in the model.

1. Prompt user for the file name of the generated file

2. Write header, doc and import statements [using information from the global attributes in the model such as the Import_Modules list]

3. For each Mode in the model, write an ODE class:
   1. Write __init__:
      • Write its signature: probe ODEs and Local_Parameters to see what variables and global parameters were used

      • Build varnames and varderivs

      • Initialize parameters (local AND global used) [use the expressions in Local_Parameters for the local parameter; in Parameters_List for the global ones]

   2. For each ODE variable, write the derivative functions [use RHS of ODEs]

4. Write the model class:
   1. Write __init__:
      • Signature: global parameters with default value [from Parameters_List]

      • Initialize global parameters; and set variables to 0.0

   2. Write the start() function:
      • Find the Start_Event and just schedule it from this function

   3. For each Event, write an event handler function definition:
      1. Write the Action_Code

      2. If there is a Mode_Transition, then
         • Assign the new Mode ODE instance to diffeqs

         • For each State_Guard connected to this new Mode, add it to zerocrossingfuncs

      3. For each Schedule connection; write an if test using Condition and put in its body the scheduling of the event which is connected to it, in Schedule_in seconds

   4. For each State_Guard, write the zero-crossing function definition [using Monitoring_fct expression]

Some validation of the model is also done while generating code. Conceptually, this should really be done at the modelling level (and not when generating code), but for a prototype, the validation was done wherever it was possible to do some (i.e. I haven't made sure that everything in the model was validated, but when I had the possibility to validate something, I used this opportunity). The only validation which is done at the modelling level automatically by AToM³ is the cardinality checking for the connections between objects; I haven't put any constraint in the meta-model of the HS formalism. Here are the few checks that are done while generating code:

• Check that the variables on the LHS of the ODE declarations (in ODEs) are in the Model_variables list;

• Check that ODEs and Local_Parameters are list of tuples;

• Verify that there is exactly one Start_Event;

• Delete unused events and inform the user about this (unused events are events which don't have any incoming connections);

4.3.2  Navigating through AToM³ data structures

• First of all, we need to get an hold on an AToM³ instance from the graph grammar; you can access it using: self.rewritingSystem.parent when you are defining a GG rule. Say we call it atom3 for the remainder of this section.

• You can now access the root of the abstract syntax graph for the model (ASG), with atom3.ASGroot. All the global attributes of the model are stored in it. You can obtain the value of an AToM³ attribute with getValue(). For example, to get the author name of the model, you can do atom3.ASGroot.Author.getValue().

• We can obtain the list of nodes in the model graph by accessing atom3.ASGroot.listNodes, which is a dictionary which contains lists of pointers to each node (indexed by node type). For example, you can access all the Mode nodes in the model with atom3.ASGroot.listNodes['Mode']. This is used to traverse the whole model. The attributes for a specify node is simply stored in this node. For example, can get the name of the first Mode with atom3.ASGroot.listNodes['Mode'][0].Name.getValue().

• To find the connections of one node, you use <node>.in_connections_ or out_connections_ which are lists of pointers to incoming connected nodes and outgoing connected nodes, respectively. Note that a relationship is also considered as a node, so you usually have to pass through two connections before getting to a new node of the same type. To check for the type of a node, you use the getClass() function. For example, here is a little code excerpt which would find the Mode connected to an event through a Mode_Transition:

  # assume event contains an Event_HS node
  for node in event.out_connections_:
      if node.getClass() == 'Mode_Transition':
          # Note that this node is a RELATIONSHIP;
          # Also, by the cardinality check of AToM3,
          # we know that only one Mode can be connected
          # (outgoing) from a Mode_Transition,
          # so we can simply take the first one in the list.
          mode = node.out_connections_[0]
          return mode
  
  

• A final important thing to know: the return values of boolean and enum types are not necessarily standard. <ATOM3Boolean>.getValue() will return (<display text>, <value>). <ATOM3Enum>.getValue() will return a tuple with the list of possible values as the first element, and the index of the current value of the enum as the second element.

4.4  How to use and install the new 522 Simulator

To simulate your HS models, you first have to generate code using the 'Code Generation' button. Then, you can load the generated file in the simulator as a model file.

Appendix E contains a list of all the important files and directory found in the HybridSystems directory.

4.5  Some observations about Hybrid Systems

Statecharts would solve this problem since it has support for concurrent states, thus avoiding the combinatorial explosion when there are independent transitions characterizing a system. The Heating Pot model using statecharts is shown in figure 9. The question that one could ask then is what would be the semantic of two concurrent modes. Well, it could simply be that the RHS of the ODEs for the same variable name would be added. For example, if the system was in the no_Warming state and in the Filling state in figure 9, then the ODE for T would be:

dT dt =   F(T) from no_Warming state  +   H(T) from Filling state

This would be meaningful if the contributions could be simply added, but this is often the case in Physics where we have the superposition principle. Apart avoiding the combinatorial explosion, this way of modelling has the merit of clearly showing the concurrent nature of two different characteristic of the system (here, the level change vs. the heating).

Another thing we have found this summer is how to use DASSL (or DASSRT) to do acausal modelling, that is, to model physical systems without assigning any causality (like is done in the modelling language Modelica). For example, we would like to model a (electric) capacitor with the equation C ·dv/dt - i = 0 without having to say if it is i (the current) or v (the voltage) which is known. Depending of what happens around the capacitor, we would like sometimes to drive with i, other times (maybe in the same simulation run) with v. This can be done in DASSL by giving as the residual equation:
RES[1] = C*YDOT[1]-Y[2]
(assuming Y[1] = v and Y[2] = i) with the following constraint which will drive the desired variable:
RES[2] = Y[f.index] - f(t)
Now, f could be a Python object which would represent which variable is driven (f.i would be set to the index of this variable in Y), and f(t) would call the appropriate driving functions. The beauty of this is that f could be changed at run-time. This would truly permits acausal modelling, and would be useful when we will want to extend the Hybrid Systems formalism with input/output support.

5  Conclusion

5.1  This summer

Regarding the modelling and simulation of Hybrid Systems, it has been shown how AToM³ could be used to build a powerful modelling environment, quickly. The 522 Simulator has been improved with a GUI, as well as with the addition of the DASRT solver. The task of extending Python with Fortran code has been very instructive, and quite useful, as the several questions I received about it can show. Several improvements for the modelling and simulation of Hybrid Systems were suggested. I will discuss them in next section.

5.2  What's next

Other improvements that could be brought are the following:

• The addition of an algebraic solver so that equations (to be solved) could be put in the action code of an event (this would be quite useful for the initialization of the state variables after a mode transition; e.g., we could demand that the variables respect conservation of energy and momentum, and the simulator would find by itself the correct values under these constraints).

• Nice documentation features so that the user could check the meaning of each variable in the simulation environment.

• Adding units/dimensions support: this would permit dimensional analysis (and verification) at the modelling level.

• Modelling with Statecharts instead of FSA (as was mentioned in section 4.5).

• Adding the support for constrained variables and the connection of components (using input/output formalism). This would permit to model driven systems in a modular way; and also to do acausal modelling.

6  Acknowledgments

A  How to Debug a Shared Library with GDB

The problem that you may face when you want to debug a shared library with GDB is that the symbol table of the shared library is not loaded until it is imported. This means that you can't set breakpoints before it is imported. But if it crashes after it gets executed, how are you supposed to debug it if you can't stop it before it crashes? The solution resides in the magical SIGINT interrupt: CTRL-C which permits to break the execution of the program at any time. So the idea is to have your program stop just after it has loaded the shared library, but before it executes them. For example, in my Python extension modules, the execution of the test modules was always done using a run method. As a result, when I imported the module, the shared library were not executed: they were only loaded. This permitted me to do the following steps to debug the shared library:

• Compile your shared library with the -g option.

• In case you won't have an obvious place where you will have time to press the CTRL-C command after the loading of the library but before its execution, put a sleep statement in the source code where you want to break it.

• Run python under gdb (with a command like gdb -args python2.2 <test_module> for Unix).

• Do run to start the execution of the import statement; if you were importing shared libraries in the test_module, now they should be loaded.

• When you reach the sleep point, send a SIGINT interrupt (CTRL-C) to switch control back to gdb.

• Look with info shared to see if the symbols have been read (look at the 'syms read' column); if they haven't been read, then you can reload them with the command share <file>.so.

• Now you can add the breakpoints as usual into your shared library, and start debugging it.

B  Source files related to the Euler solver written in C

B.1  Euler.c: source code for a Euler solver written in C

/* Written by Simon Lacoste-Julien on 2002/07/23 */
/* This is a test for Python Extension.
I'm implementing the Euler solver used in the 522 Simulator.
Interface of use:
>>> import Euler
>>> newx = Euler.step(f,x,t,delta_t) # as before!!!
Where f is an array of (derivative functions) [of the form g(x,t)], x is initial array of values,
t is time, and delta_t is step length (both are double).  That's it!
Magic of Python API... */
/* newx = Euler.step(f,x,t,dt) will simply yield:
newx = x + f(x,t) * dt
*/

#include <Python.h>

/***** 2) Module method definition ********/
static PyObject*
step(PyObject* self, PyObject* args) {

    PyObject* f_array = NULL;       /* will hold f */
    PyObject* f_call = NULL;        /* will hold f[i] (to be called) */
    PyObject* x_array= NULL;        /* will hold array x */
    PyObject* newx = NULL;          /* will hold result (array)*/
    PyObject* xObject = NULL;       /* will hold Python x[i] */
    PyObject* arglist = NULL;       /* will hold arguments for call of f[i] */
    int length;             /* size of f and x */
    int i;                  /* array index */
    double x, x2, t, delta_t;   /* x will hold x[i], x2 will hold result of computations  */

    if (PyArg_ParseTuple(args, "OOdd",
            &f_array, &x_array, &t, &delta_t) == 0) {       /* converts Python -> C */
        /* setting error message: */
        PyErr_SetString(PyExc_RuntimeError, "Wasn't able to parse args for Euler.step fct.");
        return NULL;        /* when returning null, you raise an exception */
    } else {
        Py_INCREF(f_array);
        Py_INCREF(x_array);
        if((PyList_Check(f_array) && PyList_Check(x_array)) == 0) {
            PyErr_SetString(PyExc_TypeError, "Error: f or x weren't list!");
            goto error;
        }

        length = PySequence_Length(f_array);    /* assume x_array has same length... */
        newx = PyList_New(length);              /* initializing the return array */
        Py_INCREF(newx);

        /****** Euler loop ********/
          for(i = 0; i < length; i++) {
            Py_XDECREF(f_call);
            Py_XDECREF(xObject);
            f_call = PyList_GetItem(f_array, i);
            xObject = PyList_GetItem(x_array, i);
            Py_INCREF(f_call);
            Py_INCREF(xObject);

            x = PyFloat_AsDouble(xObject);      /* set x = (double) xObject */

            /* evaluating f[i](x,t) : */
            Py_XDECREF(arglist);
            arglist = Py_BuildValue("(O,d)", x_array, t);   /* building (x,t)  NOTE: x is a vector! */
            Py_INCREF(arglist);
            x2 = PyFloat_AsDouble(PyEval_CallObject(f_call, arglist));  /* x2 = f[i](x,t) */

            /* Euler computation: */
            x2 = x + delta_t * x2;

            /* set result[i] = x2 :*/
            PyList_SetItem(newx, i, Py_BuildValue("d", x2));
        }

        /* cleaning pointers: */
        Py_XDECREF(f_array);
        Py_XDECREF(f_call);
        Py_XDECREF(x_array);
        /* not doing it on newx since we're returning it */
        Py_XDECREF(xObject);
        Py_XDECREF(arglist);

        return newx;
    }
    error:
        /* cleaning memory: */
        Py_XDECREF(f_array);
        Py_XDECREF(f_call);
        Py_XDECREF(x_array);
        Py_XDECREF(newx);
        Py_XDECREF(xObject);
        Py_XDECREF(arglist);

        return NULL;        /* this raises exception */
}

/***** 3) Registration table  **************/
/* this gives the association between the function NAMES (to be used
in Python namespace) to function pointers (addresses) giving the C implementation */

static struct PyMethodDef Euler_methods[] = {
    {"step", step, METH_VARARGS, "Implements Euler solver."},
    /* method name, C func ptr., method of calling fct, method doc */
    {NULL, NULL, 0, NULL}                    /* end of table marker */
};

/***** 4) Initialization function for the module *****/
/* this is called the first time the C module is imported in Python;
it builds the Module's attribute dictionary from the registration
table and create an entry for the C module on the sys.module table */
/* name of fct matters: should be init<moduleName> */
void initEuler() {
    (void) Py_InitModule("Euler", Euler_methods);   /* module name, table ptr */
}

B.2  setup.py: python script to build extension module Euler (ref. section 3.1.5)

from distutils.core import setup, Extension

module1 = Extension('Euler', sources = ['Euler.c'])

setup (name = 'Euler Solver',
       version = '1.0',
       description = 'Euler Solver',
       ext_modules = [module1])

B.3  test.py: a test script for the Euler module

# This is a test script for the Euler solver written in C
# To use it, simply do
# >> import test
# >> test.run()

# defining the derivative functions:
def f1(x,t):
    return x[0]*x[1]

def f2(x,t):
    return -x[0]*t

L = [f1,f2] # this is the array of derivatives (represents the ODE)

# setting the initial conditions:
x = [1,2]
t = 1.0
dt = 5.0

# importing the extension module:
import Euler # we need Euler.so or Euler.pyd in the path...

# the testing routine:
def run():
    newx = Euler.step(L,x,t,dt) # testing Euler extension!
    print "Euler.step(f,x,t,dt) with f = [x[0]*x[1], -x[0]*t], \
        x =", x, "t =", t, "and dt =", dt
    print "yields newx = ", newx
    print "This should be [11.0, -3.0] if everything works."

C  Source files related to the Euler solver written in Fortran

C.1  Euler.f: source code for the Euler solver written in Fortran

C       Written by Simon Lacoste-Julien 02/7/29
C       This is a test to implement Euler function with call-back
C       mechanism in Fortran, to extend Python...
C
C EULER(N, fun, x, t, dt, newx)
C   - N is size of array input (x)
C   - fun is the ODE function - it's really an array of function; but this will
C     be taken care of in the python code: it will be called as follows
C        fun(N,i,x,t, result) where i will be the index of which function
c           to evaluate (i starts at 0 in Python)
C        x will be vector
C        t will be scalar
C        result will hold result
C   - x is state vector
C   - t is time (double)
C   - dt is time step (double)
C
        SUBROUTINE EULER(N, fun, x, t, dt, newx)
            EXTERNAL fun
            INTEGER N
            DOUBLE PRECISION x(N), t, dt, newx(N)
            DOUBLE PRECISION res

            DO 10 i = 1, N
               CALL fun(N,i-1,x,t, res)
               newx(i) = x(i) + dt * res
 10         CONTINUE

        END SUBROUTINE EULER
 

C.2  fEuler.c: source code for the C wrapper for fEuler.f

/* Written by Simon Lacoste-Julien on 2002/07/29 */
/* This is a test for Python Extension and fortran linking
I'm implementing the Euler solver used in the 522 Simulator, but using fortran code
Interface of use:
>>> import fEuler
>>> newx = fEuler.step(f,x,t,delta_t) # almost as before!
The difference is that f is no more an array of (derivative) functions, but it
is called as f(i,x,t) where i will be the index of the derivative required...
x is initial array of values,
t is time, and delta_t is step length (both are double).  And the return value
is a simple double.
*/

#include <Python.h>

/* SUBROUTINE EULER(N,fun(N,i,x,t),x,t,dt,newx) imported from Fortran: */
extern void euler_(int* N, void (*fun)(int* N1, int* i, double x1[], double* t1, double* res),
            double x[], double* t, double* dt, double newx[]);
/* notes: 1) euler_ is in lower-case and have added '_' for the name compatibility with Fortran
          2) all the arguments are passed as pointers, since Fortran pass arguments by reference
      and not by value... */

/***** Callback function interface between C and python for the derivative function: ************/
/* CALL-BACK POINTER: */
static PyObject* funPython = NULL;  /* this will hold the function passed by Python; it is GLOBAL! */

/* CALL-BACK FUNCTION WRAPPER: */
static void callBack(int* N, int* i, double x[], double* t, double* res) {
  PyObject* arglist = NULL;
  PyObject* x_python = NULL;
  int j;

  /* arguments translation: */
  x_python = PyList_New(*N);
  Py_INCREF(x_python);
  for(j = 0; j < (*N); j++) {
    PyList_SetItem(x_python, j, Py_BuildValue("d", x[j]));
  }

  arglist = Py_BuildValue("(i,O,d)", *i, x_python, *t);    /* building (i,x,t)  NOTE: x is a vector! */
  Py_INCREF(arglist);
  *res = PyFloat_AsDouble(PyEval_CallObject(funPython, arglist)); /* *res = (double) funPython(i,x,t) */

  /* memory cleaning: */
  Py_DECREF(arglist);
  Py_DECREF(x_python);

  return;
}


/***** 2) Module method definition ********/
static PyObject*
step(PyObject* self, PyObject* args) {

  PyObject* x_python = NULL;        /* will hold array x in Python */
  PyObject* newx = NULL;          /* will hold result (array) in Python */
  PyObject* xObject = NULL;       /* will hold Python x[i] */
  double* x_array;                /* array to hold values of x in C format */
  double* workArray;              /* array to hold new values of x in Fortran call  */
  int length;             /* size of x */
  int i;                  /* array index */
  double t, dt;

  /* getting arguments from python: */
  if (PyArg_ParseTuple(args, "OOdd",
               &funPython, &x_python, &t, &dt) == 0) {       /* converts Python -> C */
    /* setting error message: */
    PyErr_SetString(PyExc_RuntimeError, "Wasn't able to parse args for Euler.step fct.");
    return NULL;        /* when returning null, you raise an exception */
  } else {
    Py_INCREF(x_python);
    Py_INCREF(funPython);
    if(PyList_Check(x_python) == 0) {
      PyErr_SetString(PyExc_TypeError, "Error: x wasn't a list!");
      goto error;
    }

    /************ converting Python to C/Fortran: **************/
    length = PySequence_Length(x_python);
    x_array = calloc(length, sizeof(double));    /* x array in C */
    workArray = calloc(length, sizeof(double));  /* to be used by Fortran code */

    /* initializing x_array: */
    for(i = 0; i < length; i++) {
      Py_XDECREF(xObject);
      xObject = PyList_GetItem(x_python, i);
      Py_INCREF(xObject);
      x_array[i] = PyFloat_AsDouble(xObject);
    }

    /* CALLING Fortran code: */
    euler_(&length, callBack, x_array, &t, &dt, workArray);

    /************* translating result to Python *****************/

    newx = PyList_New(length);              /* initializing the return array */
    Py_INCREF(newx);
    for(i = 0; i < length; i++) {
      PyList_SetItem(newx, i, Py_BuildValue("d", workArray[i]));
    }

    /* cleaning memory: */
    Py_XDECREF(x_python);
    /* not doing it on newx since we're returning it */
    Py_XDECREF(xObject);
    Py_XDECREF(funPython);
    free(x_array);
    free(workArray);

    return newx;
  }
 error:
  /* cleaning memory: */
  Py_XDECREF(x_python);
  Py_XDECREF(newx);
  Py_XDECREF(xObject);
  Py_XDECREF(funPython);

  return NULL;        /* this raises exception */
}

/***** 3) Registration table  **************/
/* this gives the association between the function NAMES (to be used
   in Python namespace) to function pointers (addresses) giving the C implementation */

static struct PyMethodDef fEuler_methods[] = {
  {"step", step, METH_VARARGS, "Call Fortran Euler Solver."},
  /* method name, C func ptr., method of calling fct, method doc */
  {NULL, NULL, 0, NULL}                    /* end of table marker */
};

/***** 4) Initialization function for the module *****/
/* this is called the first time the C module is imported in Python;
   it builds the Module's attribute dictionary from the registration
   table and create an entry for the C module on the sys.module table */
/* name of fct matters: should be init<moduleName> */
void initfEuler() {
  (void) Py_InitModule("fEuler", fEuler_methods);   /* module name, table ptr */
}

C.3  setup_fEuler.py: python script to build extension module fEuler

# written by SLJ 2002-8-28
# took model from http://www.python.org/doc/current/ext/building-on-unix.html
# To use it on my Windows XP (with mingw gcc compiler):
# python setup.py build --compiler=mingw32
# It uses the library libEuler.a

from distutils.core import setup, Extension

module1 = Extension('fEuler',
                    libraries = ['Euler'],
                    library_dirs = ['./'],
                    sources = ['fEuler.c'])

setup (name = 'fEuler',
       version = '1.0',
       description = 'This a Euler solver implemented in Fortran',
       author = 'Simon Lacoste-Julien',
       author_email = 'simonlacoste@videotron.ca',
       url = 'http://moncs.cs.mcgill.ca/MSDL/people/slacoste/',
       long_description = '''
This is really a Euler solver implemented in Fortran!
''',
       ext_modules = [module1])

C.4  test_fEuler.py: a test script for the fEuler module

# This is a test script for the Euler solver written in Fortran
# To use it, simply do
# >> import test_fEuler
# >> test_fEuler.run()

# defining the derivative functions:
def f1(x,t):
    return x[0]*x[1]

def f2(x,t):
    return -x[0]*t

L = [f1,f2] # this is the array of derivatives (represents the ODE)

# importing the extension module:
import fEuler # we need fEuler.so or fEuler.pyd in the path...

# we define here a routine similar to the one that will be
# used from the 522 Simulator
def step(f,x,t,dt):
    # we need to change the way f is called (to be consistent with
    # new interface of the solver
    def fun(i,x1,t1):
        return f[i](x1,t1)
    return fEuler.step(fun,x,t,dt)

# setting the initial conditions:
x = [1,2]
t = 1.0
dt = 5.0

# the testing routine:
def run():
    newx = step(L,x,t,dt) # testing fEuler extension!
    print "step(f,x,t,dt) with f = [x[0]*x[1], -x[0]*t], x =", x, "t =", t, \
          "and dt =", dt
    print "yields newx = ", newx
    print "This should be [11.0, -3.0] if everything works."

D  Modified code for ddassl.pyf

!%f90 -*- f90 -*-
! SLJ: 2002-08-05: added definitions of res & jac plus intent everywhere.
! SLJ: 2002-08-15: added the dimension dependance on len(y)... [CAN'T do that for call-back]
! SLJ: 2002-08-15 20h20: started tests...
!         22h45: done with JAC and RES...

! will take the convention that:
! - ipar[0] = dimension of y
! - ipar[1] = dimension of ipar (>= 3)
! - ipar[2] = dimension of rpar

python module ddassl__user__routines
    interface ddassl_user_interface
        subroutine res(t,y,yprime,delta,ires,rpar,ipar)
        double precision intent(in) :: t
            double precision dimension(ipar[0]), depend(ipar), intent(in) :: y
            double precision dimension(ipar[0]), depend(ipar), intent(in) :: yprime
            double precision dimension(ipar[0]), depend(ipar), intent(inout) :: delta
            integer intent(inout) :: ires
            double precision dimension(ipar[2]), depend(ipar), intent(in) :: rpar
            integer dimension(ipar[1]), intent(in) :: ipar
        end subroutine res
        subroutine jac(t,y,yprime,pd,cj,rpar,ipar)
            double precision intent(in) :: t
            double precision dimension(ipar[0]), depend(ipar), intent(in) :: y
            double precision dimension(ipar[0]), depend(ipar), intent(in) :: yprime
            double precision dimension(ipar[0],ipar[0]), depend(ipar), intent(inout) :: pd
            double precision intent(in) :: cj
            double precision dimension(ipar[2]), depend(ipar) intent(in) :: rpar
            integer dimension(ipar[1]), intent(in) :: ipar
        end subroutine jac
    end interface ddassl_user_interface
end python module ddassl__user__routines
python module ddassl ! in
    interface  ! in :ddassl
        subroutine ddassl(res,neq,t,y,yprime,tout,info,rtol,atol,idid,rwork,lrw,iwork,liw,
            rpar,ipar,jac) ! in :ddassl:ddassl.f
            use ddassl__user__routines
            external res
            integer intent(in) :: neq
            double precision intent(inout) :: t
            double precision dimension(*), intent(inout) :: y
            double precision dimension(*), intent(inout) :: yprime
            double precision intent(in) :: tout
            integer dimension(15), intent(in) :: info
            double precision dimension(*), intent(inout) :: rtol
            double precision dimension(*), intent(inout) :: atol
            integer intent(inout) :: idid
            double precision dimension(*), intent(inout) :: rwork
            integer, intent(in) :: lrw
            integer dimension(*), intent(inout) :: iwork
            integer, intent(in) :: liw
            double precision dimension(*), intent(in) :: rpar
            integer dimension(*), intent(in) :: ipar
            external jac
        end subroutine ddassl
    end interface
end python module ddassl

! This file was auto-generated with f2py (version:2.13.175-1250).
! See http://cens.ioc.ee/projects/f2py2e/

E  Content of the HybridSystems directory

• File HybridSystems_ER_mdl.py: contains the model of the HS formalism.

• File HybridSystems.py: contains the executable code for the modelling environment of the HS formalism.

• File HS_atom3TO522Simulator_GG_mdl.py: contains the graph grammar model to generate code for the 522 Simulator.

• File HS_atom3TO522Simulator.py: contains the executable code for the graph grammar to generate code for the 522 Simulator.

• File codeGenTools.py: a module used by the code generation algorithm.

• The other files in the directory are used by AToM³ for the modelling environment and the graph grammar.

• Directory Simulator522 - contains the files for the 522 simulator;
  • Directory Solvers_buid: contains the source code to rebuild the new solvers

  • files .so: those are the new solvers already built for Unix

  • files .pyd: those are the new solvers already built for Windows.

• Directory Simulator_minimal: contains the files for the 522 Simulator without the new solvers.

• Directory models_HS: contains examples of models in the HS formalism.

• Directory generated_models: contains examples of generated code for models in the HS formalism.

• Directory code_for_codeGen: contains archive of source code for various components of the modelling environment and the graph grammar:
  • button_copy.py Code for the button 'Copy'.

  • button_gen522Code.py Code for the button 'Generate Code'.

  • button_Start_Simulator.py Code for the button 'Launch Simulator'.

  • gen522Code.py Code for the code generation part of the code generation graph grammar.

References

[1]

Zeigler, Bernard P. et al., Theory of Modeling and Simulation, Academic Press, San Diego, 2000, QA76.9 C65 Z44 2000.

[2]

Lacoste-Julien, Simon, MSDL Talk: Intro to Hybrid Systems Modelling, School of Computer Science, McGill University, Montréal, Québec, Canada, July 2002, web site: http://moncs.cs.mcgill.ca/people/slacoste/research/talk.html

[3]

Lutz, Mark, Programming Python, O'Reilly, Sebastopol, 2001, QA76.73 P98 L881 2001.

[4]

de Lara, Juan et al., Using Meta-Modelling and Graph Grammars to create Modelling Environments, Electronics Notes in Theoretical Computer Science, 72 No.3 (2002).
link: http://www.cs.mcgill.ca/ hv/publications/02.ICGT.GT-VMT.pdf

Footnotes:

¹This is a Hybrid Systems simulator written in Python by Olivier Dubois and Eric McSween as a project during the course 308-522A in Fall 2001. See http://www.cs.mcgill.ca/ emcswe/cs522/project/index.html for their web page.

²DASSL: Differential Algebraic Systems Solver Library

³The first 'd' is to indicate that the code is using double precision numbers; sdassl.f will represent the single precision version.

⁴ This has been a major source of frustration for me during this project, as the basic commands to use f2py described in the user guide were completely obsolete. The general information about f2py given in the user guide is still correct [I think], though.

⁵If we want to be rigorous, it could happen that this chain of meta-levels would never end (and this is reminiscent of meta-theory chains that happen also in the rigorous definition of theories in logic...).

⁶developed at Universidad Autónoma de Madrid

⁷``Event_HS'' was used instead of ``Event'' for the name of this entity in order to avoid name conflict with Tkinter (used internally by AToM³).

⁸Notice the AToM³ naming convention: the name of the model being built (HybridSystems) followed by the formalism in which it is built (ER) followed by mdl.py.

⁹Here, entity means simply ``object which can be added in the model''. So it can be either an entity from ER or a relationship from ER. It is unfortunate that there can be some ambiguity in the terminology of AToM<sup>3</sup> since we have the button ``Edit entity'' which can be used to edit also relationships... It should be clear from the context whether I use entity to mean entity from ER or general entity.

¹⁰Note that AToM³displays lists (and tuple) as strings separated by |. For example, (g, 9.8) is displayed as g|9.8|.

¹¹The naming convention for a graph grammar is <formalism>_<transformation name>_GG_mdl.py.