HybridSystems

  1. Overview
  2. The elements of a hybrid system
    1. Global and local time
      1. Interpreting and converting between time frames
    2. Specifying the initial state
  3. Interfacing with the Model class
    1. Trajectories
    2. Important model methods
  4. Determining the initial sub-model
    1. Example
    2. More complex hierarchies of hybrid models
    3. Some event consistency checking comes for free
  5. Hybrid systems with only one vector field
  6. Using hybrid systems to mix maps, preset trajectories, and ODEs
    1. Using hybrid systems to implement discrete delays
  7. Using hybrid systems to represent numerically-computed maps

1. Overview

On this page we describe the use of the Model class to represent true hybrid dynamical systems, that is, dynamical systems with smooth flows that are punctuated by discrete events (optionally undergoing discrete mappings of state at these events). The vector fields either side of an event may have different dimensions. Other uses for the Model class are described in the UserDocumentation.

The focus of PyDSTool is on applications in the continuous domain, both in time and state. As a result, the approach taken to hybrid systems differs from that in many discrete event-based simulators such as [WWW] SimPy. We effectively treat the continuous parts of trajectories as more fundamental than the discrete events that may switch between different vector fields or map system states instantaneously.

PyDSTool is not an efficient general-purpose discrete event simulator (both in terms of computation and specification), but is well suited when the main focus of the model is continuous trajectories such as those defined by ODEs.

2. The elements of a hybrid system

There are three essential elements to a hybrid system model in PyDSTool: model sub-systems, events, and transition rules between the sub-systems that are applied on occurrence of the events. The specification of an initial state condition must uniquely determine which sub-model should be selected to begin determining a trajectory. This is done by evaluation of the associated events for the system at the initial condition.

A terminal event (see Events) may occur that stops the trajectory generation. At this point the transition rules for the stopped sub-model are consulted, and a new sub-model is chosen by applying those rules to the final state condition. The final state may also be mapped to a new value before becoming the initial condition for the next sub-model.

2.1. Global and local time

When used as a complete model in its own right, all references to time in a non-hybrid NonHybridModel object are to the absolute ("global") time. However, when embedded in a HybridModel object as a hybrid system, the concept of absolute and relative ("local") time will become distinct for certain purposes. Generator, HybridModel and Trajectory objects have an attribute globalt0 for these purposes.

A trajectory may represent a stereotypical behavior that is repeated often in the dynamics of a hybrid model. Rather than re-computing the same trajectory for different absolute time intervals, an independently existing trajectory whose time domain is stated in local time is more easily re-used. To commit a copy of such a model to a particular absolute time interval in a hybrid model, the globalt0 is merely set appropriately. The model classes make the handling of absolute time transparent to the end-user, and knows how to adjust queries in absolute time to the relative time of an appropriate constituent trajectory segment.

The built-in auxiliary function globalindepvar provides access to a hybrid system's absolute ("global") value of the independent variable (usually time) from within an individual model. More information about this function is provided on the page UserDocumentation.

2.1.1. Interpreting and converting between time frames

The interpretation of global time is made according to this assumption. "Global time" always refers to an external time frame, and "local time" for an object is typically set so that it works on a time interval [0, t_end] that may be independent of how that interval is interpreted external to the object. The rule is that the globalt0 attribute is itself in the global time frame. This means that a global time is converted to a local time by subtracting globalt0, and a local time is converted to a global time by adding that quantity.

2.2. Specifying the initial state

As of version 0.87 defining hybrid systems is now more systematic and robust. It used to be possible to specify to PyDSTool a sub-model that could be entered with an initial state that violates one or more of its terminal event conditions. Such an eventuality is now checked at initialization, avoiding simulation failure. This is achieved using formal boundary containment conditions, made possible using classes from MProject.py, particularly ModelInterface and Condition. Section 4 goes into detail about how this is now done.

3. Interfacing with the Model class

Some of this information is repeated from the UserDocumentation page, for completeness.

The Model and HybridModel classes serve two purposes. The first is to generalize the forms of dynamical system that can be simulated in PyDSTool, by supporting hybrid dynamical systems. The second is to provide additional utilities for model introspection and manipulation.

