Examples for CDC '25¶
Note: The results presented in the paper was generated by a WIP-version of pyspect, depending on a custom HZ library (not ZonoOpt). In favor of ZonoOpt, a significantly better implementation of zonotopes, I've updated this notebook to use it instead. The original notebook can be found in the branch
cdc25.
Case Study: Overtaking Scenario¶
The running example from the paper, illustrating basic concepts in pyspect. The task is to verify an overtaking scenario.
Imports¶
In [1]:
Copied!
from math import pi
from time import time
from functools import reduce
from contextlib import contextmanager
from pyspect import *
from math import pi
from time import time
from functools import reduce
from contextlib import contextmanager
from pyspect import *
In [2]:
Copied!
@contextmanager
def timectx(msgfunc):
"""Context manager to time a block of code."""
start = time()
yield
end = time()
print(msgfunc(end-start))
@contextmanager
def timectx(msgfunc):
"""Context manager to time a block of code."""
start = time()
yield
end = time()
print(msgfunc(end-start))
Hyperparameters¶
In [3]:
Copied!
AXES = [dict(name= 't', bounds=[ 0, 15], step=0.5, unit='s'),
dict(name= 'x', bounds=[ 0, 200], points=201, unit='m'),
dict(name= 'y', bounds=[ -4, +4], points= 15, unit='m'),
dict(name='yaw', bounds=[-pi/16, +pi/16], points= 13, unit='rad'), # pi/16 = 11.25deg
dict(name='vel', bounds=[ 15, 30], points= 51, unit='m/s')] # 15 mps ~= 54 kmph | 35 mps ~= 125 kmph
MAX_STEER = pi/64 # [rad/s] +/- 2.8 degrees/s
MAX_ACCEL = 3.5 # [mps2]
AXES = [dict(name= 't', bounds=[ 0, 15], step=0.5, unit='s'),
dict(name= 'x', bounds=[ 0, 200], points=201, unit='m'),
dict(name= 'y', bounds=[ -4, +4], points= 15, unit='m'),
dict(name='yaw', bounds=[-pi/16, +pi/16], points= 13, unit='rad'), # pi/16 = 11.25deg
dict(name='vel', bounds=[ 15, 30], points= 51, unit='m/s')] # 15 mps ~= 54 kmph | 35 mps ~= 125 kmph
MAX_STEER = pi/64 # [rad/s] +/- 2.8 degrees/s
MAX_ACCEL = 3.5 # [mps2]
Program¶
Task definition¶
Define the different regions. These are "constants" in the task specification. Note that we can define the sets lazily, such that we do not rely on a specific implementation or set representation yet. This is possible because
pyspectinternally defers evaluation until realization time.Write the specification $\varphi = \varphi_\text{env} \:\mathsf{U}\: \texttt{goal}$, where $\varphi_\text{env}, \texttt{goal} \in \mathsf{AP}$. While $\varphi_\text{env} = ( \neg p_\text{obs} \wedge ... )$ with $p_\text{obs} \in \mathsf{AP}$, and $p_\text{obs} \leftrightarrow z \in \texttt{OBS}$, we directly use the set $\texttt{OBS}$ in the specification. The proposition $\texttt{goal}$ is left symbolic and will be bound later.
Note A. This example uses functools.reduce to fold-left binary operators (e.g. AND). This is simply to reduce the code length. In this example, the semantics allow for both left- and right-associative reduction (which is not necessarily true for all downstream applications). By default, pyspect does not assume associativity for any of its operators.
Note B. We define OBS in a funny way to make it easier to write GOAL. We want to describe the goal as "in front of OBS-vehicle", hence we need to know the front of OBS. When having this, it becomes shorter to (informally) define OBS such that BEHIND_OBS < OBS < INFRONT_OBS.
In [4]:
Copied!
## (1) Define the different regions.
# All of the following objects are SetBuilder objects and will be propositions in later specs.
RIGHT_LANE = BoundedSet(y=(-3.5, -0.5))
LEFT_LANE = BoundedSet(y=(+0.5, +3.5))
BETWEEN = BoundedSet(y=(-1.5, +1.5))
STRAIGHT = BoundedSet(yaw=(-pi/18, +pi/18))
INFRONT_OBS = HalfSpaceSet(normal=[+1, -16],
offset=[40, 4],
axes=['x', 't'])
BEHIND_OBS = HalfSpaceSet(normal=[-1, +16],
offset=[0, 4],
axes=['x', 't'])
# See Note B.
OBS = reduce(Inter, [Compl(INFRONT_OBS), Compl(BEHIND_OBS), RIGHT_LANE])
GOAL = reduce(Inter, [
INFRONT_OBS,
RIGHT_LANE,
STRAIGHT,
# NOTE: For stability of HJ solution,
# target should not extend to domain boundary
BoundedSet(x=(..., 195)),
])
## (2) Write the specification.
# All of the following objects are logical specifications.
ENV = reduce(AND, [
reduce(OR, [
AND(RIGHT_LANE, STRAIGHT),
AND(LEFT_LANE, STRAIGHT),
BETWEEN,
]),
NOT(OBS),
# NOTE: For stability of HJ solution,
# target should not extend to domain boundary
BoundedSet(vel=(16, 24), yaw=(-pi/16 + pi/32, +pi/16 - pi/32)),
])
TASK = UNTIL(ENV, 'goal')
## (1) Define the different regions.
# All of the following objects are SetBuilder objects and will be propositions in later specs.
RIGHT_LANE = BoundedSet(y=(-3.5, -0.5))
LEFT_LANE = BoundedSet(y=(+0.5, +3.5))
BETWEEN = BoundedSet(y=(-1.5, +1.5))
STRAIGHT = BoundedSet(yaw=(-pi/18, +pi/18))
INFRONT_OBS = HalfSpaceSet(normal=[+1, -16],
offset=[40, 4],
axes=['x', 't'])
BEHIND_OBS = HalfSpaceSet(normal=[-1, +16],
offset=[0, 4],
axes=['x', 't'])
# See Note B.
OBS = reduce(Inter, [Compl(INFRONT_OBS), Compl(BEHIND_OBS), RIGHT_LANE])
GOAL = reduce(Inter, [
INFRONT_OBS,
RIGHT_LANE,
STRAIGHT,
# NOTE: For stability of HJ solution,
# target should not extend to domain boundary
BoundedSet(x=(..., 195)),
])
## (2) Write the specification.
# All of the following objects are logical specifications.
ENV = reduce(AND, [
reduce(OR, [
AND(RIGHT_LANE, STRAIGHT),
AND(LEFT_LANE, STRAIGHT),
BETWEEN,
]),
NOT(OBS),
# NOTE: For stability of HJ solution,
# target should not extend to domain boundary
BoundedSet(vel=(16, 24), yaw=(-pi/16 + pi/32, +pi/16 - pi/32)),
])
TASK = UNTIL(ENV, 'goal')
/tmp/ipykernel_2551/710006715.py:5: DeprecationWarning: BoundedSet is renamed to AlignedBoxSet for clarity, please use AlignedBoxSet instead of BoundedSet. RIGHT_LANE = BoundedSet(y=(-3.5, -0.5)) /tmp/ipykernel_2551/710006715.py:6: DeprecationWarning: BoundedSet is renamed to AlignedBoxSet for clarity, please use AlignedBoxSet instead of BoundedSet. LEFT_LANE = BoundedSet(y=(+0.5, +3.5)) /tmp/ipykernel_2551/710006715.py:7: DeprecationWarning: BoundedSet is renamed to AlignedBoxSet for clarity, please use AlignedBoxSet instead of BoundedSet. BETWEEN = BoundedSet(y=(-1.5, +1.5)) /tmp/ipykernel_2551/710006715.py:9: DeprecationWarning: BoundedSet is renamed to AlignedBoxSet for clarity, please use AlignedBoxSet instead of BoundedSet. STRAIGHT = BoundedSet(yaw=(-pi/18, +pi/18)) /tmp/ipykernel_2551/710006715.py:28: DeprecationWarning: BoundedSet is renamed to AlignedBoxSet for clarity, please use AlignedBoxSet instead of BoundedSet. BoundedSet(x=(..., 195)), /tmp/ipykernel_2551/710006715.py:44: DeprecationWarning: BoundedSet is renamed to AlignedBoxSet for clarity, please use AlignedBoxSet instead of BoundedSet. BoundedSet(vel=(16, 24), yaw=(-pi/16 + pi/32, +pi/16 - pi/32)),
Construct the TLT¶
Select the set of primitive TLTs. This determines the logic fragment in which the specification is interpreted. In this example, we use continuous-time LTL, i.e., LTL without the "next" operator.
Construct the temporal logic tree and bind the symbolic proposition $\texttt{goal}$ to the concrete set $\texttt{GOAL}$. The call to
.where(goal=GOAL)updates a internal proposition map $\mathsf{M}$ so that $\texttt{goal} \leftrightarrow z \in \texttt{GOAL}$.
In [5]:
Copied!
## (3) Select the set of primitive TLTs.
TLT.select(ContLTL)
## (4) Create the TLT and set the proposition 'goal'.
objective = TLT(TASK, where={'goal': GOAL})
## (3) Select the set of primitive TLTs.
TLT.select(ContLTL)
## (4) Create the TLT and set the proposition 'goal'.
objective = TLT(TASK, where={'goal': GOAL})
Initialize the implementation object¶
First now, we introduce a specific implementation. pyspect has some off-the-shelf implementations that can be imported. Note, however, pyspect does not abstract away the underlying back-end completely, as we see in the following cell.
Additional imports for the specific implementation.
TVHJImplis a wrapper class around thehj_reachabilitylibrary.pyspectincludes a forked version with convenience utilities.Initialization of the implementation is fully backend-specific. In this case, we specify the dynamics (double integrator), bounds, time horizon, and discretization. This configuration determines how sets like
GOALand operators likeUNTILwill be computed (e.g., using backward reachable sets on a grid).Realize the TLT. This is the key step where the temporal logic and reachability analysis meet:
objective.realize(impl)recursively applies the reachability and set operations defined by the selected logic fragment and implemented in the backend. The outputoutis the satisfaction set (i.e., the states from which the specification holds) in the implementation's set reprentation.Plotting. We visualize the satisfaction set using
HJImpl.plot, which displays the reachable states as a level set over position, velocity, and time. This step is backend-dependent (as value functions are used).
In [6]:
Copied!
## (5) Additional imports for the specific implementation.
from pyspect.impls.hj_reachability import TVHJImpl
from pyspect.systems.hj_reachability import *
## (6) Define the implementation of the reachability algorithm.
dynamics = dict(cls=Bicycle4D,
wheelbase=2.7,
min_accel=-MAX_ACCEL,
max_accel=+MAX_ACCEL,
min_steer=-MAX_STEER,
max_steer=+MAX_STEER)
impl = TVHJImpl(dynamics, AXES)
## (7) Run the reachability program ##
with timectx(lambda t: f"Realization with HJ took {t:.2f} seconds"):
out = objective.realize(impl)
# `out` will have the same object type that `impl` operates with.
# For TVHJImpl, `out` will be a numpy array of the gridded value function.
print(f"{type(out) = }, {out.min() = }, {out.max() = }")
## (8) Plotting
# Generate levelset for the other vehicle.
obs = OBS(impl)
print(f"{obs.min() = }, {obs.max() = }")
# NOTE: We remove a small margin for plotting (clearer boundary/separation from obs).
out = out + 0.03
impl.plot(
out,
(obs, dict(colorscale='reds')),
axes=('x', 'y', 't'),
camera_eye=impl.PLOT.EYE_ML_SW,
camera_center=dict(x=.9, y=.9, z=-.8),
layout_margin=dict(t=40, r=0, l=0, b=0),
).show()
## (5) Additional imports for the specific implementation.
from pyspect.impls.hj_reachability import TVHJImpl
from pyspect.systems.hj_reachability import *
## (6) Define the implementation of the reachability algorithm.
dynamics = dict(cls=Bicycle4D,
wheelbase=2.7,
min_accel=-MAX_ACCEL,
max_accel=+MAX_ACCEL,
min_steer=-MAX_STEER,
max_steer=+MAX_STEER)
impl = TVHJImpl(dynamics, AXES)
## (7) Run the reachability program ##
with timectx(lambda t: f"Realization with HJ took {t:.2f} seconds"):
out = objective.realize(impl)
# `out` will have the same object type that `impl` operates with.
# For TVHJImpl, `out` will be a numpy array of the gridded value function.
print(f"{type(out) = }, {out.min() = }, {out.max() = }")
## (8) Plotting
# Generate levelset for the other vehicle.
obs = OBS(impl)
print(f"{obs.min() = }, {obs.max() = }")
# NOTE: We remove a small margin for plotting (clearer boundary/separation from obs).
out = out + 0.03
impl.plot(
out,
(obs, dict(colorscale='reds')),
axes=('x', 'y', 't'),
camera_eye=impl.PLOT.EYE_ML_SW,
camera_center=dict(x=.9, y=.9, z=-.8),
layout_margin=dict(t=40, r=0, l=0, b=0),
).show()
--------------------------------------------------------------------------- ModuleNotFoundError Traceback (most recent call last) Cell In[6], line 3 1 ## (5) Additional imports for the specific implementation. 2 ----> 3 from pyspect.impls.hj_reachability import TVHJImpl 4 5 from pyspect.systems.hj_reachability import * 6 File ~/work/pyspect/pyspect/src/pyspect/impls/hj_reachability.py:32 30 import jax 31 import jax.numpy as jnp ---> 32 from jax_tqdm import scan_tqdm 33 import hj_reachability as hj 34 from hj_reachability import utils ModuleNotFoundError: No module named 'jax_tqdm'
In [7]:
Copied!
for t in [6., 8., 10., 12.]:
impl.plot(
out,
(obs, dict(colorscale='reds')),
method='bitmap',
axes=('x', 'y'),
transform_select=[('t', t)],
layout_height=300,
layout_width=600,
layout_title=f'At Time: {t = }',
).show()
for t in [6., 8., 10., 12.]:
impl.plot(
out,
(obs, dict(colorscale='reds')),
method='bitmap',
axes=('x', 'y'),
transform_select=[('t', t)],
layout_height=300,
layout_width=600,
layout_title=f'At Time: {t = }',
).show()
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[7], line 2 1 for t in [6., 8., 10., 12.]: ----> 2 impl.plot( 3 out, 4 (obs, dict(colorscale='reds')), 5 method='bitmap', NameError: name 'impl' is not defined
Specification: $\square \psi$¶
This specification requires the system to remain within a safe region $\psi$ for the entire time horizon. It encodes an invariance property, meaning the condition $\psi$ must hold at all times. The HJ backend realizes this using an avoid set computed via universal control, solving $\Box \psi$ as a special case of $\neg \Diamond \neg \psi$, while the HZ backend uses a fixed-point formulation, i.e. $\psi \land \bigcirc(\psi \land \bigcirc(\dots))$. Both methods result a valid set of initial states from which the system can always satisfy $\psi$ throughout the trajectory.
Imports¶
In [8]:
Copied!
# Implementation independent
from time import time
from contextlib import contextmanager
from tqdm import tqdm
from pyspect import *
# HJ specific
from pyspect.impls.hj_reachability import TVHJImpl
from pyspect.systems.hj_reachability import DoubleIntegrator as HJDoubleIntegrator
from pyspect.plotting.levelsets import *
# HZ specific
import zonoopt as zono
from pyspect.impls.zonoopt import ZonoOptImpl
from pyspect.impls.zonoopt import DoubleIntegrator as HZDoubleIntegrator
# Implementation independent
from time import time
from contextlib import contextmanager
from tqdm import tqdm
from pyspect import *
# HJ specific
from pyspect.impls.hj_reachability import TVHJImpl
from pyspect.systems.hj_reachability import DoubleIntegrator as HJDoubleIntegrator
from pyspect.plotting.levelsets import *
# HZ specific
import zonoopt as zono
from pyspect.impls.zonoopt import ZonoOptImpl
from pyspect.impls.zonoopt import DoubleIntegrator as HZDoubleIntegrator
--------------------------------------------------------------------------- ModuleNotFoundError Traceback (most recent call last) Cell In[8], line 5 1 # Implementation independent 2 3 from time import time 4 from contextlib import contextmanager ----> 5 from tqdm import tqdm 6 7 from pyspect import * 8 ModuleNotFoundError: No module named 'tqdm'
In [9]:
Copied!
@contextmanager
def timectx(msgfunc):
"""Context manager to time a block of code."""
start = time()
yield
end = time()
print(msgfunc(end-start))
def print_hzinfo(hz):
print(f'n = {hz.get_n()}, nGc = {hz.get_nGc()}, nGb = {hz.get_nGb()}, nC = {hz.get_nC()}')
@contextmanager
def timectx(msgfunc):
"""Context manager to time a block of code."""
start = time()
yield
end = time()
print(msgfunc(end-start))
def print_hzinfo(hz):
print(f'n = {hz.get_n()}, nGc = {hz.get_nGc()}, nGb = {hz.get_nGb()}, nC = {hz.get_nC()}')
Hyperparameters¶
In [10]:
Copied!
AXES = [dict(name='t', bounds=[ 0, 40], step=0.5, unit='s'),
dict(name='x', bounds=[-100, +100], points= 91, unit='m'),
dict(name='v', bounds=[ -20, +20], points= 91, unit='m/s')]
MAX_ACCEL = 1.0 # [mps2]
# HZ2HJ specific
TIME_STEP = AXES[0]['step']
TIME_HORIZON = AXES[0]['bounds'][1]
MIN_BOUNDS = [spec['bounds'][0] for spec in AXES if spec['name'] != 't']
MAX_BOUNDS = [spec['bounds'][1] for spec in AXES if spec['name'] != 't']
GRID_SHAPE = [spec['points'] for spec in AXES if spec['name'] != 't']
AXES = [dict(name='t', bounds=[ 0, 40], step=0.5, unit='s'),
dict(name='x', bounds=[-100, +100], points= 91, unit='m'),
dict(name='v', bounds=[ -20, +20], points= 91, unit='m/s')]
MAX_ACCEL = 1.0 # [mps2]
# HZ2HJ specific
TIME_STEP = AXES[0]['step']
TIME_HORIZON = AXES[0]['bounds'][1]
MIN_BOUNDS = [spec['bounds'][0] for spec in AXES if spec['name'] != 't']
MAX_BOUNDS = [spec['bounds'][1] for spec in AXES if spec['name'] != 't']
GRID_SHAPE = [spec['points'] for spec in AXES if spec['name'] != 't']
Program¶
In [11]:
Copied!
## SPECIFICATION
T = BoundedSet(x=(-50, +50))
phi = ALWAYS(T)
## SPECIFICATION
T = BoundedSet(x=(-50, +50))
phi = ALWAYS(T)
/tmp/ipykernel_2551/2928891190.py:3: DeprecationWarning: BoundedSet is renamed to AlignedBoxSet for clarity, please use AlignedBoxSet instead of BoundedSet. T = BoundedSet(x=(-50, +50))
Two interpretations of ALWAYS¶
In [12]:
Copied!
# Define ALWAYS through RCI set (relating to ¬◇¬ψ)
@primitive(ALWAYS('_1'))
def Always_rci(_1: 'TLTLike') -> tuple[SetBuilder, APPROXDIR]:
b1, a1 = _1._builder, _1._approx
ao = APPROXDIR.UNDER
return (
AppliedSet('rci', b1),
ao + a1 if ao * a1 == APPROXDIR.EXACT else
a1 if ao == a1 else
APPROXDIR.INVALID,
)
# Define Always through fixed-point iteration
@primitive(ALWAYS('_1'))
def Always_fp(_1: 'TLTLike') -> tuple[SetBuilder, APPROXDIR]:
N = int(TIME_HORIZON / TIME_STEP)
phi = 'psi'
for _ in range(N-1):
phi = AND('psi', NEXT(phi))
tree = TLT(phi, where={'psi': _1})
return (tree._builder, tree._approx)
# Define ALWAYS through RCI set (relating to ¬◇¬ψ)
@primitive(ALWAYS('_1'))
def Always_rci(_1: 'TLTLike') -> tuple[SetBuilder, APPROXDIR]:
b1, a1 = _1._builder, _1._approx
ao = APPROXDIR.UNDER
return (
AppliedSet('rci', b1),
ao + a1 if ao * a1 == APPROXDIR.EXACT else
a1 if ao == a1 else
APPROXDIR.INVALID,
)
# Define Always through fixed-point iteration
@primitive(ALWAYS('_1'))
def Always_fp(_1: 'TLTLike') -> tuple[SetBuilder, APPROXDIR]:
N = int(TIME_HORIZON / TIME_STEP)
phi = 'psi'
for _ in range(N-1):
phi = AND('psi', NEXT(phi))
tree = TLT(phi, where={'psi': _1})
return (tree._builder, tree._approx)
HJ Implementation¶
In [13]:
Copied!
## CONSTRUCT TLT
TLT.select(ContLTL | Always_rci)
tree = TLT(phi)
print(f"Approximation direction: {tree._approx = }")
## INITIALIZE IMPLEMENTATION
dynamics = dict(cls=HJDoubleIntegrator,
min_accel=-MAX_ACCEL,
max_accel=+MAX_ACCEL)
impl = TVHJImpl(dynamics, AXES)
## REALIZE ROOT NODE OF TLT
with timectx(lambda t: f"Realization with HJ took {t:.2f} seconds"):
out = tree.realize(impl)
print(f"{type(out) = }")
## PLOT
impl.plot(
out,
colorscale='greens',
## ## ##
method='2D',
axes=('x', 'v'),
transform_select=[('t', 0)],
# ## ## ##
# method='3D',
# axes=('x', 'v', 't'),
# camera_eye=impl.PLOT.EYE_MH_W,
)
## CONSTRUCT TLT
TLT.select(ContLTL | Always_rci)
tree = TLT(phi)
print(f"Approximation direction: {tree._approx = }")
## INITIALIZE IMPLEMENTATION
dynamics = dict(cls=HJDoubleIntegrator,
min_accel=-MAX_ACCEL,
max_accel=+MAX_ACCEL)
impl = TVHJImpl(dynamics, AXES)
## REALIZE ROOT NODE OF TLT
with timectx(lambda t: f"Realization with HJ took {t:.2f} seconds"):
out = tree.realize(impl)
print(f"{type(out) = }")
## PLOT
impl.plot(
out,
colorscale='greens',
## ## ##
method='2D',
axes=('x', 'v'),
transform_select=[('t', 0)],
# ## ## ##
# method='3D',
# axes=('x', 'v', 't'),
# camera_eye=impl.PLOT.EYE_MH_W,
)
Approximation direction: tree._approx = <APPROXDIR.UNDER: -1>
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[13], line 11 7 print(f"Approximation direction: {tree._approx = }") 8 9 ## INITIALIZE IMPLEMENTATION 10 ---> 11 dynamics = dict(cls=HJDoubleIntegrator, 12 min_accel=-MAX_ACCEL, 13 max_accel=+MAX_ACCEL) 14 NameError: name 'HJDoubleIntegrator' is not defined
HZ Implementation¶
In [14]:
Copied!
## CONSTRUCT TLT
TLT.select(DiscLTL | Always_fp)
tree = TLT(phi)
print(f"Approximation direction: {tree._approx = }")
## INITIALIZE IMPLEMENTATION
dynamics = HZDoubleIntegrator(max_accel=MAX_ACCEL, dt=TIME_STEP)
impl = ZonoOptImpl(dynamics,
AXES[1:], # exclude time axis
zono.interval_2_zono(zono.Box([-MAX_ACCEL], [+MAX_ACCEL])),
TIME_HORIZON,
time_step=TIME_STEP)
## REALIZE ROOT NODE OF TLT
with timectx(lambda t: f"Realization with HZ took {t:.2f} seconds"):
out = tree.realize(impl)
print(f"{type(out) = }")
print_hzinfo(out)
## Plot
impl.plot_fill(
out,
axes=('x', 'v'),
)
## CONSTRUCT TLT
TLT.select(DiscLTL | Always_fp)
tree = TLT(phi)
print(f"Approximation direction: {tree._approx = }")
## INITIALIZE IMPLEMENTATION
dynamics = HZDoubleIntegrator(max_accel=MAX_ACCEL, dt=TIME_STEP)
impl = ZonoOptImpl(dynamics,
AXES[1:], # exclude time axis
zono.interval_2_zono(zono.Box([-MAX_ACCEL], [+MAX_ACCEL])),
TIME_HORIZON,
time_step=TIME_STEP)
## REALIZE ROOT NODE OF TLT
with timectx(lambda t: f"Realization with HZ took {t:.2f} seconds"):
out = tree.realize(impl)
print(f"{type(out) = }")
print_hzinfo(out)
## Plot
impl.plot_fill(
out,
axes=('x', 'v'),
)
Approximation direction: tree._approx = <APPROXDIR.OVER: 1>
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[14], line 11 7 print(f"Approximation direction: {tree._approx = }") 8 9 ## INITIALIZE IMPLEMENTATION 10 ---> 11 dynamics = HZDoubleIntegrator(max_accel=MAX_ACCEL, dt=TIME_STEP) 12 impl = ZonoOptImpl(dynamics, 13 AXES[1:], # exclude time axis 14 zono.interval_2_zono(zono.Box([-MAX_ACCEL], [+MAX_ACCEL])), NameError: name 'HZDoubleIntegrator' is not defined
Specification: $\square \lozenge \psi$¶
This specification ensures that the system can reach $\psi$ infinitely often, i.e., $\psi$ must remain recurrently reachable throughout the time horizon. Formally, it is a liveness property requiring that at every point in time, it is possible to reach a state satisfying $\psi$ in the future. The HJ backend supports this through nested reachability using negation and disjunction (e.g., $\Box \Diamond \psi = \neg \Diamond \neg \Diamond \psi$). However, HZ implementations typically cannot soundly support this due to approximation mismatches may violate the logical semantics when doing the fixed-point iteration. In response, pyspect flags and rejects this combination.
Imports¶
In [15]:
Copied!
# Implementation independent
from time import time
from contextlib import contextmanager
from tqdm import tqdm
from pyspect import *
# HJ specific
from pyspect.impls.hj_reachability import TVHJImpl
from pyspect.systems.hj_reachability import DoubleIntegrator as HJDoubleIntegrator
from pyspect.plotting.levelsets import *
# HZ specific
import zonoopt as zono
from pyspect.impls.zonoopt import ZonoOptImpl
from pyspect.impls.zonoopt import DoubleIntegrator as HZDoubleIntegrator
# Implementation independent
from time import time
from contextlib import contextmanager
from tqdm import tqdm
from pyspect import *
# HJ specific
from pyspect.impls.hj_reachability import TVHJImpl
from pyspect.systems.hj_reachability import DoubleIntegrator as HJDoubleIntegrator
from pyspect.plotting.levelsets import *
# HZ specific
import zonoopt as zono
from pyspect.impls.zonoopt import ZonoOptImpl
from pyspect.impls.zonoopt import DoubleIntegrator as HZDoubleIntegrator
--------------------------------------------------------------------------- ModuleNotFoundError Traceback (most recent call last) Cell In[15], line 5 1 # Implementation independent 2 3 from time import time 4 from contextlib import contextmanager ----> 5 from tqdm import tqdm 6 7 from pyspect import * 8 ModuleNotFoundError: No module named 'tqdm'
In [16]:
Copied!
@contextmanager
def timectx(msgfunc):
"""Context manager to time a block of code."""
start = time()
yield
end = time()
print(msgfunc(end-start))
def print_hzinfo(hz):
print(f'n = {hz.get_n()}, nGc = {hz.get_nGc()}, nGb = {hz.get_nGb()}, nC = {hz.get_nC()}')
@contextmanager
def timectx(msgfunc):
"""Context manager to time a block of code."""
start = time()
yield
end = time()
print(msgfunc(end-start))
def print_hzinfo(hz):
print(f'n = {hz.get_n()}, nGc = {hz.get_nGc()}, nGb = {hz.get_nGb()}, nC = {hz.get_nC()}')
Hyperparameters¶
In [17]:
Copied!
AXES = [dict(name='t', bounds=[ 0, 40], step=0.5, unit='s'),
dict(name='x', bounds=[-100, +100], points= 91, unit='m'),
dict(name='v', bounds=[ -20, +20], points= 91, unit='m/s')]
MAX_ACCEL = 1.0 # [mps2]
# HZ2HJ specific
TIME_STEP = AXES[0]['step']
TIME_HORIZON = AXES[0]['bounds'][1]
MIN_BOUNDS = [spec['bounds'][0] for spec in AXES if spec['name'] != 't']
MAX_BOUNDS = [spec['bounds'][1] for spec in AXES if spec['name'] != 't']
GRID_SHAPE = [spec['points'] for spec in AXES if spec['name'] != 't']
AXES = [dict(name='t', bounds=[ 0, 40], step=0.5, unit='s'),
dict(name='x', bounds=[-100, +100], points= 91, unit='m'),
dict(name='v', bounds=[ -20, +20], points= 91, unit='m/s')]
MAX_ACCEL = 1.0 # [mps2]
# HZ2HJ specific
TIME_STEP = AXES[0]['step']
TIME_HORIZON = AXES[0]['bounds'][1]
MIN_BOUNDS = [spec['bounds'][0] for spec in AXES if spec['name'] != 't']
MAX_BOUNDS = [spec['bounds'][1] for spec in AXES if spec['name'] != 't']
GRID_SHAPE = [spec['points'] for spec in AXES if spec['name'] != 't']
Program¶
In [18]:
Copied!
## SPECIFICATION
T = BoundedSet(x=(-50, +50))
phi = ALWAYS(EVENTUALLY(T))
## SPECIFICATION
T = BoundedSet(x=(-50, +50))
phi = ALWAYS(EVENTUALLY(T))
/tmp/ipykernel_2551/3789796335.py:3: DeprecationWarning: BoundedSet is renamed to AlignedBoxSet for clarity, please use AlignedBoxSet instead of BoundedSet. T = BoundedSet(x=(-50, +50))
Two interpretations of ALWAYS¶
In [19]:
Copied!
# Define ALWAYS through RCI set (relating to ¬◇¬ψ)
@primitive(ALWAYS('_1'))
def Always_rci(_1: 'TLTLike') -> tuple[SetBuilder, APPROXDIR]:
b1, a1 = _1._builder, _1._approx
ao = APPROXDIR.UNDER
return (
AppliedSet('rci', b1),
ao + a1 if ao * a1 == APPROXDIR.EXACT else
a1 if ao == a1 else
APPROXDIR.INVALID,
)
# Define Always through fixed-point iteration
@primitive(ALWAYS('_1'))
def Always_fp(_1: 'TLTLike') -> tuple[SetBuilder, APPROXDIR]:
N = int(TIME_HORIZON / TIME_STEP)
phi = 'psi'
for _ in range(N-1):
phi = AND('psi', NEXT(phi))
tree = TLT(phi, where={'psi': _1})
return (tree._builder, tree._approx)
# Define ALWAYS through RCI set (relating to ¬◇¬ψ)
@primitive(ALWAYS('_1'))
def Always_rci(_1: 'TLTLike') -> tuple[SetBuilder, APPROXDIR]:
b1, a1 = _1._builder, _1._approx
ao = APPROXDIR.UNDER
return (
AppliedSet('rci', b1),
ao + a1 if ao * a1 == APPROXDIR.EXACT else
a1 if ao == a1 else
APPROXDIR.INVALID,
)
# Define Always through fixed-point iteration
@primitive(ALWAYS('_1'))
def Always_fp(_1: 'TLTLike') -> tuple[SetBuilder, APPROXDIR]:
N = int(TIME_HORIZON / TIME_STEP)
phi = 'psi'
for _ in range(N-1):
phi = AND('psi', NEXT(phi))
tree = TLT(phi, where={'psi': _1})
return (tree._builder, tree._approx)
HJ Implementation¶
In [20]:
Copied!
## CONSTRUCT TLT
TLT.select(ContLTL | Always_rci)
tree = TLT(phi)
## INITIALIZE IMPLEMENTATION
dynamics = dict(cls=HJDoubleIntegrator,
min_accel=-MAX_ACCEL,
max_accel=+MAX_ACCEL)
impl = TVHJImpl(dynamics, AXES)
## REALIZE ROOT NODE OF TLT
with timectx(lambda t: f"Realization with HJ took {t:.2f} seconds"):
out = tree.realize(impl)
## Plot
impl.plot(
out,
axes=('x', 'v', 't'),
colorscale='greens',
camera_eye=impl.PLOT.EYE_MH_W,
)
## CONSTRUCT TLT
TLT.select(ContLTL | Always_rci)
tree = TLT(phi)
## INITIALIZE IMPLEMENTATION
dynamics = dict(cls=HJDoubleIntegrator,
min_accel=-MAX_ACCEL,
max_accel=+MAX_ACCEL)
impl = TVHJImpl(dynamics, AXES)
## REALIZE ROOT NODE OF TLT
with timectx(lambda t: f"Realization with HJ took {t:.2f} seconds"):
out = tree.realize(impl)
## Plot
impl.plot(
out,
axes=('x', 'v', 't'),
colorscale='greens',
camera_eye=impl.PLOT.EYE_MH_W,
)
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[20], line 8 4 tree = TLT(phi) 5 6 ## INITIALIZE IMPLEMENTATION 7 ----> 8 dynamics = dict(cls=HJDoubleIntegrator, 9 min_accel=-MAX_ACCEL, 10 max_accel=+MAX_ACCEL) 11 NameError: name 'HJDoubleIntegrator' is not defined
HZ Implementation¶
In [21]:
Copied!
## CONSTRUCT TLT
TLT.select(DiscLTL | Always_fp)
tree = TLT(phi)
## INITIALIZE IMPLEMENTATION
dynamics = HZDoubleIntegrator(max_accel=MAX_ACCEL, dt=TIME_STEP)
impl = ZonoOptImpl(dynamics,
AXES[1:], # exclude time axis
zono.interval_2_zono(zono.Box([-MAX_ACCEL], [+MAX_ACCEL])),
TIME_HORIZON,
time_step=TIME_STEP)
## REALIZE ROOT NODE OF TLT
with timectx(lambda t: f"Realization with HZ took {t:.2f} seconds"):
out = tree.realize(impl) # NOTE: Fails due to invalid approximation direction!
print_hzinfo(out)
## Plot
impl.plot_fill(
out,
axes=('x', 'v'),
)
## CONSTRUCT TLT
TLT.select(DiscLTL | Always_fp)
tree = TLT(phi)
## INITIALIZE IMPLEMENTATION
dynamics = HZDoubleIntegrator(max_accel=MAX_ACCEL, dt=TIME_STEP)
impl = ZonoOptImpl(dynamics,
AXES[1:], # exclude time axis
zono.interval_2_zono(zono.Box([-MAX_ACCEL], [+MAX_ACCEL])),
TIME_HORIZON,
time_step=TIME_STEP)
## REALIZE ROOT NODE OF TLT
with timectx(lambda t: f"Realization with HZ took {t:.2f} seconds"):
out = tree.realize(impl) # NOTE: Fails due to invalid approximation direction!
print_hzinfo(out)
## Plot
impl.plot_fill(
out,
axes=('x', 'v'),
)
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[21], line 8 4 tree = TLT(phi) 5 6 ## INITIALIZE IMPLEMENTATION 7 ----> 8 dynamics = HZDoubleIntegrator(max_accel=MAX_ACCEL, dt=TIME_STEP) 9 impl = ZonoOptImpl(dynamics, 10 AXES[1:], # exclude time axis 11 zono.interval_2_zono(zono.Box([-MAX_ACCEL], [+MAX_ACCEL])), NameError: name 'HZDoubleIntegrator' is not defined