Lucid

Learning-Enabled Uncertainty-Aware Certification of Stochastic Dynamical Systems

E. Casablanca

O. Schön

P. Zuliani

S. Soudjani

Motivation

Reality is not ideal
Noisy sensors
External agents
Unknown/Complex dynamics
Can we trust AI?

Stochastic dynamical systems

How to learn a safety certificate for a system with unknown stochastic dynamics?

Loading diagram...

Loading diagram...

Loading diagram...

Loading diagram...

Unknown: we don't know the model

Unknown + Noisy: we don't know the model & there is random noise

Well, our aim is precisely to provide a way of learning a safety certificate for a system with unknown stochastic dynamics. Trying to abstract the problem a bit, we can imagine the position of the car in terms of x and y coordinates at a given time T. This is the state of the system. At each time step, the controller applies an input U. For example, it can turn the steering wheel, accelerate or brake. The action will, in turn, result in a change of the car's position according to some transition function F, which describes the dynamics of the system. Now, in an ideal world, if we were to model the problem mathematically, we can consider state spaces with an arbitrary number of dimensions. We would also like to have access to a well defined, and possibly well behaved, transition functions. Unfortunately, as we saw, we may not know the dynamics of the system exactly or they may be too complex to model. So, we need a way to reason about them indirectly, treating them as a black-box. Moreover, in the real world, there is randomness or noise involved in every measurement, which adds yet another layer of complexity. -->

Stochastic dynamical systems

We are interfacing with

state space $\mathbb{X} \sube \mathbb{R}^{d_x}$ ,
control input space $\mathcal{U} \sube \mathbb{R}^{d_u}$ ,
a probability measure $P_w$ over $\mathcal{W} \sube \mathbb{R}^{d_w}$ ,
a transition function $f: \mathbb{X} \times \mathcal{U} \times \mathcal{W} \to \mathbb{X}$ .

From the current state $x_t \in \mathbb{X}$ and control input $u_t \in \mathcal{U}$ , the next state $x_{t+1} \in \mathbb{X}$ is obtained as

x_{t+1} = f(x_t, u_t, w_t), \quad w_t \sim P_w

Goal

The goal is to find a lower bound on the probability that the system $\mathcal{S}$ satisfies some specification (safety) over a time horizon $T \in \mathbb{N}$ .

Barrier 3 function

Barrier certificates

Barrier certificates leverage set invariance to provide safety guarantees $^1$ .

[1] A. Salamati et al. (2021) Data-Driven Safety Verification of Stochastic Systems via Barrier Certificates

Barrier certificates

Given a state space $\mathbb{X} \subseteq \mathbb{R}^n$ , initial set $\mathbb{X}_0 \subseteq \mathbb{X}$ and unsafe set $\mathbb{X}_U \subseteq \mathbb{X}$ ,
a barrier certificate is a function $B: \mathbb{X} \to \mathbb{R}$ and constants $\gamma > \eta \ge 0$ , $c \ge 0$ such that

\begin{aligned} & \text{(i) } & & \forall x \in \mathbb{X} & & B(x) \ge 0 \newline & \text{(ii) } & & \forall x_0 \in \mathbb{X}_0 & & B(x_0) \le \eta \newline & \text{(iii) } & & \forall x_u \in \mathbb{X}_u & & B(x_u) \ge \gamma \newline & \text{(iv) } & & \forall x \in \mathbb{X} & & \mathbb{E}[B(X^+)| X = x] - B(x) \le c \end{aligned}

Barrier certificates

For any valid barrier and time horizon $T$ , the safety probability is lower bounded by $^2$

\mathbb{P}(\mathcal{S} \text{ safe}) \ge 1 - \frac{\eta + cT}{\gamma}

Barrier 3 contour

[2] H. J. Kushner (1967) Stochastic Stability and Control

Data-driven Barrier Certificates

Data-driven approach
Does not require knowledge of $f$
Learns by sampling 1-step transitions $N \in \mathbb{N}$ times

\begin{array}{c} &\{x_i , u_i , x_{i}^{+} \}_{i = 1}^{N} \newline\newline &x_i^{+} = f(x_i, u_i, w_i), \quad w_i \sim P_w \end{array}

Barrier constraints enforced via Conditional Mean Embeddings $^3$

Loading diagram...

[3] J. Park, K. Muandet (2020) A Measure-Theoretic Approach to Kernel Conditional Mean Embeddings

Spectral abstraction

$\mathcal{H}$ is infinite-dimensional
Gaussian kernel $\approx$ truncated Fourier series expansion (feature map)
Constraint relaxation via trigonometric bounds on lattices
Can be used in an linear optimization $^4$