A model instance holds information that determines which sub-model is used to calculate a trajectory segment, given the initial state. When a terminal event occurs during the calculation, the model uses a rule to determine which sub-model to use next.

3.1. Trajectories

Hybrid trajectories are stored inside the Model object they are associated with. (The implementation of Trajectory objects is described on the page TechDocumentation.) Hybrid trajectories can be accessed much the same way as regular trajectory objects when treated as continuous curves. However, hybrid trajectories can also be treated as discrete mappings, from one "terminal event" to the next. A hybrid trajectory contains information about the event parameters used during its calculation. It can be referenced as my_model.trajectories[<traj_name>].

To access a hybrid trajectory as a map, we call the Model instance that contains it with the additional option 'asmap=True': my_model(<trajname>, <integer>, asmap=True), where the integer represents the number of the partitions to access. These integers may range from 0 to the number of partitions minus one. The number of partitions can be found by calling my_model.trajectories[<traj_name>].numPartitions.

The event parameters are stored in a copy of the EventStruct object associated with each sub-model used to create the hybrid trajectory in a dictionary attribute 'modelEventStructs', accessed in the following way: my_model.trajectories[<traj_name>].modelEventStructs[<gen_name>].

3.2. Important model methods

When you have built a model, use help(Model.<methodname>) to read the documentation string for any of these methods.

The query method is particularly important for large models built using the ModelSpec syntax, where there are hierarchically-defined components with composite names (e.g. 'cell2.syn.s') separated by dots. With this method the user can return information about the Model using this syntax for names. (Internally, these names are converted to strings that use '_' for separators instead of dots.)

4. Determining the initial sub-model

The approach taken in PyDSTool is that when a hybrid system is declared in which there is more than one vector field (i.e. more than one sub-model), the hybrid solver must be able to automatically determine the correct sub-model to use at the initial time for a given initial condition (IC). There need to be a set of self-consistent conditions placed on the initial state in order for the hybrid solver to determine which sub-model to use at the starting time -- the user cannot explicitly name the initial sub-model.

Primarily, these conditions come from evaluating the terminal events associated with each sub-model at the initial condition provided. The HybridModel class will only evaluate the terminal events that have been set active (which is their default state if the user has not specified). The reason the user is not allowed to suggest an initial sub-model is twofold. Firstly, the user might not specify a sub-model that is consistent with the way the terminal events partition the state space. Therefore, regardless of how the initial sub-model is determined, this consistency check ought to be made. Secondly, if this consistency check remains incomplete because the initial sub-model choice is under-specified, then we are at the mercy of human input to initiate hybrid solving. How will the hybrid solver know what to do if it is running unsupervised, and the initial condition is generated automatically by an algorithm that doesn't know about the hybrid system definition that will be using that IC? Forcing the setup of the hybrid system to be logically self-consistent without the arbitrary input of a user to set the initial sub-model guarantees some safety in running the hybrid solver unsupervised.

How does this problem arise in practice? For one thing, directionless threshold-crossing events (dircode=0) do not partition the state space in a way that helps to determine in which sub-domain an initial condition lies. Additionally, even for a set of events that all have a fixed direction of threshold crossing, the set still may not uniquely determine a sub-domain. This may happen for several reasons: (1) because it is natural to place ICs right on an event surface for some sub-model "G", leaving the eligibility of "G" indeterminate; (2) because it is most efficient to specify only sufficient conditions for termination of a trajectory computation with the current vector field: an exhaustive set of necessary and sufficient termination conditions would be required to guarantee that an IC is uniquely placed in a sub-domain of the state space; (3) even if all those conditions were listed, some terminal events are set to occur not when a function of the system state passes a threshold, but when a certain amount of time has elapsed.

If the initial sub-model is under-specified then you will get an error saying that the system has found too many eligible sub-models to begin computing a solution. If case (1) holds then the only chance for the system to determine the correct initial sub-model is if all the other sub-model can be eliminated for the IC. This can never happen if there are non-state events (i.e. when case (3) holds). Artificially introduced "indicator" variables are a "trick" solution to this problem that avoids forcing terminal events to be specified in an unnatural way. Version 0.87 allows directions to be specified as part of the boundary definitions so that initial conditions can include derivative information that point in to or out of the event boundaries, enabling resolution of this problem.

