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]:
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:
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]:
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:
The states of the DFA, \(Q\), correspond to the elements of \(\mathbb{B}^N\):
\(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\):
The terminals are the states for whose valuations the guard evaluates to true.
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:
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\):
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) \}\):
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:
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)\):
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:
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|}}\):
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))