Remodeling a Flower System

John Eivind Helset

2021-02-24 03:55

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 HA 2:

\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: See Structured Transition Guards.

Simulating

The HA is simulated in python with hysj and IDA, using the original model parameters4 5. 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].

Glossary

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[0] + 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[0]]) 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[0][0],X[1][0]),
          ir.strict_inequality(ir.negation(X[0][0]),X[1][0]) ]

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

  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[0]],
      make_constraint = lambda m,l,s: region[other_location(l)[0]]))
  #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[0],
    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][0]) for t in trajectories for s in t[1]]) for i in variables ]
  margin = 0.1

  axes.set_xlim((-limit[0]-margin,limit[0]+margin))
  axes.set_ylim((-limit[1]-margin,limit[1]+margin))

  axes.grid(True)
  axes.set_xlabel(f'${latex_continuous_variables[0][0]}$')
  axes.set_ylabel(f'${latex_continuous_variables[1][0]}$')

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

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

  def legend_labels():
    R = [ f'$ray_{r}$' for r,_ in enumerate(problem.rays) ]
    L = [ f'$location_{l[0]}$' 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",
              borderaxespad=0, ncol=len(problem.rays) + len(locations))

  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[1][0].index()] = 0.0
  for i in range(len(problem.rays)):
    xdata = []
    ydata = []
    for x in grid[0]:
      calculator[problem.continuous.variables.dependent[0][0].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[0])):
        quiver_data[-1][0].append([])
        quiver_data[-1][1].append([])
        for b in range(len(quiver_grid[0][a])):
          for k in variables:
            calculator[problem.continuous.variables.dependent[k][0].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[0],
                         quiver_grid[1],
                         quiver_data[0][0],
                         quiver_data[0][1],
                         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][1]))] 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][0][0],[]).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][0],quiver_data[current_location][1])

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

      xdata = [ s.dependent[0][0] for k,s in S ]
      ydata = [ s.dependent[1][0] for k,s in S ]
      if S[0][0] == 0:
        artists.append(axes.plot(xdata[0],ydata[0],marker='.',color=color)[0])
      artists.append(axes.plot(xdata,ydata,color=color)[0])
    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)}')