Indicator variables are artificial state variables (usually integers, e.g. 0 and 1) that you force to lie in a singular domain (i.e. the singleton set containing that integer), so that you get additional domain constraints without adding terminal events to your system. This information usually represents something about the system's recent history, in terms of prior terminal events. Note that although adding these variables artificially means that they are ancillary to the main variables, you should not make these critical indicator variables "auxiliary" variables (either explicitly using the auxvars initialization key or by omission from a vars initialization key), otherwise the values will not be tested and updated properly. However, the Model formalism has a secondary way to mark certain variables as ancillary: the forceIntVars(<var_list>) method labels variables as "internal", while the remaining variables are "external", i.e. externally observable. This will have more significance when performing Model estimation (described later).

In order to use an indicator variable, simply set its value as part of the initial conditions: only the sub-models for which that value corresponds to the singular domain set for that sub-model will match, thereby providing you sufficient conditions for determining a unique initial sub-model.

4.1. Example

In PyDSTool/tests/IF_model_test.py, an integrate-and-fire neuron model with a square-pulse spike is implemented as a hybrid system with two sub-systems: 'leak' and 'spike'. A spike threshold event occurs when the variable V increases through a value of -60. The sub-threshold 'leak' ODE system is then replaced by the 'spike' explicit function system that terminates after a set duration given by the parameter 'splen'. This terminal event specification means that an initial condition cannot be placed uniquely in one of the two sub-systems. The excited variable has been introduced to record the most recent history of whether spike threshold was passed.

Here is the specification for this model. 

all_model_names = ['leak', 'spike']

# 'excited' is an internal variable of the model, and is used to
# ensure that the.compute() method can determine which Generator
# to start the calculation with
leak_event_args = {'name': 'threshold',
                   'eventtol': 1e-3,
                   'eventdelay': 1e-5,
                   'starttime': 0,
                   'active': True,
                   'term': True,
                   'precise': True}
leak_thresh_ev = Events.makePythonStateZeroCrossEvent('V', -60, 1,
                                                 leak_event_args)
leak_args = {'pars': {'I': 1.3, 'gl': 0.1, 'vl': -67},
              'xdomain': {'V': [-100,50], 'excited': 0},
              'xtype': {'excited': int},
              'varspecs': {'V':       "I - gl*(V-vl)",
                              'excited': "0"
                              },
              'algparams': {'init_step': 0.02},
              'events': leak_thresh_ev,
              'abseps': 1.e-7,
              'name': 'leak'}

ics = {'V': -80., 'excited': 0.}

DS_leak = embed(Generator.Vode_ODEsystem(leak_args), icdict=ics, tdata=[0, 30])
feat = binary_feature('is_excited', pars={'varname':'excited'})
DS_l_condition = condition({feat: False})
DS_leak_MI = intModelInterface(DS_leak, DS_l_condition)

# spike length parameter 'splen' must be contained within 'tdomain' in
# order to get a fully-formed square-pulse `spike`
spike_args = {'tdomain': [0.0, 1.5],
                'varspecs': {'V': "if(t<splen,50,-95)", 'excited': "1."},
                'pars': {'splen': 0.75},
                'xdomain': {'V': [-97, 51], 'excited': 1},
                'xtype': {'excited': int},
                'name': 'spike'}
DS_spike = embed(Generator.ExplicitFnGen(spike_args), icdict=ics, tdata=[0, 30])

DS_s_condition = condition({feat: True})
DS_spike_MI = intModelInterface(DS_spike, DS_s_condition)

# test discrete state mapping that is used at changes of vector field
# after terminal events
# must set excited to 0 in order to meet bounds requirements of IF gen
epmapping = EvMapping({"xdict['V']": "xdict['V']+15",
                       "xdict['excited']": "0"})

# build model object from individual sub-models
DS_leak_info = makeModelInfoEntry(DS_leak_MI, all_model_names,
                                  [('threshold', 'spike')])
DS_spike_info = makeModelInfoEntry(DS_spike_MI, all_model_names,
                                 [('time', ('leak', epmapping))])