Barrier fourier

[4] O. Schön, Z. Zhong, S. Soudjani (2025) Kernel-Based Learning of Safety Barriers

Lucid

Learning-Enabled Uncertainty-Aware Certification of Stochastic Dynamical Systems $^5$

C++ for performance and control
Python3 interface for usability
Modular and extensible
Python package (pylucid)
Docker image (ghcr.io/tendto/lucid)
Live demo on the browser

[5] LUCID GitHub Repository

Features

Tool	Guarantees	Data Driven	Stochastic Dyn.	Non-poly Dyn.	Stat. correct.
LUCID	✓	✓	✓	✓	✓
TRUST (2025)	✓	✓	✗	✗	✓
FOSSIL (2024)	✓	✗	✗	✓	`NA`
OMNISAFE (2024)	✗	✓	✓	✓	✗
NPINTERVAL (2023)	✓	✗	✓/✗	✓	`NA`

Main components

Lucid flow

Execution flow

Lucid flow

Architecture

Loading diagram...

Example

Define the system and sample some points.

N = 1000
# Transition function f
f = lambda x: ...

# State space
X_bounds = RectSet((-3, -2), (2.5, 1))
# Initial set
X_i = MultiSet(
    RectSet((1, -0.5), (2, 0.5)),
    RectSet((-1.8, -0.1), (-1.2, 0.1)),
    RectSet((-1.4, -0.5), (-1.2, 0.1)),
)
# Unsafe set
X_u = MultiSet(
    RectSet((0.4, 0.1), (0.6, 0.5)),
    RectSet((0.4, 0.1), (0.8, 0.3))
)

x_samples = X_bounds.sample(N)
# We simulate the next state samples by applying f
xp_samples = f(x_samples)

Loading diagram...

Example

Train the estimator to predict the feature map.

# Create a truncated Fourier feature map
feature_map = LinearTruncatedFourierFeatureMap(
    num_frequencies=7, # Include the 0th frequency
    sigma_f=1.0,
    sigma_l=[0.05, 0.15],
    X_bounds=X_bounds,
)

# Apply the feature map to the samples
xp_features = feature_map(xp_samples)

# Train the estimator
estimator = KernelRidgeRegressor(
    kernel=GaussianKernel(sigma_l=[0.320, 0.157],
                          sigma_f=1.0),
    regularization_constant=1e-8,
)
estimator.fit(x_samples, xp_features)

Loading diagram...

Example

Synthesize the barrier certificate.

# Synthesize the barrier certificate
b = FourierBarrierCertificate(T=5, gamma=1.0)
success = b.synthesize(
    X_bounds=X_bounds,
    X_init=X_init,
    X_unsafe=X_unsafe,
    estimator=estimator,
    feature_map=feature_map,
    lattice_resolution=330,
    parameters=FourierBarrierCertificateParameters(
        set_scaling=0.025,
    ),
)

# [info] Computing LP coefficients
# [info] Optimizing
# [info] Solution found
# [info] Safety probability is 68.08%

Loading diagram...

Live demo

Future work

Tackle Non-Linear optimization directly.
Temporal logic specifications, not just safety.
Support more kernels, feature maps and estimators.
Extension to control synthesis.
Improved hyperparameter tuning.

References

Thank you!

Barrier certificates

\begin{aligned} & \text{(i) } & & \forall x \in \mathbb{X} & & B(x) \ge 0 \newline & \text{(ii) } & & \forall x_0 \in \mathbb{X}_0 & & B(x_0) \le \eta \newline & \text{(iii) } & & \forall x_u \in \mathbb{X}_u & & B(x_u) \ge \gamma \newline & \text{(iv) } & & \forall x \in \mathbb{X} & & \mathbb{E}[B(X^+)| X = x] - B(x) \le c \end{aligned}

Conditional mean embeddings

