Discrete Reachability

Chapter 1: Graph Formulation

A discrete system is a directed graph G = (V, E). The vertices V are the system states. An edge (u, v) ∈ E means the system can transition from state u to state v in one step. Some edges may carry weights (probabilities, costs).

Adjacency Representations

How you store the graph determines the cost of reachability queries:

For reachability analysis, we almost always want both forward and backward edges, so we store two adjacency lists: one for outneighbors, one for inneighbors. This doubles the storage but makes both directions equally fast.

Transition Systems vs. Graphs

In verification, the graph often represents a transition system: a tuple (S, S₀, T) where S is the state set, S₀ ⊆ S is the set of initial states, and T ⊆ S × S is the transition relation. States may have labels (e.g., "safe", "failure"). The reachability question becomes: starting from S₀, can the system reach any state labeled "failure"?

python
# Simple directed graph using adjacency lists
class DiGraph:
    def __init__(self, n):
        self.n = n
        self.out_adj = [[] for _ in range(n)]
        self.in_adj  = [[] for _ in range(n)]

    def add_edge(self, u, v):
        self.out_adj[u].append(v)
        self.in_adj[v].append(u)

Chapter 2: Forward Reachability

The most natural reachability question: starting from a set of initial states S₀, which states can the system eventually reach? This is forward reachability, and it is solved by breadth-first search (BFS).

BFS explores states layer by layer. Layer 0 is S₀ itself. Layer 1 is everything reachable in one step from S₀. Layer 2 is everything reachable in two steps. And so on, until no new states are discovered.

When R_k+1 = R_k, we have found the full forward reachable set. The algorithm terminates in at most |V| steps (since each step adds at least one new state, or we stop).

python
from collections import deque

def forward_reachable(graph, initial_states):
    """BFS from a set of initial states. Returns reachable set."""
    visited = set(initial_states)
    queue = deque(initial_states)
    while queue:
        u = queue.popleft()
        for v in graph.out_adj[u]:
            if v not in visited:
                visited.add(v)
                queue.append(v)
    return visited

Chapter 3: Backward Reachability

Forward reachability answers: "what can the system reach?" Backward reachability answers the dual question: "what can reach the target?" Given a set of target states T (perhaps failure states), backward reachability computes the set of all states from which T is eventually reachable.

This is BFS on the reverse graph — following edges backward (using inneighbors instead of outneighbors). The set B_k contains all states that can reach T in at most k steps.

Combining Forward and Backward

The most powerful analysis combines both directions:

If R ∩ B is empty, the system is safe: no initial condition can lead to the failure state. If R ∩ B is non-empty, every state in the intersection lies on a path from S₀ to T — a potential failure trajectory.

Try the simulation in Chapter 2: switch to "Backward" mode and watch BFS expand from the goal. The backward reachable set at layer k contains all cells that can reach the goal in exactly k steps. Compare the forward and backward wavefronts — they meet where paths connect start to goal.

python
def backward_reachable(graph, target_states):
    """BFS on reverse graph from target states."""
    visited = set(target_states)
    queue = deque(target_states)
    while queue:
        v = queue.popleft()
        for u in graph.in_adj[v]:  # inneighbors!
            if u not in visited:
                visited.add(u)
                queue.append(u)
    return visited

# Safety check
R = forward_reachable(G, S0)
B = backward_reachable(G, T_fail)
if R & B:
    print("UNSAFE: failure is reachable")
else:
    print("SAFE: no path from S0 to failure")

Chapter 4: Counterexample Search

Knowing that a failure state is reachable is only half the story. We also want a counterexample: a concrete path from an initial state to the failure state. This path explains how the system can fail and is essential for debugging.

A counterexample is a sequence of states s₀, s₁, ..., s_k where s₀ ∈ S₀, s_k ∈ T, and each (s_i, s_i+1) ∈ E. This is simply a path in the graph from the initial set to the target set.

Search Strategies

BFS finds the shortest counterexample (fewest transitions). This is valuable because shorter failure paths are often more plausible and easier to understand.

DFS uses less memory (O(depth) vs O(branching factor^depth)) and may find a counterexample faster if one exists deep in the graph. But the path it finds may not be the shortest.

A* is optimal for weighted graphs. If edges have costs (e.g., negative log-probabilities), A* with an admissible heuristic finds the minimum-cost path — the most likely failure trajectory.

Extracting the Path

During BFS or DFS, maintain a parent pointer for each discovered state: parent[v] = u means "v was first reached from u." Once the search reaches a target state t, trace back from t through the parent pointers to reconstruct the full path.

python
def find_counterexample(graph, s0, target):
    """BFS shortest path from s0 to any state in target."""
    parent = {s0: None}
    queue = deque([s0])
    while queue:
        u = queue.popleft()
        if u in target:
            # Trace back
            path = []
            while u is not None:
                path.append(u)
                u = parent[u]
            return path[::-1]
        for v in graph.out_adj[u]:
            if v not in parent:
                parent[v] = u
                queue.append(v)
    return None  # no path exists

Chapter 5: Boolean Satisfiability (SAT)

For small graphs, BFS and DFS work perfectly. But what about systems with 2¹⁰⁰ states? Explicit enumeration is hopeless. SAT-based methods exploit the fact that many large systems have compact symbolic descriptions — the state can be encoded as a vector of Boolean variables, and the transition relation as a logical formula.

The Encoding

