# Remodeling a Flower System

Here we reformulate the piecewise affine flower system referred to in [Hedlund99] as a Lazy Hybrid Automaton and simulate it as such with hysj.

## Remodeling

The PWA flower system is restated here1 as:

\begin{equation*} \dot{x}(t) = \begin{cases} A_0x(t), \quad x^2_0(t) - x^2_1(t) \geq 0 \\ A_1x(t), \quad x^2_0(t) - x^2_1(t) < 0\end{cases} \end{equation*}
\begin{equation*} A_0 = \begin{bmatrix} \epsilon & \alpha \omega \\ -\omega & -\epsilon \end{bmatrix} \end{equation*}
\begin{equation*} A_1 = \begin{bmatrix} -\epsilon & \omega \\ -\alpha\omega & -\epsilon \end{bmatrix} \end{equation*}

Which is remodeled as the HA2:

\begin{equation*} F_i = A_i x \end{equation*}
\begin{equation*} G_{ii^\prime} = x^2_i - x^2_{i^\prime} \end{equation*}
\begin{equation*} E_i = \{ l_{i^\prime} \} \end{equation*}

The constraints on the guards of the HA are here only that they are expressible as logical expressions over functions whose roots can be detected by the solver used for simulation3. The difference in squares is here used directly, but one could also f.ex. have formulated the guard $$x^2_{i^\prime} - x^2_i < 0$$ as $$(x_i > x_{i^\prime}) \land (-x_i < x_{i^\prime})$$.

1

Using null-indexing.

2

Using subscripts $$F_i$$ instead of $$F(l_i)$$.

3

## Simulating

The HA is simulated in python with hysj and IDA, using the original model parameters45. The simulation is illustrated with matplotlib and numpy. The rays used to define the regions of the state-space in [Hedlund99] are added for illustrative purposes. Roots are drawn with dots. The phase plane of the currently active location is drawn in the background. Code listing below.

4

Ignoring time management.

5

The bootstrapping process outlined in Structured Transition Guards ensures that one does not need to specify the correct initial discrete location since the HA is deterministic in the sense of Kowalewski et al. in [Lunze09].

PWA

Piecewise affine

HA

Hybrid automaton

## Bibliography

Hedlund99(1,2)

"Sven Hedlund", "Computational Methods for Hybrid Systems"

Hindmarsh2005

"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"

Lunze09

"Jan Lunze", "Handbook of Hybrid Systems Control"

## Code

import hysj

continuous = hysj.simulators.continuous
discrete   = hysj.simulators.discrete
hybrid     = hysj.simulators.hybrid

regions   = range(2)
variables = range(2)
order     = 1

locations = [(0,),(1,)]
def other_location(l):
return ((l + 1)%2,)

symbol = 0

def make_problem():

ir = hysj.mathematics.Program()

alpha,omega,epsilon = list(map(ir.constant,[5.0,1.0,0.1]))
zero = ir.zero()

A = [ [ [ ir.negation(epsilon), ir.multiplication([alpha,omega]) ],
[ ir.negation(omega),   ir.negation(epsilon)             ] ],
[ [ ir.negation(epsilon),                           omega                ],
[ ir.multiplication([ir.negation(alpha), omega]), ir.negation(epsilon) ] ] ]

t = ir.symbol()
X = [ ir.variable(independent_variable = t,order = 1,info = i) for i in variables ]

def f(region,variable):
return ir.addition([ir.multiplication([a,x]) for a,x in zip(A[region][variable],X) ])

F = [ [ f(i,j) for j in variables ] for i in regions ]

equations = [ [ ir.equality(X[j][-1],F[i][j]) for j in variables ] for i in regions ]

ray = [ ir.strict_inequality(X,X),
ir.strict_inequality(ir.negation(X),X) ]

abs_greater = [ ir.assertion(ir.strict_inequality(ir.power(X[i],ir.constant(2)),
ir.power(X[(i+1)%2],ir.constant(2))))
for i in variables ]
region = [ abs_greater,
ir.complement(abs_greater) ]

problem = hysj.simulators.hybrid.Problem(
program = ir,
discrete = hysj.simulators.discrete.Problem(
make_symbols = lambda l: [symbol],
make_target  = lambda l,s: other_location(l),
is_immediate = lambda l,s: False),
continuous = hysj.simulators.hybrid.ContinuousProblem(
variables       = hysj.simulators.continuous.Variables(independent = t,dependent = X),
make_equations  = lambda m,l:   equations[l],
make_constraint = lambda m,l,s: region[other_location(l)]))
#NOTE: for convenience - jeh
problem.activity_function = F
problem.rays              = ray
return problem

#NOTE: fixed-step simulation of problem for t e [0,20.0]
def simulate(problem,initial_continuous_valuation):
from copy import deepcopy

t_begin,t_end,t_count = (0.0,20.0,2000+1)
t_delta  = (t_end - t_begin) / (t_count - 1)

initial_valuation = hybrid.Valuation(
discrete   = locations,
continuous = continuous.Valuation(
independent = t_begin,
dependent  = initial_continuous_valuation))

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

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