modelInfoDict = makeModelInfo([DS_leak_info, DS_spike_info])

IFmodel = Model.HybridModel({'name': 'IF_fit', 'modelInfo': modelInfoDict})

Here, the excited variable is set to exist in the singleton domain {0} for the leak phase, and {1} for the spike phase. The initial condition icdict = {'V': -80., 'excited': 0}  will only satisfy the domain test for the excited variable in the leak phase.

An example of case (2) occurs in PyDSTool/tests/SLIP2D.py, in which a spring loaded inverted pendulum model of a leg in 2D is implemented. The hybrid model consists of two phases, a 'stance' and a 'flight' phase. The flight phase is entered when the leg reaches full extension (zeta=1). The initial velocities of the state variables in the file are pointing upwards, but the events are only given in terms of position variables. It is also natural to specify an IC for the system where zeta=1, which is on the event boundary for the stance phase Generator. So while it is intuitively obvious that the flight phase should be the initial Generator in this situation, the problem is formally underspecified. A very small number could be added to the position variables such that zeta would be microscopically greater than 1, but this is inelegant. The hybrid solver can resolve the situation if a variable incontact is introduced, and used with the binary_feature class. This is implemented almost identically to the excited variable in the Integrate-and-Fire model.

4.2. More complex hierarchies of hybrid models

Although each of the two sub-models in the Integrate-and-Fire example above consists of a single Generator object embedded into a NonHybridModel object, we could in principle build hybrid models out of components that are themselves hybrid models. Such is the flexibility of this specification scheme, which comes at the expense of a bit of additional syntactic baggage.

4.3. Some event consistency checking comes for free

As a final note, if you make a mistake in the setting up of the events, such that a given IC does not meet any of the pre-terminal event criteria specified, then you will get an error saying that there were no sub-models eligible for beginning the trajectory computation.

5. Hybrid systems with only one vector field

Some hybrid systems involve a single vector field with events that cause discrete state changes. In PyDSTool such events are implemented as terminal events that have associated state mappings to make the discrete changes. These events must either point back to the same vector field (using the associated sub-model's name) or the name 'terminate' to end the trajectory solution altogether.

6. Using hybrid systems to mix maps, preset trajectories, and ODEs

The event transition functions of a hybrid system Model object can be used to represent sophisticated forms of mapping that can be interspersed with ODE vector fields. By the flexible class hierarchy of Generators, a mixture of types of Generator can be present in a hybrid system. For instance, part of a vector field can be specified by an ODE, and another by a preset time-course (e.g., using an InterpolateTable Generator), such as that in /tests/IF_squarespike_model.py, in which a threshold event in an ODE model of an excitable medium triggers a "square-pulse" spike of pre-determined duration to occur, before returning control over the trajectory computation to the ODE. A related example is shown above.

6.1. Using hybrid systems to implement discrete delays

One application of event manipulation by transition functions is to implement a form of discrete delay in the coupling between ODE variables. This is currently in development, and so is not a feature that is easy to use.

The delay of an event can be used to effectively implement discretely-delayed dirac delta function pulses (where the pulses correspond to discrete events). For instance, see the delta-function pulse-coupling between Integrate-and-Fire neurons in /tests/IF_delaynet.py. Note that the effect of the pulse (event) need not merely be an instantaneous change in state. The event can be used to trigger longer-lasting changes in state (e.g. see /tests/IF_delaynet_syn.py).

7. Using hybrid systems to represent numerically-computed maps

Consider the example of a 2D spring-loaded inverted pendulum (SLIP) model of a leg in /tests/SLIP_2D_maps.py. This particular mechanical model has an unstable limit cycle, which is identified using a simple shooting method and a SciPy wrapper around Minpack routines in /tests/SLIP_2D_pdc.py. Of interest in the analysis of this model is the map from one foot-fall impact to the next, which in general can only be computed numerically. However, the power of dynamic typing and the interactive nature of Python means that a Python function can be written that behaves exactly like the required mapping, although in fact it involves a numerical integration of the hybrid system until the appropriate event occurs.

This is an example of the power of Python's dynamic typing and the PyDSTool classes in order to embed numerical algorithms inside each other.

last edited 2009-02-18 02:10:30 by RobClewley