Skip to main content

Structured Transition Guards

The definitions of hybrid automata that I have come across define the transition guards as subsets of the continouous state space. Here I propose a more structured definition; one suitable for computation using differential-algebraic equation system solvers with root-finding capabilities, like IDA from [Hindmarsh2005].

A miniature HA

We start from the definition in [Schaft2000]:

\begin{equation*} \begin{gathered} H = (L,X,A,W,E,Inv,Act) \\ e = (l,a,Guard_{ll^\prime},Jump_{ll^\prime},l^\prime) \\ Guard_{ll^\prime} \subset X \end{gathered} \end{equation*}

Only some of these components are useful here. We throw out unneeded pieces: symbols (\(A\), \(a\)), invariants (\(Inv\)), and jumps (\(Jump_{ll^\prime}\)). To add insult to injury we rename its activity-function (\(Act\)) to \(\mathbf{F}\), and guard (\(Guard_{ll^\prime}\)) to \(g\). Lets assume that the DAEs are regular 1 for good measure. A HA with \(M\) discrete variables and \(N\) continuous variables is then:

\begin{equation*} \begin{gathered} H^M_N = (L,X,E,\mathbf{F}) \\ L = \Z^M \\ X = \R^N \times \R^N \times \R\\ E \subset L \times \mathcal{P}(X) \times L \\ \mathbf{F} \colon L \to X \to \R^N \end{gathered} \end{equation*}

Simulating HAs means advancing \(?(P,x) \in L \times X\) through interleavings of continuous and discrete problems. In the continuous problem we are given an initial continuous state, an activity and a set of guards. After advancing the solution until one of the guards become active, we are rewarded with a discrete problem. In the discrete problem we are given an initial discrete state, zero or more symbols (corresponding f.ex. to an active guard), and compute the next discrete state in which a new continuous problem is waiting for us. The scope of this post streches from the start of solving a continuous problem all the way up to and including, detecting that a guard is active.

1

Not overdetermined, nor underdetermined.

Logical expressions and DFAs

We shift our focus to the booleans, \(\mathbb{B} = \{ \bot, \top \}\), and expressions over boolean variables \(?b \in \mathbb{B}^N\). These expressions, \(g \colon \mathbb{B}^N \to \mathbb{B}\), involve the usual suspects; \(\land\), \(\lor\), \(\neg\). We intend to use these expressions to represent guards, which let us know when a continuous problem is solved. To this end we construct a DFA that traces changes in the truth value of \(g\) caused by changes in \(?b\). We begin from the definition of a DFA in [Hopcroft2006]:

\begin{equation*} \begin{gathered} D = (Q,\Sigma,\delta,q_0,F) \\ \delta \colon Q \times \Sigma \to Q \end{gathered} \tag{automaton} \end{equation*}

We ignore the initial state \(q_0\), and instead use a variable \(?q \in Q\), separate from the definition of the DFA (\(?q^0\) being analogous to \(q_0\)). We let \(\top_i\) denote a change of \(?b^{t}_i = \bot\) to \(?b^{t+1}_i = \top\), and vice versa from \(\bot\) to \(\top\). These symbols define the alphabet:

\begin{equation*} \Sigma = \displaystyle \bigcup_{i \in N} \{ \top_i, \bot_i \} \tag{alphabet} \end{equation*}

The states of the DFA, \(Q\), correspond to the elements of \(\mathbb{B}^N\):

\begin{equation*} Q = \{ q_0, \ldots, q_{2^N - 1} \} \tag{states} \end{equation*}

\(q_i\) is mapped to \(b \in \mathbb{B}^N\) by representing \(i\) in base 2 and identifying \(b_j\) with the \(j\)th bit of \(i\) (\(0 \to \bot\), \(1 \to \top\)). This mapping is denoted \(h\):

\begin{equation*} h \colon Q \to \mathbb{B}^N \end{equation*}

The terminals are the states for whose valuations the guard evaluates to true.

\begin{equation*} F = \{ q \in Q \mid (g \circ h)(q) = \top \} \tag{terminals} \end{equation*}