k_\mathbb{X}(x,x') := \sigma_f^2 \exp\left( -\frac{1}{2} (x-x')^\top \Sigma^{-1} (x-x') \right),\quad \Sigma:=\text{diag}(\sigma_l)^{2}

\mu_{k_\mathbb{Y}|k_\mathbb{X}}(p)(\cdot) := \mathbb{E}_{p}[\phi_\mathbb{Y}(Y)\mid X=\cdot]

\mathbb{E}_{p}[f(Y)\mid X=x] = \langle f,{\mu}_{k_\mathbb{Y}|k_\mathbb{X}}(p)(x)\rangle_{\mathcal{H}_{k_\mathbb{Y}}}

\hat{\mu}_{k_\mathbb{Y}|k_\mathbb{X}}^N(\cdot) := k_{\hat{X}}(\cdot)^\top\left[ K_{\hat{X}}^N+N\lambda I_N\right]^{-1\!} \phi_\mathbb{Y}(\hat{Y}_N)

\mathbb{E}_{p}[f(Y)\mid X=x] \approx k_{\hat{X}}^N(x)^\top\left[ k_{\hat{X}}^N+N \lambda I_N\right]^{-1\!} f(\hat{Y}_N)

LP formulation

\tiny{ \begin{alignedat}{3} & \min \quad & & \eta + cT, & & \\ & \text{subject to}\quad & & \check{\boldsymbol{B}}_{\tilde{N}}^{\mathbb{X}_0}\leq\phi_M(x_0^{(i)})^\top b\leq\hat{\eta}, \quad & & i=1,\ldots,N_0, \\ & & & \hat{\gamma}\leq\phi_M(x_u^{(i)})^\top b\leq\hat{\boldsymbol{B}}_{\tilde{N}}^{\mathbb{X}_U}, \quad & & i=1,\ldots,N_U, \\ & & & \check{\boldsymbol{B}}_\Delta^{\mathbb{X}}\leq\phi_M(x^{(i)})^\top(Hb - b) \leq \hat{\Delta}, \quad & & i=1,\ldots,N, \\ & & & \hat{\xi}\leq\phi_M(x^{(i)})^\top b\leq\hat{\boldsymbol{B}}_{\tilde{N}}^{{\mathbb{X}}}, \quad & & i=1,\ldots,N, \\ & & & \phi_M(x^{(i)})^\top b\leq\hat{\boldsymbol{B}}_{\tilde{N}}^{\tilde{\mathbb{X}}\setminus\mathbb{X}_0}, \quad & & i=1,\ldots,\tilde{N}-N_0, \\ & & & \check{\boldsymbol{B}}_{\tilde{N}}^{\tilde{\mathbb{X}}\setminus\mathbb{X}_U}\leq\phi_M(x^{(i)})^\top b, \quad & & i=1,\ldots,\tilde{N}-N_U, \\ & & & \check{\boldsymbol{B}}_{\tilde{N}}^{\tilde{\mathbb{X}}\setminus\mathbb{X}}\leq\phi_M(x^{(i)})^\top b, \quad & & i=1,\ldots,\tilde{N}-N, \\ & & & \phi_M(x^{(i)})^\top (Hb - b)\leq\hat{\boldsymbol{B}}_\Delta^{\tilde{\mathbb{X}}\setminus\mathbb{X}}, \quad & & i=1,\ldots,\tilde{N}-N, \\ & & & c\geq 0,\,1>\eta\geq 0,\, b\in\R^{2M+1} & & \end{alignedat} }

Example

Define the system and sample some points.

N = 1000
# Transition function f
f = lambda x: ...

# State space
X_bounds = RectSet((-3, -2), (2.5, 1))
# Initial set
X_i = MultiSet(
    RectSet((1, -0.5), (2, 0.5)),
    RectSet((-1.8, -0.1), (-1.2, 0.1)),
    RectSet((-1.4, -0.5), (-1.2, 0.1)),
)
# Unsafe set
X_u = MultiSet(
    RectSet((0.4, 0.1), (0.6, 0.5)),
    RectSet((0.4, 0.1), (0.8, 0.3))
)

x_samples = X_bounds.sample(N)
# We simulate the next state samples by applying f
xp_samples = f(x_samples)

Loading diagram...

Example

Train the estimator to predict the feature map.

# Create a truncated Fourier feature map
feature_map = LinearTruncatedFourierFeatureMap(
    num_frequencies=7, # Include the 0th frequency
    sigma_f=1.0,
    sigma_l=[0.05, 0.15],
    X_bounds=X_bounds,
)

# Apply the feature map to the samples
xp_features = feature_map(xp_samples)

# Train the estimator
estimator = KernelRidgeRegressor(
    kernel=GaussianKernel(sigma_l=[0.320, 0.157],
                          sigma_f=1.0),
    regularization_constant=1e-8,
)
estimator.fit(x_samples, xp_features)

Loading diagram...

Example

Synthesize the barrier certificate.

# Synthesize the barrier certificate
b = FourierBarrierCertificate(T=5, gamma=1.0)
success = b.synthesize(
    X_bounds=X_bounds,
    X_init=X_init,
    X_unsafe=X_unsafe,
    estimator=estimator,
    feature_map=feature_map,
    lattice_resolution=330,
    parameters=FourierBarrierCertificateParameters(
        set_scaling=0.025,
    ),
)

# [info] Computing LP coefficients
# [info] Optimizing
# [info] Solution found
# [info] Safety probability is 68.08%

Loading diagram...