Reachability of discrete polynomial systems with boxes using Bernstein methods

[mforets]

CC: @roccaa

Introduction

Given a polynomial function f : R^n -> R, Bernstein methods allow to overapproximate the range of f(X) where X is, eg. a box. The method involves computing Bernstein coefficients, either using the implicit form (tensor), or through explicit matrix operations. See BernsteinExpansions.jl for further details.

We would like to write a reachability algorithm for polynomial ODEs, first in discrete time, then in continuous time. This approach has been described in the literature, see for example:

• http://www-verimag.imag.fr/~tdang/Papers/RCJournal2012.pdf
• https://arxiv.org/pdf/1607.02200.pdf
• https://people.eecs.berkeley.edu/~tommasodreossi/papers/hscc2016.pdf

The reachability method for a discrete dynamical system is described as follows. Given x_{n+1} = g(xn) = xn + dt*f(xn) we transform it to set representation, and to compute the successor we determine the enclosure of the image set by computing the Bernstein coeffs of the RHS function where the feasible set is the current set.

Using boxes:

x <= c ?-> search for c s.t x<= max(y) | y \in Xn+1 = {y, y = g(z), z \in Xn}

Given a direction D:

Dx<= max(Dy) | y \in Xn+1 = {y, y = g(z), z \in Xn}

Example

Let: x_{n+1} = x_n^2 - x_n, n >= 0, x_0 \in [0.99, 1.01]. We proceed as follows:

    1. Transform to set representation,  X_0 = [0.99, 1.01]. 2 directions -1 and 1 .
    2. For j = 0, 1, ..., N-1 :
        X_{j+1} = [max( - (x^2 -x) | x\in X_j) ,   max( (x^2 -x) | x\in X_j)  ]

    3. Return X_0, X_1, .., X_N

In continuous time, one has to consider the discretization error, e, and use a bloated initial set X_0' = Box(conv(X_0,X_1) \oplus e' ). In main iteration step, for X_{j+1} we take the M-sum with e.

API

The API would require several changes:

• the option algorithm will play a more important role inside solve.jl
• push the steps Transformation, Sparse/dense matrix conversion further inside reach algorithms for affine systems
• ??

---

Notes

(these are possible features to be analyzed, not for the preliminary implementation)

• (perf) if you have directions +1 and -1, you only need to compute the Bersntein coeffs only once.
• (generalization) this generalizes to n > 1 since directions will be [1,0, 0, ...] etc
• (accuracy) splitting
• (perf) different methods to compute Bernstein coeffs

[schillic]

This sounds like you will observe a wrapping effect. Any way around that?

[roccaa]

Yes wrapping effect is the biggest problem. On this technique there is no magic solution to my best knowledge. We can always use the classical trick: splitting, adaptive time steps, better set enclosure/ dynamic templates, CEGAR procedure etc...

[schillic]

push the steps Transformation, Sparse/dense matrix conversion further inside reach algorithms for affine systems

Yes, what about we introduce a new type for distinguishing different algorithms (like AnalysisMode) and then splitting the solve function into just calling generic functions like this:

function solve(mode::AnalysisMode, system, options)
    init = initialize(mode, system, options)
    res = compute(mode, init, system, options)
    output = project(res, system, options)
end

Then we can write specialized functions by dispatching on mode.

It is not clear to me at the moment if there are any synergies between different modes. The way I wrote it above, a three-liner does not really help much. So another idea is that we directly write a dispatched version of solve for each type of problem. This would be more flexible.