The transition function, \(\delta\), is defined in two parts. \(\delta_{move}\) reacts to changes in valuation, while \(\delta_{stay}\) absorb valuations already present in the current state. In what follows the specific type of change is abstracted away by \(\Delta_j \in \{ \top_j, \bot_j \}\). Note that the automaton is deterministic by construction:

\begin{equation*} \begin{aligned} \delta_{move}(q_i,\Delta_j) &= q_{i \oplus 2^j} \\ \delta_{stay}(q_i,\Delta_j) &= q_i \end{aligned} \end{equation*}
\begin{equation*} \delta(q_i,\Delta_j) = \begin{cases} h(q_i)_j = \Delta \to \delta_{move}(q_i,\Delta_j) \\ h(q_i)_j \neq \Delta \to \delta_{stay}(q_i,\Delta_j) \end{cases} \tag{transitions} \end{equation*}

Evaluating \(g(b)\) can be done by checking if the DFA accepts the symbol sequence \(\{ (?b_j)_j \ldots \}\). The truth-value of \(g\) is equivalent to the terminality of the final valuation \(?q = q_f\). Due to the way the DFA is constructed it does not matter which state is chosen as \(?q^0\):

\begin{equation*} g(b) = \begin{cases} q_f \in T \to \top \\ q_f \notin T \to \bot\end{cases} \end{equation*}

Responding to a change in the valuation of \(b_j\), one would simply continue from \(?q\), consume the corresponding symbol and see if \(?q\) now is in an accepting state. By consuming a symbol, \(\Delta_j\) we mean computing a new valuation of \(q\), \(?q^{t+1} = \delta(?q^{t},\Delta_j)\).

Simulating structured guards

Recall that [Schaft2000] defined the guard as a subset of \(X\). We will restructure it to a form suitable for computation by defining the guard as a logical expression over a variable \(?b \in \mathbb{B}^{|G|}\), and identify \(?b\) with the truth-values of implicit algebraic inequalities \(R = \{ r(x) \}\):

\begin{equation*} b_i(x) = \begin{cases} r_i(x) < 0 \to \top \\ r_i(x) \nless 0 \to \bot \end{cases} \end{equation*}

We abuse notation and write \(g(x)\) for \(g(b_i(x) \ldots)\). The guard \(g \colon X \to \mathbb{B}\) describes the subset \(\{ x \in X \mid g(x) = \top \}\), and the transitions of \(H^M_N\) are now:

\begin{equation*} \begin{aligned} E &\subset L \times \mathbf{G} \times L \\ \mathbf{G} &= \{ X \to \mathbb{B} \} \end{aligned} \end{equation*}

Let \(R_i = \{ r(x) \}\) denote the inequalities implicit in \(g_i\), and \(R = \bigcup_{i \in |G|} R_i\) the inequalities implicit in all of \(G\). \(R\) are the root-functions provided to our root-finding solver, which in turn outputs the positive and negative roots of \(R\), \(\{ +j, -j \; | \; j \in |R| \}\) as the continuous problem is being solved. These roots correspond to changes in the truth values of \(r_j(x) < 0\), and thus to changes in \(?b_j\), and thus (!) to symbols in \(D(G)\).

We now have all the pieces to outline an algorithm for solving a continuous problem with structured guards. We construct the guard automata, and define a guard variable \(?q \in \prod_{i \in |G|} Q_i\), with \(?q^0 = \{ q_{i0} \; \mid \; i \in |G| \}\). We then compute and consume the roots implicit in the initial state, and continue feeding it the roots found during simulation. For every root \(\Delta_j\), \(q\) is advanced in lockstep \(?q^{t+1} = \delta^\ast(?q^{t},\Delta_j)\):

\begin{equation*} \delta^\ast_i(q_i,\Delta_j) = \begin{cases} \exists r_k \in R_i \mid r_k = r_j \to \delta_i(q_i,\Delta_k) \\ \nexists r_k \in R_i \mid r_k = r_j \to q_i \end{cases} \end{equation*}
\begin{equation*} \delta^\ast(q,\Delta_j) = \{ \delta^\ast_i(q_i,\Delta_j) \mid i \in |G| \} \end{equation*}