Suppose the system has n Boolean state variables x₁, ..., x_n. A state is a truth assignment to these variables (2ⁿ possible states). The transition relation T(x, x') is a Boolean formula relating the current state x to the next state x'.

The k-step reachability question becomes: is there a sequence of states x⁽⁰⁾, x⁽¹⁾, ..., x^(k) such that x⁽⁰⁾ ∈ S₀, x^(k) ∈ T, and T(x⁽ⁱ⁾, x⁽ⁱ⁺¹⁾) holds for all i? This is a satisfiability problem:

This formula, called the bounded model checking (BMC) formula, is satisfiable if and only if there exists a failure path of length k. A SAT solver determines satisfiability and, if satisfiable, returns a satisfying assignment — which is exactly the counterexample.

Incremental Deepening

We don't know the minimum failure path length in advance. So we try k = 0, 1, 2, ... increasing the bound until either the formula is satisfiable (failure found) or we reach a completeness threshold. The completeness threshold is the maximum k such that if no failure exists at depth k, none exists at any depth. For simple reachability, this is the diameter of the graph (longest shortest path).

Chapter 6: Probabilistic Occupancy

So far we have asked binary questions: is a state reachable or not? But if transitions have probabilities, we can ask a richer question: what is the probability of being in each state at time k?

A Markov chain over n states is described by a transition matrix P, where P_ij = Pr(next state = j | current state = i). The occupancy distribution μ_k is a row vector giving the probability of being in each state at time k:

Starting from μ₀ (which puts all probability on the initial states), each multiplication by P propagates the distribution one time step. After k steps, μ_k = μ₀ P^k.

Worked Example

Consider a 3-state system with transition matrix:

Starting from state 0: μ₀ = [1, 0, 0].

μ₁ = [1, 0, 0] · P = [0.5, 0.5, 0.0]

μ₂ = [0.5, 0.5, 0.0] · P = [0.25, 0.40, 0.35]

State 2 is an absorbing state (P₂₂ = 1). Its probability grows monotonically. Eventually, all probability mass concentrates there: μ_∞ = [0, 0, 1].

Chapter 7: Finite & Infinite Horizon

The reachability probability question comes in two flavors, depending on whether we impose a time limit.

Finite Horizon

Finite-horizon reachability: what is the probability of reaching a target set T within k time steps? Define:

This satisfies a recursive formula. If s ∈ T, then p_k(s) = 1 for all k ≥ 0. Otherwise:

Starting from p₀(s) = 1 if s ∈ T, 0 otherwise, we compute p₁, p₂, ... by propagating backward through the transition matrix. After k iterations, p_k(s) gives the exact k-step reachability probability for each state.

Infinite Horizon

Infinite-horizon reachability: what is the probability of ever reaching T? Define:

This is the limit of p_k(s) as k → ∞. If T is reachable with positive probability at each step, the sequence p_k(s) is non-decreasing and bounded by 1, so it converges. In practice, we iterate the recursive formula until convergence (|p_k+1 − p_k| < ε).

Matrix Formulation

Partition the states into target states T and non-target states S \ T. The sub-matrix Q = P restricted to transitions among non-target states governs the dynamics before reaching T. The infinite-horizon reachability probability for non-target states is:

where r_s = ∑_{t ∈ T} P(s, t) is the one-step probability of transitioning into T from state s. The matrix (I − Q)⁻¹ is the fundamental matrix of the absorbing Markov chain — its (i,j) entry counts the expected number of times the process visits state j starting from state i before being absorbed.

Chapter 8: Discrete State Abstractions

What if the system is continuous but we want to use discrete reachability tools? The answer: discretize. Partition the continuous state space into finitely many cells. Build a graph where each cell is a node, and there is an edge from cell A to cell B if the continuous dynamics can take any state in A to some state in B.

Constructing the Abstract Graph

Step 2 is where the nonlinear reachability methods from Chapter 9 enter. For each cell, we use interval arithmetic (or Taylor models) to overapproximate the image under the dynamics. The edges are conservative: we add an edge whenever the overapproximated image might overlap a neighboring cell.

Refinement

If the abstract graph says a failure state is reachable but we suspect this is spurious, we can refine: split cells in the abstract path into smaller cells and re-check. If the path survives refinement, it is likely genuine. If it disappears, it was spurious. This counterexample-guided abstraction refinement (CEGAR) loop is the backbone of many verification tools.

Chapter 9: Resolution Tradeoffs

Every method in this chapter lives on a spectrum: finer resolution = tighter results but more computation. This final chapter maps out the tradeoffs and points forward.

The Resolution Ladder

When to Use What

Small state space (|V| < 10⁶): Use explicit BFS. It is simple, exact, and fast. Store the graph in memory and run forward/backward reachability directly.

Large structured state space (|V| > 10⁶, compact description): Use SAT/BMC. Encode the transitions as Boolean formulas and let the SAT solver handle the search. This is the sweet spot for hardware and protocol verification.

Continuous system, low dimension (n ≤ 5): Use discrete abstraction with adaptive refinement. The cell count is manageable, and the graph tools from this chapter apply directly.

Continuous system, high dimension (n > 10): Return to Chapter 9's nonlinear methods (interval arithmetic, Taylor models) or use sampling-based falsification (Chapters 4-5). Discrete abstraction is impractical due to the curse of dimensionality.

Connections

This chapter closes the reachability trilogy (Chapters 8-10). Each chapter handles a different flavor of dynamics:

The reachability methods provide pre-deployment guarantees: before the system runs, we prove it cannot reach unsafe states. Chapters 11 and 12 address what happens during deployment: monitoring the system in real time and explaining its decisions when things go wrong.