trajectories = []
print('simulating...')
while event := solution():
if event == discrete.Event.vertex:
trajectories.append((tuple(solution.discrete().state().location()),[]))
if event == continuous.Event.step or \
event == continuous.Event.init or \
event == continuous.Event.root:
trajectories[-1][-1].append(deepcopy(solution.continuous().state().valuation))
if event == continuous.Event.fail:
raise Exception('simulation failed')
print('done...')
return trajectories

def animate(problem,trajectories):

import numpy
import matplotlib.pyplot as plot
from matplotlib import lines
from matplotlib.animation import FuncAnimation as Animation
from matplotlib.animation import writers as animation_writers
import hysj.latex

print('animating...')

def symbol_labeler(s,d):
return {tuple(problem.continuous.variables.independent.index()): 't'}.get(tuple(s.index()),f'x_{d}')

latex_renderer = hysj.latex.Renderer(problem.program,symbol_labeler)
latex_continuous_variables = [ list(map(latex_renderer,X)) for X in problem.continuous.variables.dependent ]

figure = plot.figure()
axes   = plot.axes()

legend_lines  = []
legend_labels = []

limit  = [ max([abs(s.dependent[i]) for t in trajectories for s in t]) for i in variables ]
margin = 0.1

axes.set_xlim((-limit-margin,limit+margin))
axes.set_ylim((-limit-margin,limit+margin))

axes.grid(True)
axes.set_xlabel(f'${latex_continuous_variables}$')
axes.set_ylabel(f'${latex_continuous_variables}$')

def colors(i):
return ['#377eb8', '#ee6e00', '#4daf4a',
'#f781bf', '#a65628', '#984ea3',
'#999999', '#e41a1c', '#dede00'][i]

def legend_lines():
R = [ lines.Line2D(,,color = colors(len(regions) + r),linestyle = '--')
for r,_ in enumerate(problem.rays) ]
L = [ lines.Line2D(,,color = colors(l)) for l in locations ]
return R + L

def legend_labels():
R = [ f'$ray_{r}$' for r,_ in enumerate(problem.rays) ]
L = [ f'$location_{l}$' for l in locations ]
return R + L

axes.legend(legend_lines(),
legend_labels(),
bbox_to_anchor=(0,1.02,1,0.2), loc="lower center",

grid       = [ numpy.linspace(-limit[i],limit[i],num = 20) for i in variables ]
calculator = hysj.calculators.Basic(problem.program)

#NOTE: compute the rays - jeh
calculator[problem.continuous.variables.dependent.index()] = 0.0
for i in range(len(problem.rays)):
xdata = []
ydata = []
for x in grid:
calculator[problem.continuous.variables.dependent.index()] = x
xdata.append(x)
ydata.append(calculator(problem.rays[i]))
plot.plot(xdata,ydata,color = colors(len(regions) + i),linestyle = '--')

#NOTE: compute quiver of phase plane - jeh
def make_quiver():
quiver_grid = numpy.meshgrid(*grid)
quiver_data = []
for i in regions:
quiver_data.append(([],[]))
for a in range(len(quiver_grid)):
quiver_data[-1].append([])
quiver_data[-1].append([])
for b in range(len(quiver_grid[a])):
for k in variables:
calculator[problem.continuous.variables.dependent[k].index()] = quiver_grid[k][a][b]
for k in variables:
quiver_data[-1][k][-1].append(calculator(problem.activity_function[i][k]))

quiver = plot.quiver(quiver_grid,
quiver_grid,
quiver_data,
quiver_data,
pivot='mid',
color=colors(0),width=0.003)
return quiver_data,quiver

quiver_data,quiver = make_quiver()

#NOTE: set up animation - jeh
state_coords = sum([[(i,j) for j in range(len(trajectories[i]))] for i in range(len(trajectories))],[])
frame_step   = 10
frame_count  = len(state_coords) // frame_step

def render(frame_index):
coords = {}
for i in range(max(frame_step * (frame_index - 1),0),min(frame_step*frame_index + 1,len(state_coords))):
j,k = state_coords[i]
coords.setdefault(trajectories[j],[]).append((k,trajectories[j][-1][k]))

if not coords:
return []

artists = []

current_location = max(coords.keys())

quiver.set_color(colors(current_location))
quiver.set_UVC(quiver_data[current_location],quiver_data[current_location])

for l,S in coords.items():
color = colors(l)

xdata = [ s.dependent for k,s in S ]
ydata = [ s.dependent for k,s in S ]
if S == 0:
artists.append(axes.plot(xdata,ydata,marker='.',color=color))
artists.append(axes.plot(xdata,ydata,color=color))
return artists

Animation(figure,render,frames = frame_count,interval = 1,blit = True)\
.save(
'remodeling-a-flower-system.mp4',
writer = animation_writers['ffmpeg'](fps = 30),
dpi = 100.0)
print('done...')

try:
problem      = make_problem()
trajectories = simulate(
problem = problem,
initial_continuous_valuation = [[1.0,0.0],[0.0,0.0]])
animate(problem,trajectories)
except Exception as e:
print(f'exception: {str(e)}')