When any \(q_i\) is in an accepting state, it means a guard is active, and that the continuous problem has been solved. The algorithm is illustrated with a diagram:

G Solving continuous problems with structured guardsvb Beginv0 Construct guards and initialize variablevb->v0 v5 Compute initial rootsv0->v5 v1 Advance solutionv2 Root found?v1->v2 v2->v1 nov3 Advance guardsv2->v3 yesv4 In accepting state?v3->v4 v4->v1 nove Endv4->ve    yesv5->v4

Example system

Lest do some simple simulations. Recall that when simulating a HA, an active guard marks the end of a continuous problem and the beginning of a discrete one. Here we instead continue solving the continuous problem and advancing the guard automata, as the purpose is to illustrate the correlation between continuous roots and the change in state of guard DFAs.

To this end we define a class of hybrid automata with sinusoid dynamics, \(S(A,b) \colon \R^{+^{|A|}}_{\ast} \times \R \to \mathbf{H^1_{|A|}}\):

\begin{equation*} \begin{gathered} S(A,b) = (\{0\},X,\{(0,g,0)\},F)\\ X = \R^N \times \R^N \times \R\\ g = \bigwedge\limits_{i = 0}^{N} ?(b - x_i < 0) \\ F = \{ \dot{x}_i - sin(a_i \cdot t) \} \end{gathered} \end{equation*}

The system is modelled and simulated with hysj and IDA from [Hindmarsh2005]. The illustrations are done with matplotlib from [Hunter2007]. The plot shows the parameter, \(b\), and the continuous state, \(x\), drawn on top of the guard automaton for \(g\). \(\delta_{stay}\) is not drawn. The active state \(?q\) has a non-white facecolor, with terminals \(F\) and green, and non-terminals \(Q \setminus F\) in red. Code listing below.

\(S(\{ 1.0 \},0.5)\):

\(S(\{ 0.75, 1.25 \},0.5)\):

\(S(\{ 0.75, 1.0, 1.25\},0.5)\):

Problems

This approach is only appropriate for HAs whose guards involve few inequalities, due to the explosion of the order and size of the guard automata.

We defined the guard with respect to (strict) inequalities. A root finding solver reports roots, \(r(x) = 0\). A positive root would imply \(r(x) < 0\), but a negative root does not imply \(r(x) \nless 0\). The guard is not really active yet, and might not ever be if \(\dot{x} = 0\). In addition one might want to support nonstrict inequalities and equalities.

The formulation of guards as DFAs suggests a defininiton of non-determinism of the HA in terms of the existence of a location for which the product of the guard-DFAs are non-deterministic. This was not pursued here.

Glossary

DFA

Deterministic finite-state automaton

DAE

Differential-algebraic equation

HA

Hybrid automaton

Bibliography

Schaft2000(1,2)

"Arjan van der Schaft and Hans Schumacher", "An introduction to hybrid dynamical systems"

Hunter2007

"Hunter, J. D.","Matplotlib: A 2D graphics environment"

Hindmarsh2005(1,2)

"Hindmarsh, Alan C and Brown, Peter N and Grant, Keith E and Lee, Steven L and Serban, Radu and Shumaker, Dan E and Woodward, Carol S","SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers"

Hopcroft2006

"John E Hopcroft, Rajeev Motwani, Jeffrey D Ullman", "Introduction to Automata Theory, Languages, and Computation"

Code

import hysj
import numpy
from matplotlib import pyplot as plot

def make_problem(A,b):
  import math

  continuous = hysj.simulators.continuous

  order = 1
  N = len(A)

  ir = hysj.mathematics.Program()
  t = ir.symbol()

  X = [ ir.variable(t,order,i) for i in range(N) ]
  F = [ ir.sine(ir.multiplication([ir.constant(a),t])) for a in A ]
  E = [ ir.equality(x[-1],f) for x,f in zip(X,F) ]
  R = [ ir.strict_inequality(ir.constant(b),x[0]) for x in X ]
  G = [ ir.conjunction([ir.assertion(r) for r in R]) ]

  t_begin  = 0.0
  t_size   = 5
  t_end    = t_size * math.pi
  t_count  = t_size * 256
  t_delta  = (t_end - t_begin) / (t_count - 1)

  t_space = numpy.linspace(t_begin,t_end,num = t_count)

  problem = continuous.Problem(program = ir,
                               variables    = continuous.Variables(independent = t,
                                                                   dependent   = X),
                               equations    = E,
                               constraints  = G)

  X_init = [ (order+1)*[0] for i in range(N) ]
  initial_valuation = continuous.Valuation(independent = t_begin,
                                           dependent   = X_init)

  config = continuous.Config(tolerance = continuous.Tolerance(relative = 1.0e-4,
                                                              absolute = 1.0e-3),
                             step = continuous.Step(mode  = continuous.StepMode.fixed,
                                                    delta = t_delta),
                             stop = t_end)

  return problem,initial_valuation,config

def simulate(A,b):
  from copy import deepcopy

  continuous = hysj.simulators.continuous

  problem,initial_valuation,config = make_problem(A,b)

  solution = continuous.make_solution(problem           = problem,
                                      initial_valuation = initial_valuation,
                                      config            = config)

  states = []
  print('simulating...')
  while event := solution():
    state = solution.state()
    if event in [continuous.Event.init,
                 continuous.Event.step,
                 continuous.Event.stop]:
      states.append(deepcopy(state))
    if event == continuous.Event.stop:
      break
    if event == continuous.Event.fail:
      raise Exception('simulation failed')
  print('done...')
  dfa = solution.constraints().graphs[0]
  return states,dfa

def make_dfa_pos(dfa,xy,wh):
  import networkx

  box_f = 3 / 4
  wh = [ e * box_f for e in wh ]
  xy = [ xy[i] + wh[i] * (1 - box_f)/2 for i in range(len(xy)) ]

  nx = networkx.DiGraph()
  for v in dfa.vertices():
    nx.add_node(v.value)

  nx_edges = []
  for e in dfa.edges():
    v_in,v_out = dfa(e).ports
    nx.add_edge(v_in.value,v_out.value)
    nx_edges.append((v_in.value,v_out.value))

  nx_pos = list(networkx.drawing.nx_agraph.graphviz_layout(nx).values())
  nx_xmin = min([x for x,y in nx_pos])
  nx_xmax = max([x for x,y in nx_pos])
  nx_ymin = min([y for x,y in nx_pos])
  nx_ymax = max([y for x,y in nx_pos])
  nx_xext = nx_xmax - nx_xmin
  nx_yext = nx_ymax - nx_ymin

  def mapx(x):
    return ((x - nx_xmin) / nx_xext) * wh[0] + xy[0]
  def mapy(y):
    return ((y - nx_ymin) / nx_yext) * wh[1] + xy[1]

  return [ (mapx(x),mapy(y)) for x,y in nx_pos ]


dfa_vertex_radius = 0.04

def make_dfa_vertex_patches(dfa,pos):
  from matplotlib.patches import Circle
  return [ Circle([x,y],
                  radius = dfa_vertex_radius,
                  fill = True,
                  linewidth = 1.0,
                  facecolor = 'white',
                  edgecolor = 'black') for x,y in pos ]

def make_dfa_edge_patches(dfa,pos):
  from matplotlib.patches import FancyArrowPatch
  patches = []
  for e in dfa.edges():
    v_in,v_out = dfa(e).ports
    if v_in == v_out:
      continue

    r = dfa_vertex_radius / 2.0
    a = numpy.array(pos[v_in.value])
    b = numpy.array(pos[v_out.value])

    c = b - a
    n = c / numpy.linalg.norm(c)
    b = b - n * (1.1*r)
    patch = FancyArrowPatch(posA = a,
                            posB = b,
                            linewidth = 1.0,
                            arrowstyle='-|>,head_length=4.0,head_width=3.0',
                            edgecolor = 'black',
                            linestyle = (0,(1,3)),
                            facecolor = 'black',
                            connectionstyle='arc3,rad=0.5')

    patches.append(patch)
  return patches

def make_dfa_patches(dfa,xy,wh):
  pos = make_dfa_pos(dfa,xy,wh)
  return make_dfa_vertex_patches(dfa,pos),make_dfa_edge_patches(dfa,pos)

def animate(A,b,states,dfa):
  print('animating...')
  from matplotlib.animation import FuncAnimation as Animation
  from matplotlib.animation import writers as animation_writers
  N = len(A)

  xdata = [ [s.valuation.dependent[i][1] for s in states ] for i in range(N) ]
  ydata = [ [s.valuation.dependent[i][0] for s in states ] for i in range(N) ]
  vdata = [ s.constraints[0] if s.constraints else None for s in states ]


  figure,axes = plot.subplots(1,1)

  axes.set_aspect("equal", adjustable="datalim")
  axes.set_xlabel('$x_i$')
  axes.set_ylabel('$\dot{x_i}$')


  xmin,xmax = min([x for X in xdata for x in X]),max([x for X in xdata for x in X])
  ymin,ymax = min([y for Y in ydata for y in Y]),max([y for Y in ydata for y in Y])
  xext = xmax - xmin
  yext = ymax - ymin

  axes_margin = 0.1

  axes.set_xlim([xmin - axes_margin,xmax + axes_margin])
  axes.set_ylim([ymin - axes_margin,ymax + axes_margin])

  dfa_vertex_patches,dfa_edge_patches = make_dfa_patches(dfa,[xmin,ymin],[xext,yext])
  for p in dfa_edge_patches + dfa_vertex_patches:
    axes.add_patch(p)

  lines  = [ axes.plot([],[],label=f'$x_{{{i}}}$')[0] for i in range(N) ]
  points = [ axes.plot([],[],color=lines[i].get_color(),marker='o')[0] for i in range(N) ]

  axes.axhline(y = b,color='red',label='b')

  axes.legend(bbox_to_anchor=(0,1.02,1,0.2), loc="lower center",
              borderaxespad=0, ncol=(N + 1))

  frame_step   = 2
  frame_count  = len(states) // frame_step
  frame_window = frame_step * 10

  def draw_X(first,last):
    for i in range(N):
      lines[i].set_data(xdata[i][first:last],ydata[i][first:last])
      points[i].set_data(xdata[i][last:last+1],ydata[i][last:last+1])

  def draw_dfa(last):

    def vertex(i):
      if i < 0 or i >= len(states):
        return None
      return vdata[i]

    def leave(v):
      dfa_vertex_patches[v.value].set_facecolor('white')

    def enter(v):
      dfa_vertex_patches[v.value].set_facecolor('green' if dfa(v)() else 'red')

    v_prev = vertex(last - frame_step)
    v_next = vertex(last)
    if v_prev:
      leave(v_prev)
    enter(v_next)

  def animate(frame_index):
    last  = min((frame_index + 1) * frame_step,len(states) - 1)
    first = max(0,last - frame_window * frame_step)

    draw_X(first,last)
    draw_dfa(last)
    return lines + points

  animation = Animation(figure,
                        animate,
                        frames = frame_count,
                        interval = 1,
                        blit = True)

  animation.save(f'structured-transition-guards-{N}.mp4',
                 writer = animation_writers['ffmpeg'](fps = 30),
                 dpi = 100.0)
  print('done')

def make_animation(b,A):
  states,dfa = simulate(A,b)
  animate(A,b,states,dfa)

try:
  make_animation(A = [1.0          ],b = 0.5)
  make_animation(A = [0.75,1.25    ],b = 0.5)
  make_animation(A = [0.75,1.0,1.25],b = 0.5)
except Exception as e:
  print(str(e))